Script:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import os
import glob
import logging
from typing import Tuple, List, Dict, Any
import pandas as pd
import geopandas as gpd
import numpy as np
from tqdm import tqdm
from colorama import init, Fore, Style
from shapely.geometry import LineString, Point, MultiLineString
import matplotlib.pyplot as plt
import folium
from datetime import datetime
from scipy import stats
import math

# Initialize colorama
init(autoreset=True)

# Configure logging
logging.basicConfig(
    filename='geoprocessor.log',
    filemode='a',
    format='%(asctime)s - %(levelname)s - %(message)s',
    level=logging.INFO
)

# Constants
REQUIRED_COLUMNS = ['Time', 'Image', 'X', 'Y', 'Z', 'Roll', 'Pitch', 'Heading']
OUTPUT_FOLDER = 'analysis_output'
MAX_SHAPEFILE_FIELD_LENGTH = 10  # ESRI Shapefile field name limit

def setup_output_folders(base_path: str) -> Dict[str, str]:
    """
    Creates necessary output folders for analysis results.

    Args:
        base_path (str): The base directory for output.

    Returns:
        Dict[str, str]: A dictionary containing paths to specific output folders.
    """
    folders = {
        'plots': os.path.join(base_path, 'plots'),
        'reports': os.path.join(base_path, 'reports'),
        'maps': os.path.join(base_path, 'maps'),
        'stats': os.path.join(base_path, 'stats')
    }
    
    for folder in folders.values():
        os.makedirs(folder, exist_ok=True)
    
    logging.info(f"Output folders set up at {base_path}.")
    return folders

def get_input(prompt: str, default: Any = None, cast_type: type = str) -> Any:
    """
    Prompts the user for input with an optional default value and casts it to the specified type.

    Args:
        prompt (str): The input prompt to display.
        default (Any, optional): The default value if no input is provided. Defaults to None.
        cast_type (type, optional): The type to cast the input to. Defaults to str.

    Returns:
        Any: The user input cast to the specified type.
    """
    while True:
        try:
            if default is not None:
                user_input = input(f"{prompt} [{default}]: ").strip()
                if user_input == "":
                    user_input = default
            else:
                user_input = input(f"{prompt}: ").strip()
            return cast_type(user_input)
        except ValueError:
            print(Fore.RED + f"Invalid input. Please enter a value of type {cast_type.__name__}." + Style.RESET_ALL)
            logging.warning(f"Invalid input for prompt '{prompt}'. Expected type {cast_type.__name__}.")

def merge_files(input_folder: str, file_extension: str = 'txt') -> Tuple[pd.DataFrame, int]:
    """
    Merges multiple CSV or TXT files from a specified folder into a single DataFrame.

    Args:
        input_folder (str): Directory containing the files to merge.
        file_extension (str, optional): Extension of files to merge ('csv' or 'txt'). Defaults to 'txt'.

    Returns:
        Tuple[pd.DataFrame, int]: The merged DataFrame and the count of files merged.

    Raises:
        FileNotFoundError: If no files with the specified extension are found.
        ValueError: If no valid data is found for merging.
    """
    file_pattern = os.path.join(input_folder, f'*.{file_extension}')
    file_list = glob.glob(file_pattern)

    if not file_list:
        raise FileNotFoundError(f"No files with extension '{file_extension}' found in '{input_folder}'.")

    df_list = []
    logging.info(f"Starting to merge {len(file_list)} files with extension '{file_extension}' from '{input_folder}'.")

    for file in tqdm(file_list, desc="Merging files", unit="file"):
        try:
            if file_extension == 'csv':
                df = pd.read_csv(file)
            elif file_extension == 'txt':
                df = pd.read_csv(file, delimiter='\t')
            else:
                raise ValueError("Unsupported file extension. Use 'csv' or 'txt'.")

            # Validate required columns
            missing_cols = [col for col in REQUIRED_COLUMNS if col not in df.columns]
            if missing_cols:
                logging.warning(f"File {file} missing columns: {missing_cols}")
                continue

            df['Source_File'] = os.path.basename(file)
            df_list.append(df)

        except Exception as e:
            logging.error(f"Error reading {file}: {e}")
            print(Fore.RED + f"Error reading {file}: {e}" + Style.RESET_ALL)

    if not df_list:
        raise ValueError("No valid data found for merging.")

    merged_df = pd.concat(df_list, ignore_index=True)
    merged_df.drop_duplicates(inplace=True)

    # Convert Time to datetime if it's not already
    if not pd.api.types.is_datetime64_any_dtype(merged_df['Time']):
        try:
            merged_df['Time'] = pd.to_datetime(merged_df['Time'], unit='s')
            logging.info("Converted 'Time' column to datetime.")
        except Exception as e:
            logging.warning(f"Unable to convert Time to datetime: {e}")

    logging.info(f"Merging completed. Files merged: {len(file_list)}. Rows: {len(merged_df)}.")
    return merged_df, len(file_list)

def dataframe_to_geodataframe(df: pd.DataFrame, lon_column: str, lat_column: str, crs: str = 'EPSG:4326') -> gpd.GeoDataFrame:
    """
    Converts a pandas DataFrame into a GeoDataFrame with Point geometries.

    Args:
        df (pd.DataFrame): Input DataFrame containing coordinate columns.
        lon_column (str): Name of the longitude column.
        lat_column (str): Name of the latitude column.
        crs (str, optional): Coordinate Reference System. Defaults to 'EPSG:4326'.

    Returns:
        gpd.GeoDataFrame: GeoDataFrame with Point geometries.

    Raises:
        KeyError: If coordinate columns are missing.
    """
    try:
        geometry = [Point(xy) for xy in zip(df[lon_column], df[lat_column])]
    except KeyError as e:
        logging.error(f"Missing coordinate column: {e}")
        raise KeyError(f"Missing coordinate column: {e}")

    gdf = gpd.GeoDataFrame(df, geometry=geometry)
    gdf.set_crs(crs, inplace=True)
    logging.info(f"DataFrame converted to GeoDataFrame with CRS '{crs}'. Geometries: {len(gdf)}.")
    return gdf

def assign_path_ids(gdf_sorted: gpd.GeoDataFrame, gap_threshold: float) -> Tuple[gpd.GeoDataFrame, pd.DataFrame]:
    """
    Assigns PathID based on temporal gaps and calculates gaps.

    Args:
        gdf_sorted (gpd.GeoDataFrame): GeoDataFrame sorted by time.
        gap_threshold (float): Threshold in seconds to determine a new PathID.

    Returns:
        Tuple[gpd.GeoDataFrame, pd.DataFrame]: GeoDataFrame with PathID assigned and DataFrame of gaps.
    """
    gdf_sorted = gdf_sorted.copy()

    # Convert 'Time' to numeric if necessary
    if not pd.api.types.is_numeric_dtype(gdf_sorted['Time']):
        if pd.api.types.is_datetime64_any_dtype(gdf_sorted['Time']):
            # Convert datetime to UNIX timestamp in seconds
            gdf_sorted['Time_numeric'] = gdf_sorted['Time'].astype(np.int64) // 10**9
        else:
            # Attempt to convert to numeric, forcing errors to NaN
            gdf_sorted['Time_numeric'] = pd.to_numeric(gdf_sorted['Time'], errors='coerce')
            if gdf_sorted['Time_numeric'].isnull().any():
                logging.warning("Some 'Time' values could not be converted to numeric.")

    # Calculate gaps
    gdf_sorted['Previous_Time'] = gdf_sorted['Time_numeric'].shift(1)
    gdf_sorted['Gap'] = gdf_sorted['Time_numeric'] - gdf_sorted['Previous_Time']

    # Assign PathID
    gdf_sorted['PathID'] = (gdf_sorted['Gap'] > gap_threshold).cumsum() + 1

    # Create gaps dataframe
    gaps_df = gdf_sorted[['Gap']].dropna().reset_index().rename(columns={'index': 'Entry_Index'})

    # Clean up
    gdf_sorted.drop(['Time_numeric'], axis=1, inplace=True)

    logging.info("PathID assigned and gaps calculated.")
    return gdf_sorted, gaps_df

