BE-LightningReport/utm_ed50_to_wgs84_converter.py

#!/usr/bin/env python3
"""
UTM ED50 to WGS84 Coordinate Converter

This script converts UTM coordinates in 6-degree zones from ED50 (European Datum 1950)
reference system to WGS84 format.

Requirements:
    pip install pyproj pandas

Usage:
    python utm_ed50_to_wgs84_converter.py input_file.csv output_file.csv
    python utm_ed50_to_wgs84_converter.py --interactive
"""

import argparse
import pandas as pd
import pyproj
from typing import Tuple, List, Optional
import sys


class UTMED50ToWGS84Converter:
    """Converter for UTM ED50 coordinates to WGS84."""
    
    def __init__(self):
        """Initialize the converter with ED50 and WGS84 projections."""
        self.ed50_utm_projections = {}
        self.wgs84_utm_projections = {}
        
    def get_ed50_utm_projection(self, zone: int, northern: bool = True) -> pyproj.Proj:
        """Get ED50 UTM projection for a specific zone."""
        key = (zone, northern)
        if key not in self.ed50_utm_projections:
            hemisphere = 'N' if northern else 'S'
            proj_string = f"+proj=utm +zone={zone} +ellps=intl +towgs84=-87,-98,-121,0,0,0,0 +units=m +no_defs"
            self.ed50_utm_projections[key] = pyproj.Proj(proj_string)
        return self.ed50_utm_projections[key]
    
    def get_wgs84_utm_projection(self, zone: int, northern: bool = True) -> pyproj.Proj:
        """Get WGS84 UTM projection for a specific zone."""
        key = (zone, northern)
        if key not in self.wgs84_utm_projections:
            hemisphere = 'N' if northern else 'S'
            proj_string = f"+proj=utm +zone={zone} +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
            self.wgs84_utm_projections[key] = pyproj.Proj(proj_string)
        return self.wgs84_utm_projections[key]
    
    def convert_single_point(self, easting: float, northing: float, zone: int, 
                           northern: bool = True) -> Tuple[float, float]:
        """
        Convert a single UTM ED50 point to WGS84 lat/lon.
        
        Args:
            easting: UTM easting coordinate in meters
            northing: UTM northing coordinate in meters
            zone: UTM zone number (1-60)
            northern: True if in northern hemisphere, False if southern
            
        Returns:
            Tuple of (latitude, longitude) in WGS84 decimal degrees
        """
        if not (1 <= zone <= 60):
            raise ValueError(f"Invalid UTM zone: {zone}. Must be between 1 and 60.")
        
        # Create ED50 UTM projection
        ed50_proj = self.get_ed50_utm_projection(zone, northern)
        
        # Create WGS84 lat/lon projection
        wgs84_latlon = pyproj.Proj('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
        
        # Create transformer from ED50 UTM to WGS84 lat/lon
        transformer = pyproj.Transformer.from_proj(ed50_proj, wgs84_latlon)
        
        # Transform coordinates directly to WGS84 lat/lon
        wgs84_lon, wgs84_lat = transformer.transform(easting, northing)
        
        return wgs84_lat, wgs84_lon
    
    def convert_dataframe(self, df: pd.DataFrame, easting_col: str, northing_col: str, 
                         zone_col: str, northern_col: Optional[str] = None,
                         northern_default: bool = True) -> pd.DataFrame:
        """
        Convert UTM ED50 coordinates in a DataFrame to WGS84.
        
        Args:
            df: Input DataFrame with UTM coordinates
            easting_col: Column name for easting coordinates
            northing_col: Column name for northing coordinates
            zone_col: Column name for UTM zone
            northern_col: Column name for northern hemisphere flag (optional)
            northern_default: Default value for northern hemisphere if column not provided
            
        Returns:
            DataFrame with additional WGS84 lat/lon columns
        """
        result_df = df.copy()
        result_df['wgs84_lat'] = None
        result_df['wgs84_lon'] = None
        
        for idx, row in df.iterrows():
            try:
                easting = float(row[easting_col])
                northing = float(row[northing_col])
                zone = int(row[zone_col])
                
                if northern_col and northern_col in row:
                    northern = bool(row[northern_col])
                else:
                    northern = northern_default
                
                lat, lon = self.convert_single_point(easting, northing, zone, northern)
                result_df.at[idx, 'wgs84_lat'] = lat
                result_df.at[idx, 'wgs84_lon'] = lon
                
            except (ValueError, KeyError) as e:
                print(f"Warning: Could not convert row {idx}: {e}")
                result_df.at[idx, 'wgs84_lat'] = None
                result_df.at[idx, 'wgs84_lon'] = None
        
        return result_df


def interactive_mode():
    """Run the converter in interactive mode."""
    converter = UTMED50ToWGS84Converter()
    
    print("UTM ED50 to WGS84 Coordinate Converter")
    print("=" * 40)
    print("Enter coordinates (type 'quit' to exit)")
    print()
    
    while True:
        try:
            # Get input
            easting_input = input("Enter easting (meters): ").strip()
            if easting_input.lower() == 'quit':
                break
            
            northing_input = input("Enter northing (meters): ").strip()
            if northing_input.lower() == 'quit':
                break
            
            zone_input = input("Enter UTM zone (1-60): ").strip()
            if zone_input.lower() == 'quit':
                break
            
            hemisphere_input = input("Enter hemisphere (N/S) [default: N]: ").strip().upper()
            if hemisphere_input == 'QUIT':
                break
            
            # Parse inputs
            easting = float(easting_input)
            northing = float(northing_input)
            zone = int(zone_input)
            northern = hemisphere_input != 'S' if hemisphere_input else True
            
            # Convert
            lat, lon = converter.convert_single_point(easting, northing, zone, northern)
            
            print(f"\nWGS84 Coordinates:")
            print(f"Latitude:  {lat:.8f}°")
            print(f"Longitude: {lon:.8f}°")
            print("-" * 40)
            
        except ValueError as e:
            print(f"Error: {e}")
            print("Please enter valid numeric values.")
        except KeyboardInterrupt:
            print("\nExiting...")
            break


def main():
    """Main function to handle command line arguments and file processing."""
    parser = argparse.ArgumentParser(
        description="Convert UTM ED50 coordinates to WGS84 format",
        formatter_class=argparse.RawDescriptionHelpFormatter,
        epilog="""
Examples:
  # Interactive mode
  python utm_ed50_to_wgs84_converter.py --interactive
  
  # Convert CSV file
  python utm_ed50_to_wgs84_converter.py input.csv output.csv
  
  # Convert with custom column names
  python utm_ed50_to_wgs84_converter.py input.csv output.csv --easting-col X --northing-col Y --zone-col ZONE
        """
    )
    
    parser.add_argument('input_file', nargs='?', help='Input CSV file with UTM coordinates')
    parser.add_argument('output_file', nargs='?', help='Output CSV file for WGS84 coordinates')
    parser.add_argument('--interactive', '-i', action='store_true', 
                       help='Run in interactive mode')
    parser.add_argument('--easting-col', default='easting', 
                       help='Column name for easting coordinates (default: easting)')
    parser.add_argument('--northing-col', default='northing', 
                       help='Column name for northing coordinates (default: northing)')
    parser.add_argument('--zone-col', default='zone', 
                       help='Column name for UTM zone (default: zone)')
    parser.add_argument('--northern-col', 
                       help='Column name for northern hemisphere flag (optional)')
    parser.add_argument('--northern-default', action='store_true', default=True,
                       help='Default value for northern hemisphere (default: True)')
    
    args = parser.parse_args()
    
    if args.interactive:
        interactive_mode()
        return
    
    if not args.input_file or not args.output_file:
        parser.error("Both input_file and output_file are required when not in interactive mode")
    
    try:
        # Read input file
        print(f"Reading input file: {args.input_file}")
        df = pd.read_csv(args.input_file)
        
        # Check required columns
        required_cols = [args.easting_col, args.northing_col, args.zone_col]
        missing_cols = [col for col in required_cols if col not in df.columns]
        if missing_cols:
            print(f"Error: Missing required columns: {missing_cols}")
            print(f"Available columns: {list(df.columns)}")
            sys.exit(1)
        
        # Convert coordinates
        print("Converting coordinates...")
        converter = UTMED50ToWGS84Converter()
        result_df = converter.convert_dataframe(
            df, 
            args.easting_col, 
            args.northing_col, 
            args.zone_col,
            args.northern_col,
            args.northern_default
        )
        
        # Save output
        print(f"Saving output file: {args.output_file}")
        result_df.to_csv(args.output_file, index=False)
        
        # Print summary
        total_rows = len(result_df)
        successful_conversions = result_df['wgs84_lat'].notna().sum()
        print(f"\nConversion complete!")
        print(f"Total rows: {total_rows}")
        print(f"Successful conversions: {successful_conversions}")
        print(f"Failed conversions: {total_rows - successful_conversions}")
        
    except FileNotFoundError:
        print(f"Error: Input file '{args.input_file}' not found")
        sys.exit(1)
    except Exception as e:
        print(f"Error: {e}")
        sys.exit(1)


if __name__ == "__main__":
    main()