def calculate_path_metrics(gdf: gpd.GeoDataFrame) -> Dict[str, pd.DataFrame]:
    """
    Calculates comprehensive path metrics.

    Args:
        gdf (gpd.GeoDataFrame): GeoDataFrame with PathID assigned.

    Returns:
        Dict[str, pd.DataFrame]: Dictionary containing path metrics and statistics.
    """
    try:
        if 'PathID' not in gdf.columns:
            raise ValueError("Column 'PathID' not found. Run assign_path_ids first.")

        # Identify valid paths (at least 2 points)
        path_counts = gdf.groupby('PathID').size()
        valid_path_ids = path_counts[path_counts >= 2].index.tolist()
        logging.info(f"Found {len(valid_path_ids)} valid paths out of {len(path_counts)} total paths.")

        paths_list = []
        geometries = []

        for path_id in tqdm(valid_path_ids, desc="Processing paths"):
            path_data = gdf[gdf['PathID'] == path_id].copy().sort_values('Time')

            # Extract valid points
            valid_points = [(point.x, point.y) for point in path_data.geometry if isinstance(point, Point)]
            if len(valid_points) < 2:
                logging.warning(f"Insufficient valid points for PathID {path_id}: {len(valid_points)} points")
                continue

            line = LineString(valid_points)
            if not line.is_valid:
                logging.warning(f"Invalid LineString for PathID {path_id}")
                continue

            length_km = line.length / 1000  # Assuming units are meters

            # Calculate metrics
            duration = (path_data['Time'].max() - path_data['Time'].min()).total_seconds()
            avg_speed = length_km / (duration / 3600) if duration > 0 else 0
            points_per_km = len(path_data) / length_km if length_km > 0 else 0
            photos_per_km = path_data['Image'].nunique() / length_km if length_km > 0 else 0

            heading_changes = np.abs(np.diff(path_data['Heading'].values))
            heading_changes = np.where(heading_changes > 180, 360 - heading_changes, heading_changes)
            mean_heading_change = np.mean(heading_changes) if len(heading_changes) > 0 else 0

            roll_std = path_data['Roll'].std()
            pitch_std = path_data['Pitch'].std()
            camera_stability = np.sqrt(roll_std ** 2 + pitch_std ** 2) if not np.isnan(roll_std) and not np.isnan(pitch_std) else 0

            path_dict = {
                'PathID': path_id,
                'length_km': round(length_km, 3),
                'duration_s': round(duration, 2),
                'avg_speed_kmh': round(avg_speed, 2),
                'points_per_km': round(points_per_km, 1),
                'photos_per_km': round(photos_per_km, 1),
                'mean_heading_change': round(mean_heading_change, 2),
                'camera_stability': round(camera_stability, 3),
                'roll_std': round(roll_std, 3) if not np.isnan(roll_std) else 0,
                'pitch_std': round(pitch_std, 3) if not np.isnan(pitch_std) else 0
            }

            paths_list.append(path_dict)
            geometries.append(line)

        if not paths_list:
            raise ValueError("No valid paths could be created.")

        # Create DataFrame and GeoDataFrame
        paths_df = pd.DataFrame(paths_list)
        
        # Shorten column names to fit ESRI Shapefile limits
        paths_df = paths_df.rename(columns={
            'duration_s': 'dur_s',
            'avg_speed_kmh': 'avg_spd',
            'points_per_km': 'pts_per_km',
            'photos_per_km': 'phs_per_km',
            'mean_heading_change': 'mean_hd',
            'camera_stability': 'cam_sta',
            'roll_std': 'roll_std',
            'pitch_std': 'pitch_std'
        })

        paths_gdf = gpd.GeoDataFrame(
            paths_df,
            geometry=geometries,
            crs=gdf.crs
        )
        logging.info(f"Created GeoDataFrame with {len(paths_gdf)} paths.")

        # Calculate statistics
        path_stats = gdf[gdf['PathID'].isin(valid_path_ids)].groupby('PathID').agg(
            count_points=('PathID', 'count'),
            min_time=('Time', 'min'),
            max_time=('Time', 'max'),
            min_z=('Z', 'min'),
            max_z=('Z', 'max'),
            mean_z=('Z', 'mean'),
            std_z=('Z', 'std'),
            count_images=('Image', 'nunique'),
            mean_roll=('Roll', 'mean'),
            std_roll=('Roll', 'std'),
            mean_pitch=('Pitch', 'mean'),
            std_pitch=('Pitch', 'std'),
            mean_heading=('Heading', 'mean'),
            std_heading=('Heading', 'std')
        ).reset_index()

        # Shorten column names in path_stats
        path_stats = path_stats.rename(columns={
            'count_points': 'cnt_pts',
            'min_time': 'min_tm',
            'max_time': 'max_tm',
            'min_z': 'min_z',
            'max_z': 'max_z',
            'mean_z': 'mean_z',
            'std_z': 'std_z',
            'count_images': 'cnt_imgs',
            'mean_roll': 'mean_rl',
            'std_roll': 'std_rl',
            'mean_pitch': 'mean_pt',
            'std_pitch': 'std_pt',
            'mean_heading': 'mean_hd',
            'std_heading': 'std_hd'
        })

        logging.info("Path metrics and statistics calculated.")
        return {
            'paths': paths_gdf,
            'path_stats': path_stats
        }

def analyze_temporal_patterns(gdf: gpd.GeoDataFrame) -> Dict[str, pd.DataFrame]:
    """
    Analyzes temporal patterns in the data.

    Args:
        gdf (gpd.GeoDataFrame): GeoDataFrame with temporal data.

    Returns:
        Dict[str, pd.DataFrame]: Dictionary containing temporal metrics.
    """
    temporal_metrics = {}
    
    # Calculate time gaps between consecutive points
    gdf_sorted = gdf.sort_values('Time').copy()
    gdf_sorted['time_gap'] = gdf_sorted['Time'].diff().dt.total_seconds()
    
    # Time gaps statistics
    gap_stats = {
        'min_gap': round(gdf_sorted['time_gap'].min(), 3),
        'max_gap': round(gdf_sorted['time_gap'].max(), 3),
        'mean_gap': round(gdf_sorted['time_gap'].mean(), 3),
        'median_gap': round(gdf_sorted['time_gap'].median(), 3),
        'std_gap': round(gdf_sorted['time_gap'].std(), 3)
    }
    temporal_metrics['gap_stats'] = pd.DataFrame([gap_stats])
    
    # Temporal coverage
    temporal_coverage = {
        'start_time': gdf['Time'].min(),
        'end_time': gdf['Time'].max(),
        'total_duration_hours': round((gdf['Time'].max() - gdf['Time'].min()).total_seconds() / 3600, 2),
        'active_time_hours': round(gdf_sorted['time_gap'][gdf_sorted['time_gap'] < 300].sum() / 3600, 2)  # Excluding gaps > 5 min
    }
    temporal_metrics['coverage'] = pd.DataFrame([temporal_coverage])
    
    # Time distribution by hour
    gdf['hour'] = gdf['Time'].dt.hour
    time_distribution = gdf.groupby('hour').size().reset_index(name='count')
    temporal_metrics['time_distribution'] = time_distribution
    
    logging.info("Temporal analysis completed.")
    return temporal_metrics

def analyze_spatial_patterns(gdf: gpd.GeoDataFrame) -> Dict[str, pd.DataFrame]:
    """
    Analyzes spatial patterns and coverage.

    Args:
        gdf (gpd.GeoDataFrame): GeoDataFrame with spatial data.

    Returns:
        Dict[str, pd.DataFrame]: Dictionary containing spatial metrics.
    """
    spatial_metrics = {}
    try:
        # Create paths as LineStrings with validation
        paths_list = []
        for path_id in gdf['PathID'].unique():
            path_points = gdf[gdf['PathID'] == path_id].geometry.tolist()
            if len(path_points) >= 2:
                coords = [(p.x, p.y) for p in path_points if isinstance(p, Point)]
                path_line = LineString(coords)
                if path_line.is_valid:
                    paths_list.append(path_line)

        if not paths_list:
            raise ValueError("No valid paths could be created for spatial analysis.")

        paths_gdf = gpd.GeoDataFrame(geometry=paths_list, crs=gdf.crs)

        # Calculate total area covered (convex hull)
        try:
            union = paths_gdf.union_all()  # Updated to union_all()
            if union:
                hull = union.convex_hull
                spatial_metrics['total_area'] = pd.DataFrame([{
                    'area_km2': round(hull.area / 1e6, 3),  # Assuming units are meters
                    'perimeter_km': round(hull.length / 1000, 3)
                }])
            else:
                spatial_metrics['total_area'] = pd.DataFrame([{'area_km2': 0, 'perimeter_km': 0}])
            logging.info("Total area and perimeter calculated.")
        except Exception as e:
            logging.error(f"Error calculating convex hull: {e}")
            spatial_metrics['total_area'] = pd.DataFrame([{'area_km2': 0, 'perimeter_km': 0}])

        # Analyze path intersections
        intersections = []
        for i, path1 in enumerate(paths_list):
            for j, path2 in enumerate(paths_list[i + 1:], i + 1):
                if path1.intersects(path2):
                    intersection = path1.intersection(path2)
                    if not intersection.is_empty:
                        # Count intersection points
                        if intersection.geom_type == 'Point':
                            num_points = 1
                        elif intersection.geom_type == 'MultiPoint':
                            num_points = len(intersection)
                        else:
                            num_points = 0  # For LineString intersections, not counting
                        intersections.append({
                            'path1': i + 1,
                            'path2': j + 1,
                            'intersection_points': num_points
                        })

        spatial_metrics['intersections'] = pd.DataFrame(intersections) if intersections else pd.DataFrame()

        # Calculate path density
        total_length = paths_gdf.length.sum() / 1000  # km
        if 'hull' in locals() and hull.area > 0:
            area_km2 = hull.area / 1e6
            path_density = round(total_length / area_km2, 3)  # km/km²
        else:
            path_density = 0
        spatial_metrics['density'] = pd.DataFrame([{
            'total_path_length_km': round(total_length, 3),
            'path_density_km_per_km2': path_density
        }])

        logging.info("Spatial analysis completed.")
        return spatial_metrics

    except Exception as e:
        logging.error(f"Error in spatial analysis: {e}")
        # Return empty metrics rather than failing
        return {
            'total_area': pd.DataFrame(),
            'intersections': pd.DataFrame(),
            'density': pd.DataFrame()
        }

def calculate_overlap_analysis(gdf: gpd.GeoDataFrame, distance_threshold: float = 5.0) -> Dict[str, pd.DataFrame]:
    """
    Analyzes path overlaps and parallel runs.

    Args:
        gdf (gpd.GeoDataFrame): GeoDataFrame with spatial data.
        distance_threshold (float, optional): Distance threshold in meters for overlap. Defaults to 5.0.

    Returns:
        Dict[str, pd.DataFrame]: Dictionary containing overlap metrics.
    """
    overlap_metrics = {}
    try:
        # Create paths with validation
        paths_list = []
        for path_id in gdf['PathID'].unique():
            path_points = gdf[gdf['PathID'] == path_id].geometry.tolist()
            if len(path_points) >= 2:
                coords = [(p.x, p.y) for p in path_points if isinstance(p, Point)]
                path_line = LineString(coords)
                if path_line.is_valid:
                    paths_list.append(path_line)

        if not paths_list:
            logging.warning("No valid paths for overlap analysis.")
            return {
                'overlaps': pd.DataFrame(),
                'summary': pd.DataFrame([{
                    'total_overlapping_paths': 0,
                    'total_overlap_length_m': 0,
                    'mean_overlap_percentage': 0,
                    'max_overlap_percentage': 0
                }])
            }

        paths_gdf = gpd.GeoDataFrame(geometry=paths_list, crs=gdf.crs)
        buffered_paths = paths_gdf.geometry.buffer(distance_threshold)

        overlaps = []
        for i, path1 in enumerate(buffered_paths):
            for j, path2 in enumerate(buffered_paths[i + 1:], i + 1):
                if path1.intersects(path2):
                    intersection = path1.intersection(path2)
                    if not intersection.is_empty:
                        if intersection.geom_type in ['Point', 'MultiPoint']:
                            overlap_length = 0  # Points do not have length
                        elif intersection.geom_type in ['LineString', 'MultiLineString']:
                            overlap_length = intersection.length
                        else:
                            overlap_length = 0
                        overlap_percentage = (overlap_length / min(paths_gdf.geometry[i].length, paths_gdf.geometry[j].length)) * 100 if min(paths_gdf.geometry[i].length, paths_gdf.geometry[j].length) > 0 else 0
                        overlaps.append({
                            'path1': i + 1,
                            'path2': j + 1,
                            'overlap_length_m': round(overlap_length, 2),
                            'overlap_percentage': round(overlap_percentage, 2)
                        })

        overlap_metrics['overlaps'] = pd.DataFrame(overlaps) if overlaps else pd.DataFrame()

        if not overlaps:
            summary = {
                'total_overlapping_paths': 0,
                'total_overlap_length_m': 0,
                'mean_overlap_percentage': 0,
                'max_overlap_percentage': 0
            }
        else:
            overlapping_paths = set([x['path1'] for x in overlaps] + [x['path2'] for x in overlaps])
            summary = {
                'total_overlapping_paths': len(overlapping_paths),
                'total_overlap_length_m': round(sum(x['overlap_length_m'] for x in overlaps), 2),
                'mean_overlap_percentage': round(np.mean([x['overlap_percentage'] for x in overlaps]), 2),
                'max_overlap_percentage': round(max([x['overlap_percentage'] for x in overlaps]), 2)
            }

        overlap_metrics['summary'] = pd.DataFrame([summary])

        logging.info("Overlap analysis completed.")
        return overlap_metrics

    except Exception as e:
        logging.error(f"Error in overlap analysis: {e}")
        return {
            'overlaps': pd.DataFrame(),
            'summary': pd.DataFrame([{
                'total_overlapping_paths': 0,
                'total_overlap_length_m': 0,
                'mean_overlap_percentage': 0,
                'max_overlap_percentage': 0
            }])
        }

def analyze_elevation_profile(gdf: gpd.GeoDataFrame) -> Dict[str, pd.DataFrame]:
    """
    Analyzes elevation profiles and terrain characteristics.

    Args:
        gdf (gpd.GeoDataFrame): GeoDataFrame with elevation data.

    Returns:
        Dict[str, pd.DataFrame]: Dictionary containing elevation metrics.
    """
    elevation_metrics = {}
    
    def calculate_slope(points: List[Point]) -> np.ndarray:
        """Calculate slope between consecutive points."""
        distances = [p1.distance(p2) for p1, p2 in zip(points[:-1], points[1:])]
        elevation_changes = np.diff([p.z for p in points])
        slopes = np.degrees(np.arctan2(elevation_changes, distances))
        return slopes

    # Calculate per-path elevation metrics
    elevation_stats = []
    for path_id in tqdm(gdf['PathID'].unique(), desc="Calculating elevation profiles"):
        path_data = gdf[gdf['PathID'] == path_id].copy().sort_values('Time')

        # Extract 3D points
        points = [Point(row['X'], row['Y'], row['Z']) for _, row in path_data.iterrows() if isinstance(row['X'], (int, float))]

        if len(points) > 1:
            # Calculate slopes
            slopes = calculate_slope(points)

            # Calculate elevation gain/loss
            elevation_changes = np.diff([p.z for p in points])
            elevation_gain = np.sum(elevation_changes[elevation_changes > 0])
            elevation_loss = abs(np.sum(elevation_changes[elevation_changes < 0]))

            elevation_stats.append({
                'PathID': path_id,
                'min_elevation': round(min(p.z for p in points), 2),
                'max_elevation': round(max(p.z for p in points), 2),
                'elevation_range': round(max(p.z for p in points) - min(p.z for p in points), 2),
                'mean_elevation': round(np.mean([p.z for p in points]), 2),
                'elevation_gain_m': round(elevation_gain, 2),
                'elevation_loss_m': round(elevation_loss, 2),
                'mean_slope_deg': round(np.mean(slopes), 2) if len(slopes) > 0 else 0,
                'max_slope_deg': round(np.max(np.abs(slopes)), 2) if len(slopes) > 0 else 0,
                'slope_std_deg': round(np.std(slopes), 2) if len(slopes) > 0 else 0
            })

    elevation_metrics['path_elevation'] = pd.DataFrame(elevation_stats)
    
    # Calculate overall elevation statistics
    if not elevation_metrics['path_elevation'].empty:
        elevation_metrics['summary'] = pd.DataFrame([{
            'total_elevation_gain_m': round(elevation_metrics['path_elevation']['elevation_gain_m'].sum(), 2),
            'total_elevation_loss_m': round(elevation_metrics['path_elevation']['elevation_loss_m'].sum(), 2),
            'max_elevation_m': round(elevation_metrics['path_elevation']['max_elevation'].max(), 2),
            'min_elevation_m': round(elevation_metrics['path_elevation']['min_elevation'].min(), 2),
            'mean_elevation_m': round(elevation_metrics['path_elevation']['mean_elevation'].mean(), 2),
            'max_slope_deg': round(elevation_metrics['path_elevation']['max_slope_deg'].max(), 2),
            'mean_slope_deg': round(elevation_metrics['path_elevation']['mean_slope_deg'].mean(), 2)
        }])
    else:
        elevation_metrics['summary'] = pd.DataFrame([{
            'total_elevation_gain_m': 0,
            'total_elevation_loss_m': 0,
            'max_elevation_m': 0,
            'min_elevation_m': 0,
            'mean_elevation_m': 0,
            'max_slope_deg': 0,
            'mean_slope_deg': 0
        }])

    logging.info("Elevation analysis completed.")
    return elevation_metrics

def analyze_photo_coverage(gdf: gpd.GeoDataFrame) -> Dict[str, pd.DataFrame]:
    """
    Analyzes photo distribution and coverage.

    Args:
        gdf (gpd.GeoDataFrame): GeoDataFrame with photo data.

    Returns:
        Dict[str, pd.DataFrame]: Dictionary containing photo coverage metrics.
    """
    photo_metrics = {}
    
    # Analyze photo spacing
    photo_stats = []
    for path_id in tqdm(gdf['PathID'].unique(), desc="Analyzing photo coverage"):
        path_data = gdf[gdf['PathID'] == path_id].copy().sort_values('Time')
        
        if len(path_data) > 1:
            # Calculate distances between consecutive photos
            points = path_data.geometry.tolist()
            distances = [p1.distance(p2) for p1, p2 in zip(points[:-1], points[1:])]
            
            # Extract photo file information
            photos = path_data['Image'].tolist()
            photo_folders = [os.path.dirname(p) for p in photos]
            
            photo_stats.append({
                'PathID': path_id,
                'num_photos': len(photos),
                'mean_spacing_m': round(np.mean(distances), 2) if distances else 0,
                'max_spacing_m': round(max(distances), 2) if distances else 0,
                'min_spacing_m': round(min(distances), 2) if distances else 0,
                'spacing_std_m': round(np.std(distances), 2) if distances else 0,
                'unique_folders': len(set(photo_folders)),
                'first_photo': photos[0],
                'last_photo': photos[-1]
            })

    photo_metrics['path_photos'] = pd.DataFrame(photo_stats)
    
    # Calculate overall photo statistics
    total_photos = gdf['Image'].nunique()
    total_folders = gdf['Image'].apply(os.path.dirname).nunique()
    
    # Calculate geographical bounds in meters for photo density
    bounds = gdf.total_bounds  # [minx, miny, maxx, maxy]
    if gdf.crs.is_geographic:
        # Reproject to a metric CRS for accurate area calculation
        gdf_projected = gdf.to_crs(epsg=3857)  # Web Mercator
        area_m2 = (gdf_projected.total_bounds[2] - gdf_projected.total_bounds[0]) * (gdf_projected.total_bounds[3] - gdf_projected.total_bounds[1])
    else:
        area_m2 = (bounds[2] - bounds[0]) * (bounds[3] - bounds[1])

    photo_density = round(total_photos / area_m2, 6) if area_m2 > 0 else 0  # photos per m²

    photo_metrics['summary'] = pd.DataFrame([{
        'total_photos': total_photos,
        'total_folders': total_folders,
        'mean_photos_per_path': round(photo_metrics['path_photos']['num_photos'].mean(), 1) if not photo_metrics['path_photos'].empty else 0,
        'mean_photo_spacing_m': round(photo_metrics['path_photos']['mean_spacing_m'].mean(), 2) if not photo_metrics['path_photos'].empty else 0,
        'overall_photo_density_photos_per_m2': photo_density
    }])

    logging.info("Photo coverage analysis completed.")
    return photo_metrics

def analyze_camera_orientation(gdf: gpd.GeoDataFrame) -> Dict[str, pd.DataFrame]:
    """
    Analyzes camera orientation and stability.

    Args:
        gdf (gpd.GeoDataFrame): GeoDataFrame with camera orientation data.

    Returns:
        Dict[str, pd.DataFrame]: Dictionary containing camera orientation metrics.
    """
    camera_metrics = {}
    
    # Calculate per-path camera metrics
    camera_stats = []
    for path_id in tqdm(gdf['PathID'].unique(), desc="Analyzing camera orientation"):
        path_data = gdf[gdf['PathID'] == path_id].copy().sort_values('Time')
        
        # Calculate heading changes
        heading_changes = np.abs(np.diff(path_data['Heading'].values))
        heading_changes = np.where(heading_changes > 180, 360 - heading_changes, heading_changes)
        
        # Calculate stability metrics
        roll_std = path_data['Roll'].std()
        pitch_std = path_data['Pitch'].std()
        stability_index = np.sqrt(roll_std ** 2 + pitch_std ** 2) if not np.isnan(roll_std) and not np.isnan(pitch_std) else 0
        
        camera_stats.append({
            'PathID': path_id,
            'mean_roll': round(path_data['Roll'].mean(), 3) if not path_data['Roll'].isnull().all() else 0,
            'roll_std': round(roll_std, 3) if not np.isnan(roll_std) else 0,
            'mean_pitch': round(path_data['Pitch'].mean(), 3) if not path_data['Pitch'].isnull().all() else 0,
            'pitch_std': round(pitch_std, 3) if not np.isnan(pitch_std) else 0,
            'mean_heading': round(path_data['Heading'].mean(), 2) if not path_data['Heading'].isnull().all() else 0,
            'heading_std': round(path_data['Heading'].std(), 2) if not path_data['Heading'].isnull().all() else 0,
            'max_heading_change': round(max(heading_changes), 2) if len(heading_changes) > 0 else 0,
            'mean_heading_change': round(np.mean(heading_changes), 2) if len(heading_changes) > 0 else 0,
            'sharp_turns': sum(heading_changes > 45) if len(heading_changes) > 0 else 0,
            'stability_index': round(stability_index, 3)
        })

    camera_metrics['path_camera'] = pd.DataFrame(camera_stats)
    
    # Calculate overall camera stability metrics
    if not camera_metrics['path_camera'].empty:
        camera_metrics['summary'] = pd.DataFrame([{
            'overall_roll_std': round(gdf['Roll'].std(), 3) if not gdf['Roll'].isnull().all() else 0,
            'overall_pitch_std': round(gdf['Pitch'].std(), 3) if not gdf['Pitch'].isnull().all() else 0,
            'total_sharp_turns': camera_metrics['path_camera']['sharp_turns'].sum(),
            'mean_stability_index': round(camera_metrics['path_camera']['stability_index'].mean(), 3) if not camera_metrics['path_camera'].empty else 0,
            'worst_stability_index': round(camera_metrics['path_camera']['stability_index'].max(), 3) if not camera_metrics['path_camera'].empty else 0
        }])
    else:
        camera_metrics['summary'] = pd.DataFrame([{
            'overall_roll_std': 0,
            'overall_pitch_std': 0,
            'total_sharp_turns': 0,
            'mean_stability_index': 0,
            'worst_stability_index': 0
        }])

    logging.info("Camera orientation analysis completed.")
    return camera_metrics

def calculate_path_complexity(gdf: gpd.GeoDataFrame) -> Dict[str, pd.DataFrame]:
    """
    Calculates path complexity metrics including curvature and turn analysis.

    Args:
        gdf (gpd.GeoDataFrame): GeoDataFrame with path data.

    Returns:
        Dict[str, pd.DataFrame]: Dictionary containing path complexity metrics.
    """
    complexity_metrics = {}
    
    path_complexity = []
    for path_id in tqdm(gdf['PathID'].unique(), desc="Calculating path complexity"):
        path_data = gdf[gdf['PathID'] == path_id].copy().sort_values('Time')
        
        if len(path_data) > 2:
            # Calculate curvature based on heading changes
            headings = path_data['Heading'].values
            heading_changes = np.abs(np.diff(headings))
            heading_changes = np.where(heading_changes > 180, 360 - heading_changes, heading_changes)

            # Calculate distances between points
            points = path_data.geometry.tolist()
            distances = [p1.distance(p2) for p1, p2 in zip(points[:-1], points[1:])]
            total_distance = sum(distances) if distances else 0

            # Path sinuosity
            line = LineString([(p.x, p.y) for p in points if isinstance(p, Point)])
            direct_distance = Point(points[0].x, points[0].y).distance(Point(points[-1].x, points[-1].y))
            path_sinuosity = round(line.length / direct_distance, 3) if direct_distance > 0 else 1

            # Handle division by zero for turn density
            turn_density = round(sum(heading_changes > 10) / total_distance, 3) if total_distance > 0 else 0

            # Combined curvature index
            curvature_index = round(np.mean(heading_changes) * np.std(heading_changes), 3) if len(heading_changes) > 0 else 0

            path_complexity.append({
                'PathID': path_id,
                'total_turns': sum(heading_changes > 10),
                'sharp_turns': sum(heading_changes > 45),
                'mean_turn_angle': round(np.mean(heading_changes), 2) if len(heading_changes) > 0 else 0,
                'max_turn_angle': round(np.max(heading_changes), 2) if len(heading_changes) > 0 else 0,
                'turn_density_turns_per_meter': turn_density,
                'curvature_index': curvature_index,
                'path_sinuosity': path_sinuosity
            })

    complexity_metrics['path_complexity'] = pd.DataFrame(path_complexity)
    
    # Calculate summary statistics
    if not complexity_metrics['path_complexity'].empty:
        complexity_metrics['summary'] = pd.DataFrame([{
            'total_sharp_turns': complexity_metrics['path_complexity']['sharp_turns'].sum(),
            'mean_sinuosity': round(complexity_metrics['path_complexity']['path_sinuosity'].mean(), 3),
            'max_sinuosity': round(complexity_metrics['path_complexity']['path_sinuosity'].max(), 3),
            'mean_curvature_index': round(complexity_metrics['path_complexity']['curvature_index'].mean(), 3)
        }])
    else:
        complexity_metrics['summary'] = pd.DataFrame([{
            'total_sharp_turns': 0,
            'mean_sinuosity': 0,
            'max_sinuosity': 0,
            'mean_curvature_index': 0
        }])

    logging.info("Path complexity analysis completed.")
    return complexity_metrics

def create_visualizations(gdf: gpd.GeoDataFrame, metrics: Dict[str, Dict[str, pd.DataFrame]], output_folder: str):
    """
    Creates comprehensive visualizations of the analysis results.

    Args:
        gdf (gpd.GeoDataFrame): GeoDataFrame with path data.
        metrics (Dict[str, Dict[str, pd.DataFrame]]): Dictionary containing all calculated metrics.
        output_folder (str): Base output folder path.
    """
    try:
        plots_folder = os.path.join(output_folder, 'plots')
        os.makedirs(plots_folder, exist_ok=True)

        # Set style for all plots
        plt.style.use('default')

        def save_plot(fig, name):
            try:
                plt.tight_layout()
                fig.savefig(os.path.join(plots_folder, f"{name}.png"), dpi=300, bbox_inches='tight')
                plt.close(fig)
            except Exception as e:
                logging.error(f"Error saving plot {name}: {e}")
                print(f"Error saving plot {name}: {e}")

        # Verify metrics structure
        if not isinstance(metrics, dict):
            raise ValueError("Metrics must be a dictionary.")

        # Extract paths data
        paths_data = metrics.get('paths', {}).get('paths', None)
        if not isinstance(paths_data, gpd.GeoDataFrame):
            raise ValueError("Paths metrics not found or invalid.")

        # 1. Path Length Distribution
        try:
            fig, ax = plt.subplots(figsize=(10, 6))
            paths_data['length_km'].hist(bins=30, ax=ax, color='skyblue')
            ax.set_title('Path Length Distribution')
            ax.set_xlabel('Length (km)')
            ax.set_ylabel('Frequency')
            save_plot(fig, 'path_length_distribution')
        except Exception as e:
            logging.error(f"Error creating path length distribution: {e}")
            print(f"Error creating path length distribution: {e}")

        # 2. Elevation Profile per Path
        try:
            fig, ax = plt.subplots(figsize=(12, 6))
            paths_to_plot = min(5, paths_data.shape[0])
            for _, row in paths_data.head(paths_to_plot).iterrows():
                path_id = row['PathID']
                path = gdf[gdf['PathID'] == path_id].sort_values('Time')
                distances = [0]
                for i in range(1, len(path)):
                    distances.append(distances[-1] + path.iloc[i-1].geometry.distance(path.iloc[i].geometry))
                ax.plot(distances, path['Z'].values, label=f'Path {path_id}')
            ax.set_title('Elevation Profiles')
            ax.set_xlabel('Distance (m)')
            ax.set_ylabel('Elevation (m)')
            ax.legend()
            save_plot(fig, 'elevation_profiles')
        except Exception as e:
            logging.error(f"Error creating elevation profiles: {e}")
            print(f"Error creating elevation profiles: {e}")

        # 3. Camera Stability Plot
        try:
            if 'camera' in metrics and 'path_camera' in metrics['camera']:
                fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))
                camera_data = metrics['camera']['path_camera']
                ax1.scatter(camera_data['PathID'], camera_data['roll_std'], label='Roll STD', alpha=0.7)
                ax1.scatter(camera_data['PathID'], camera_data['pitch_std'], label='Pitch STD', alpha=0.7)
                ax1.set_title('Camera Stability per Path')
                ax1.set_xlabel('PathID')
                ax1.set_ylabel('Standard Deviation (degrees)')
                ax1.legend()

                ax2.scatter(camera_data['PathID'], camera_data['stability_index'], color='green', alpha=0.7)
                ax2.set_title('Overall Stability Index per Path')
                ax2.set_xlabel('PathID')
                ax2.set_ylabel('Stability Index')
                save_plot(fig, 'camera_stability')
        except Exception as e:
            logging.error(f"Error creating camera stability plot: {e}")
            print(f"Error creating camera stability plot: {e}")

        # 4. Photo Spacing Analysis
        try:
            if 'photos' in metrics and 'path_photos' in metrics['photos']:
                fig, ax = plt.subplots(figsize=(10, 6))
                photo_data = metrics['photos']['path_photos']
                photo_data['mean_spacing_m'].plot.box(ax=ax, color='lightgreen')
                ax.set_title('Photo Spacing Distribution')
                ax.set_ylabel('Mean Spacing (m)')
                save_plot(fig, 'photo_spacing')
        except Exception as e:
            logging.error(f"Error creating photo spacing plot: {e}")
            print(f"Error creating photo spacing plot: {e}")

        # 5. Speed Profile
        try:
            fig, ax = plt.subplots(figsize=(12, 6))
            paths_data.plot.scatter(x='PathID', y='avg_spd', ax=ax, color='coral')
            ax.set_title('Average Speed per Path')
            ax.set_xlabel('PathID')
            ax.set_ylabel('Average Speed (km/h)')
            save_plot(fig, 'speed_profile')
        except Exception as e:
            logging.error(f"Error creating speed profile: {e}")
            print(f"Error creating speed profile: {e}")

        # 6. Interactive Map
        try:
            gdf_wgs84 = gdf.to_crs('EPSG:4326')
            paths_wgs84 = metrics.get('paths', {}).get('paths', None).to_crs('EPSG:4326') if metrics.get('paths', {}).get('paths', None) is not None else None

            if paths_wgs84 is None or paths_wgs84.empty:
                raise ValueError("No path data available for creating the interactive map.")

            center_lat = gdf_wgs84.geometry.centroid.y.mean()
            center_lon = gdf_wgs84.geometry.centroid.x.mean()

            m = folium.Map(location=[center_lat, center_lon], zoom_start=14)

            for _, path in paths_wgs84.iterrows():
                coordinates = [(p.y, p.x) for p in path.geometry.coords]
                folium.PolyLine(
                    coordinates,
                    weight=2,
                    color=f'#{np.random.randint(0, 16777215):06x}',
                    popup=f'Path {path.PathID}<br>Length: {path.length_km:.2f} km'
                ).add_to(m)

            # Highlight overlapping areas if available
            if 'overlaps' in metrics and 'overlaps' in metrics['overlaps'] and not metrics['overlaps']['overlaps'].empty:
                for _, overlap in metrics['overlaps']['overlaps'].iterrows():
                    try:
                        path1 = paths_wgs84.iloc[overlap['path1'] - 1].geometry
                        path2 = paths_wgs84.iloc[overlap['path2'] - 1].geometry
                        intersection = path1.intersection(path2)
                        if not intersection.is_empty:
                            if intersection.geom_type == 'Point':
                                folium.Marker([intersection.y, intersection.x], icon=folium.Icon(color='red')).add_to(m)
                            elif intersection.geom_type == 'MultiPoint':
                                for point in intersection:
                                    folium.Marker([point.y, point.x], icon=folium.Icon(color='red')).add_to(m)
                            elif intersection.geom_type in ['LineString', 'MultiLineString']:
                                folium.PolyLine(
                                    [(p.y, p.x) for p in intersection.coords],
                                    color='red',
                                    weight=4,
                                    opacity=0.7
                                ).add_to(m)
                    except IndexError:
                        logging.warning(f"Invalid PathID in overlaps: {overlap['path1']} or {overlap['path2']}")
                        continue

            m.save(os.path.join(plots_folder, 'interactive_map.html'))
        except Exception as e:
            logging.error(f"Error creating interactive map: {e}")
            print(f"Error creating interactive map: {e}")

        # 7. Path Complexity Analysis
        try:
            if 'complexity' in metrics and 'path_complexity' in metrics['complexity']:
                fig, ax = plt.subplots(figsize=(12, 6))
                complexity_data = metrics['complexity']['path_complexity']
                scatter = ax.scatter(
                    complexity_data['path_sinuosity'],
                    complexity_data['curvature_index'],
                    c=complexity_data['sharp_turns'],
                    cmap='viridis',
                    alpha=0.7
                )
                cbar = plt.colorbar(scatter, ax=ax)
                cbar.set_label('Number of Sharp Turns')
                ax.set_title('Path Complexity Analysis')
                ax.set_xlabel('Sinuosity')
                ax.set_ylabel('Curvature Index')
                save_plot(fig, 'path_complexity')
        except Exception as e:
            logging.error(f"Error creating path complexity plot: {e}")
            print(f"Error creating path complexity plot: {e}")

        # 8. Time Gap Analysis
        try:
            fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))
            time_gaps = gdf.sort_values('Time')['time_gap']
            filtered_gaps = time_gaps[time_gaps < np.percentile(time_gaps.dropna(), 95)]

            filtered_gaps.hist(bins=50, ax=ax1, color='lightblue')
            ax1.set_title('Temporal Gap Distribution (95th Percentile)')
            ax1.set_xlabel('Gap (seconds)')
            ax1.set_ylabel('Frequency')

            ax2.boxplot(filtered_gaps, vert=False)
            ax2.set_title('Temporal Gap Box Plot')
            ax2.set_xlabel('Gap (seconds)')
            save_plot(fig, 'time_gaps')
        except Exception as e:
            logging.error(f"Error creating time gap analysis plots: {e}")
            print(f"Error creating time gap analysis plots: {e}")

        # 9. Dashboard (Example Placeholder)
        try:
            fig = plt.figure(figsize=(20, 20))
            gs = fig.add_gridspec(4, 4)

            # Path Length Distribution on Dashboard
            ax1 = fig.add_subplot(gs[0, :2])
            paths_data['length_km'].hist(bins=30, ax=ax1, color='skyblue')
            ax1.set_title('Path Length Distribution')
            ax1.set_xlabel('Length (km)')
            ax1.set_ylabel('Frequency')

            # Elevation Profile on Dashboard
            ax2 = fig.add_subplot(gs[0, 2:])
            for _, row in paths_data.head(3).iterrows():
                path_id = row['PathID']
                path = gdf[gdf['PathID'] == path_id].sort_values('Time')
                distances = [0]
                for i in range(1, len(path)):
                    distances.append(distances[-1] + path.iloc[i-1].geometry.distance(path.iloc[i].geometry))
                ax2.plot(distances, path['Z'].values, label=f'Path {path_id}')
            ax2.set_title('Elevation Profiles (Top 3 Paths)')
            ax2.set_xlabel('Distance (m)')
            ax2.set_ylabel('Elevation (m)')
            ax2.legend()

            # Add more dashboard components as needed

            plt.tight_layout()
            save_plot(fig, 'dashboard')
        except Exception as e:
            logging.error(f"Error creating dashboard: {e}")
            print(f"Error creating dashboard: {e}")

        logging.info(f"Visualizations saved in {plots_folder}")
        print(Fore.GREEN + f"Visualizations saved in {plots_folder}" + Style.RESET_ALL)

def generate_final_report(
    input_shapefile: str,
    output_shp: str,
    output_kml: str,
    metrics: Dict[str, Dict[str, pd.DataFrame]],
    output_folder: str
) -> None:
    """
    Generates a comprehensive final report with all analysis results.

    Args:
        input_shapefile (str): Path to the input shapefile (if applicable).
        output_shp (str): Path to the output shapefile.
        output_kml (str): Path to the output KML file.
        metrics (Dict[str, Dict[str, pd.DataFrame]]): Dictionary containing all calculated metrics.
        output_folder (str): Base output folder path.
    """
    try:
        report_lines = []
        
        def add_section(title: str, content: List[str]):
            report_lines.extend([
                f"\n{'='*20} {title} {'='*20}\n",
                *content,
                "\n"
            ])

        # Basic Information
        basic_info = [
            f"Analysis Date and Time: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}",
            f"Input File: {input_shapefile if input_shapefile else 'Multiple CSV/TXT files'}",
            f"Output Shapefile: {output_shp}",
            f"Output KML: {output_kml}",
        ]
        add_section("Basic Information", basic_info)

        # Path Statistics
        paths_data = metrics.get('paths', {}).get('paths', None)
        if paths_data is not None:
            path_stats = [
                f"Total number of paths: {len(paths_data)}",
                f"Total length: {paths_data['length_km'].sum():.3f} km",
                f"Average length: {paths_data['length_km'].mean():.3f} km",
                f"Minimum length: {paths_data['length_km'].min():.3f} km",
                f"Maximum length: {paths_data['length_km'].max():.3f} km",
                f"Average points per km: {paths_data['pts_per_km'].mean():.1f} points/km"
            ]
            add_section("Path Statistics", path_stats)

        # Temporal Analysis
        temporal_data = metrics.get('temporal', {}).get('coverage', pd.DataFrame())
        if not temporal_data.empty:
            temporal_stats = [
                f"Total acquisition duration: {temporal_data['total_duration_hours'].iloc[0]:.2f} hours",
                f"Active time: {temporal_data['active_time_hours'].iloc[0]:.2f} hours",
                f"Acquisition start time: {temporal_data['start_time'].iloc[0].strftime('%Y-%m-%d %H:%M:%S')}",
                f"Acquisition end time: {temporal_data['end_time'].iloc[0].strftime('%Y-%m-%d %H:%M:%S')}"
            ]
            add_section("Temporal Analysis", temporal_stats)

        # Gap Analysis
        gap_data = metrics.get('temporal', {}).get('gap_stats', pd.DataFrame())
        if not gap_data.empty:
            gap_stats = [
                f"Minimum gap: {gap_data['min_gap'].iloc[0]:.3f} seconds",
                f"Maximum gap: {gap_data['max_gap'].iloc[0]:.3f} seconds",
                f"Mean gap: {gap_data['mean_gap'].iloc[0]:.3f} seconds",
                f"Median gap: {gap_data['median_gap'].iloc[0]:.3f} seconds",
                f"Standard Deviation gap: {gap_data['std_gap'].iloc[0]:.3f} seconds"
            ]
            add_section("Gap Analysis", gap_stats)

        # Elevation Analysis
        elevation_data = metrics.get('elevation', {}).get('summary', pd.DataFrame())
        if not elevation_data.empty:
            elevation_stats = [
                f"Minimum elevation: {elevation_data['min_elevation'].iloc[0]:.2f} m",
                f"Maximum elevation: {elevation_data['max_elevation'].iloc[0]:.2f} m",
                f"Elevation range: {elevation_data['elevation_range'].iloc[0]:.2f} m",
                f"Average elevation: {elevation_data['mean_elevation'].iloc[0]:.2f} m",
                f"Total elevation gain: {elevation_data['total_elevation_gain_m'].iloc[0]:.2f} m",
                f"Total elevation loss: {elevation_data['total_elevation_loss_m'].iloc[0]:.2f} m",
                f"Mean slope: {elevation_data['mean_slope_deg'].iloc[0]:.2f}°",
                f"Maximum slope: {elevation_data['max_slope_deg'].iloc[0]:.2f}°"
            ]
            add_section("Elevation Analysis", elevation_stats)

        # Photo Coverage
        photo_data = metrics.get('photos', {}).get('summary', pd.DataFrame())
        if not photo_data.empty:
            photo_stats = [
                f"Total number of photos: {photo_data['total_photos'].iloc[0]}",
                f"Number of folders: {photo_data['total_folders'].iloc[0]}",
                f"Average photos per path: {photo_data['mean_photos_per_path'].iloc[0]:.1f}",
                f"Average photo spacing: {photo_data['mean_photo_spacing_m'].iloc[0]:.2f} m",
                f"Overall photo density: {photo_data['overall_photo_density_photos_per_m2'].iloc[0]:.6f} photos/m²"
            ]
            add_section("Photo Coverage", photo_stats)

        # Camera Orientation
        camera_data = metrics.get('camera', {}).get('summary', pd.DataFrame())
        if not camera_data.empty:
            camera_stats = [
                f"Overall Roll STD: {camera_data['overall_roll_std'].iloc[0]:.3f}°",
                f"Overall Pitch STD: {camera_data['overall_pitch_std'].iloc[0]:.3f}°",
                f"Total sharp turns: {camera_data['total_sharp_turns'].iloc[0]}",
                f"Average stability index: {camera_data['mean_stability_index'].iloc[0]:.3f}",
                f"Worst stability index: {camera_data['worst_stability_index'].iloc[0]:.3f}"
            ]
            add_section("Camera Orientation", camera_stats)

        # Path Complexity
        complexity_data = metrics.get('complexity', {}).get('summary', pd.DataFrame())
        if not complexity_data.empty:
            complexity_stats = [
                f"Total sharp turns: {complexity_data['total_sharp_turns'].iloc[0]}",
                f"Mean sinuosity: {complexity_data['mean_sinuosity'].iloc[0]:.3f}",
                f"Maximum sinuosity: {complexity_data['max_sinuosity'].iloc[0]:.3f}",
                f"Mean curvature index: {complexity_data['mean_curvature_index'].iloc[0]:.3f}"
            ]
            add_section("Path Complexity", complexity_stats)

        # Overlap Analysis
        overlaps_present = 'overlaps' in metrics and 'summary' in metrics['overlaps']
        if overlaps_present:
            overlap_data = metrics['overlaps']['summary']
            if not overlap_data.empty:
                overlap_stats = [
                    f"Number of overlapping paths: {overlap_data['total_overlapping_paths'].iloc[0]}",
                    f"Total overlap length: {overlap_data['total_overlap_length_m'].iloc[0]:.2f} m",
                    f"Average overlap percentage: {overlap_data['mean_overlap_percentage'].iloc[0]:.2f}%",
                    f"Maximum overlap percentage: {overlap_data['max_overlap_percentage'].iloc[0]:.2f}%"
                ]
                add_section("Overlap Analysis", overlap_stats)

        # Detailed Path Statistics
        if paths_data is not None:
            detailed_stats = []
            for _, path in paths_data.iterrows():
                path_details = [
                    f"\nPath {int(path['PathID'])}:",
                    f"  Length: {path['length_km']:.3f} km",
                    f"  Duration: {path['dur_s']:.2f} s",
                    f"  Average Speed: {path['avg_spd']:.2f} km/h",
                    f"  Points per km: {path['pts_per_km']:.1f}",
                    f"  Photos per km: {path['phs_per_km']:.1f}",
                    f"  Camera Stability Index: {path['cam_sta']:.3f}"
                ]
                detailed_stats.extend(path_details)
            add_section("Detailed Path Statistics", detailed_stats)

        # Save report
        report_file = os.path.join(output_folder, 'reports', 'final_report.txt')
        os.makedirs(os.path.dirname(report_file), exist_ok=True)
        
        with open(report_file, 'w', encoding='utf-8') as f:
            f.write('\n'.join(report_lines))

        logging.info(f"Final report saved at {report_file}")
        print(Fore.GREEN + f"Final report saved at {report_file}" + Style.RESET_ALL)

    except Exception as e:
        logging.error(f"Error generating final report: {e}")
        print(Fore.RED + f"Error generating final report: {e}" + Style.RESET_ALL)

def create_crs_prompt(crs_prompt: str = "Coordinate Reference System (e.g., 32632)", default_crs: str = "3857") -> str:
    """
    Creates a simplified CRS prompt allowing users to input just the EPSG code without 'EPSG:'.

    Args:
        crs_prompt (str, optional): The input prompt to display. Defaults to "Coordinate Reference System (e.g., 32632)".
        default_crs (str, optional): The default EPSG code. Defaults to "3857".

    Returns:
        str: The formatted CRS string.
    """
    user_input = get_input(crs_prompt, default_crs, cast_type=str).strip()
    if not user_input.upper().startswith('EPSG:'):
        user_input = f"EPSG:{user_input}"
    return user_input.upper()

def main():
    """
    Main function executing the complete workflow.
    """
    print(Fore.BLUE + "\n=== Mobile Mapping System Data Analyzer ===\n" + Style.RESET_ALL)

    try:
        # Setup output folders
        output_base = os.path.join(os.getcwd(), OUTPUT_FOLDER)
        folders = setup_output_folders(output_base)

        # Get processing mode
        create_shapefile = get_input(
            "Do you want to create a shapefile from CSV/TXT files? (y/n)",
            default='y',
            cast_type=str
        ).lower()

        if create_shapefile == 'y':
            # Get input folder
            input_folder = get_input(
                "Input files folder",
                default=None,
                cast_type=str
            ).strip('"').strip("'")

            if not os.path.isdir(input_folder):
                raise ValueError("Input folder does not exist.")

            # Get file extension
            file_extension = get_input(
                "File extension (csv/txt)",
                default='txt',
                cast_type=str
            ).lower()

            if file_extension not in ['csv', 'txt']:
                raise ValueError("Invalid file extension. Use 'csv' or 'txt'.")

            # Get CRS
            crs = create_crs_prompt()

            # Process files
            merged_df, file_count = merge_files(input_folder, file_extension)
            gdf = dataframe_to_geodataframe(merged_df, 'X', 'Y', crs)

            print("Available columns:", gdf.columns.tolist())

            # Get gap threshold and assign paths
            gap_threshold = get_input(
                "Gap threshold (seconds)",
                default=3.0,
                cast_type=float
            )

            logging.info(f"Sorting GeoDataFrame by Time.")
            gdf_sorted = gdf.sort_values('Time')
            gdf_sorted, gaps_df = assign_path_ids(gdf_sorted, gap_threshold)

            print("\nAfter assigning PathIDs:")
            print("Unique PathIDs:", sorted(gdf_sorted['PathID'].unique()))
            print("Number of paths:", gdf_sorted['PathID'].nunique())

            # Set output paths
            folder_name = os.path.basename(os.path.normpath(input_folder))
            output_base_name = get_input(
                "Output file name (without extension)",
                default=f"{folder_name}_processed",
                cast_type=str
            )
            output_shp = os.path.join(folders['stats'], f"{output_base_name}.shp")
            output_kml = os.path.join(folders['stats'], f"{output_base_name}.kml")

            try:
                # Calculate all metrics
                print("\nCalculating metrics...")
                metrics = {}

                # Path Metrics
                path_metrics = calculate_path_metrics(gdf_sorted)
                metrics.update(path_metrics)

                if 'paths' in path_metrics:
                    # Save Shapefile and KML
                    print("\nSaving outputs...")
                    # Handle field name truncation by renaming columns in DataFrame
                    paths_gdf = path_metrics['paths']
                    paths_gdf.to_file(output_shp)
                    paths_gdf.to_crs(epsg=4326).to_file(output_kml, driver='KML')
                    logging.info(f"Shapefile saved at {output_shp}")
                    logging.info(f"KML file saved at {output_kml}")
                    print(Fore.GREEN + f"Files saved: {output_shp} and {output_kml}" + Style.RESET_ALL)

                # Additional Metrics
                print("\nCalculating additional metrics...")
                metrics['temporal'] = analyze_temporal_patterns(gdf_sorted)
                metrics['spatial'] = analyze_spatial_patterns(gdf_sorted)
                metrics['elevation'] = analyze_elevation_profile(gdf_sorted)
                metrics['photos'] = analyze_photo_coverage(gdf_sorted)
                metrics['camera'] = analyze_camera_orientation(gdf_sorted)
                metrics['complexity'] = calculate_path_complexity(gdf_sorted)
                metrics['overlaps'] = calculate_overlap_analysis(gdf_sorted)

                # Enrich GeoDataFrame with gaps data if needed
                metrics['gaps'] = gaps_df

                # Create visualizations
                print("\nCreating visualizations...")
                create_visualizations(gdf_sorted, metrics, output_base_name)

                # Generate final report
                print("\nGenerating final report...")
                generate_final_report(None, output_shp, output_kml, metrics, output_base_name)

                print(Fore.GREEN + "\nAnalysis completed successfully!" + Style.RESET_ALL)

            except Exception as process_error:
                logging.error(f"Error during processing: {process_error}")
                print(Fore.RED + f"Error during processing: {process_error}" + Style.RESET_ALL)
                raise

    except Exception as e:
        logging.error(f"Execution error: {e}")
        print(Fore.RED + f"Error: {e}" + Style.RESET_ALL)

if __name__ == "__main__":
    main()