Add coordinate system detection for Argus GeoTIFF imagery

[Copilot]

Argus GeoTIFF imagery uses NC State Plane (NAD83) coordinates, not lon/lat. The existing State Plane → Lon/Lat conversion is correct, but there was no automated detection or validation of coordinate systems, leading to confusion.

Changes

• Added detect_geotiff_crs() function - Analyzes coordinate ranges to return 'state_plane', 'lonlat', or 'unknown'

• Enhanced get_geotiff_extent() documentation - Clarified that it returns coordinates in the GeoTIFF's native CRS with explicit warnings and conversion examples

• Updated example scripts - Added CRS detection, conditional conversion, and extent validation with warnings for suspicious ranges

• Added unit tests - 4 tests covering State Plane, lon/lat, and unknown CRS detection scenarios

Usage

from murgtools.getdata import get_geotiff_extent, detect_geotiff_crs
from murgtools.utils import geoprocess as gp

# Detect and convert appropriately
crs_type = detect_geotiff_crs('argus_image.tif')

if crs_type == 'state_plane':
    extent = get_geotiff_extent('argus_image.tif')
    ll = gp.FRFcoord(extent[0], extent[2], coordType='ncsp')
    ur = gp.FRFcoord(extent[1], extent[3], coordType='ncsp')
    lonlat_extent = [ll['Lon'], ur['Lon'], ll['Lat'], ur['Lat']]
elif crs_type == 'lonlat':
    lonlat_extent = get_geotiff_extent('argus_image.tif')

The State Plane → Lon/Lat conversion logic is preserved unchanged - this adds detection and validation layers only.

Original prompt

Fix Argus Imagery Coordinate System Detection and Proper Conversion

Problem Analysis

Based on diagnostic results from examples/diagnose_coordinates.py, we've confirmed:

1. ✅ Argus GeoTIFF imagery uses NC State Plane (NAD83) coordinates
2. ✅ The State Plane → Lon/Lat conversion is working correctly
3. ✅ Extents do overlap after conversion

However, PR #48 was created reporting alignment issues. The diagnostic shows the coordinate conversion itself is correct, so the issue must be in how the imagery is being displayed or used in examples/test_wave_and_imagery.py.

Root Cause

Looking at examples/test_wave_and_imagery.py lines 130-144, the code:

1. Extracts State Plane extent from Argus GeoTIFF
2. Converts corners to lon/lat
3. Builds extent for matplotlib

The problem is likely in the extent format or plotting logic, not the coordinate conversion itself.

Required Changes

1. Keep the State Plane Conversion (It's Correct!)

The conversion from State Plane to lon/lat in test_wave_and_imagery.py is correct and must be kept. Do NOT remove it as PR #48 suggests.

2. Fix the Real Issue

The actual problem appears to be in how the extent is being used or displayed. Review test_wave_and_imagery.py lines 145-175 to ensure:

• Matplotlib extent parameter format is correct: [left, right, bottom, top] in data coordinates
• Image aspect ratio is being handled properly
• Coordinate reference systems are consistent across all plotting elements

3. Add Coordinate System Detection

Add a helper function to detect and validate coordinate systems:

def detect_geotiff_crs(filepath):
    """Detect coordinate reference system of a GeoTIFF.
    
    Returns:
        str: 'state_plane', 'lonlat', or 'unknown'
    """
    import tifffile
    
    extent = get_geotiff_extent(filepath)
    left, right, bottom, top = extent
    
    # NC State Plane NAD83 typical ranges at FRF
    if (800000 < left < 1000000 and 200000 < bottom < 300000):
        return 'state_plane'
    # Geographic coordinates typical ranges
    elif (-180 <= left <= 180 and -90 <= bottom <= 90):
        return 'lonlat'
    else:
        return 'unknown'

4. Update test_wave_and_imagery.py

In examples/test_wave_and_imagery.py:

1. Keep the State Plane → Lon/Lat conversion (lines 130-144)
2. Add CRS detection to confirm assumptions
3. Add validation that extent is reasonable before plotting
4. Add better error messages when coordinates don't make sense

Example fix around line 130:

# Extract extent from the GeoTIFF
sp_extent = get_geotiff_extent(tmp_path)
print(f"  Raw GeoTIFF extent: {sp_extent}")

# Detect coordinate system
crs_type = detect_geotiff_crs(tmp_path)
print(f"  Detected CRS: {crs_type}")

if crs_type == 'state_plane':
    # Convert State Plane corners to lat/lon (THIS IS CORRECT!)
    sp_left, sp_right, sp_bottom, sp_top = sp_extent
    
    # Convert corners from State Plane to lat/lon
    ll_corner = gp.FRFcoord(sp_left, sp_bottom, coordType='ncsp')  # SW corner
    ur_corner = gp.FRFcoord(sp_right, sp_top, coordType='ncsp')    # NE corner
    
    # Build lat/lon extent [left, right, bottom, top] as [lon_min, lon_max, lat_min, lat_max]
    argus_extent = [ll_corner['Lon'], ur_corner['Lon'], ll_corner['Lat'], ur_corner['Lat']]
    print(f"  Converted to lon/lat: {argus_extent}")
    
elif crs_type == 'lonlat':
    # Already in lon/lat, use directly
    argus_extent = sp_extent
    print(f"  Using lon/lat extent directly: {argus_extent}")
    
else:
    print(f"  ✗ Unknown CRS type for extent: {sp_extent}")
    argus_extent = None

# Validate extent is reasonable
if argus_extent:
    lon_range = argus_extent[1] - argus_extent[0]
    lat_range = argus_extent[3] - argus_extent[2]
    
    if lon_range > 1 or lat_range > 1:  # More than 1 degree is suspicious
        print(f"  ⚠ Warning: Extent spans {lon_range:.3f}° lon, {lat_range:.3f}° lat")
        print(f"  ⚠ This seems too large, coordinates may be incorrect")

5. Add Documentation

Add a docstring to get_geotiff_extent() explaining coordinate systems:

def get_geotiff_extent(filepath):
    """Extract matplotlib extent from GeoTIFF using tifffile.

    Parses GeoTIFF tags (ModelTiepointTag and ModelPixelScaleTag) to compute
    geographic bounds suitable for matplotlib imshow extent parameter.
    
    ⚠️ IMPORTANT: This function returns coordinates in whatever CRS the GeoTIFF
    uses. For Argus imagery at FRF, this is typically NC State Plane (NAD83).
    You MUST convert to lon/lat before overlaying with satellite imagery!

    Args:
        filepath (str): Path to GeoTIFF file.

    Returns:
        list: [left, right, bottom, top] extent in the GeoTIFF's native CRS.
            For Argus imagery: NC State Plane meters (easting, northing)
            For satellite imagery: Geographic degrees (lon, lat)

    Example:
        >>> extent = get_geotiff_extent('argus_image.tif')
    ...

</details>



<!-- START COPILOT CODING AGENT SUFFIX -->

*This pull request was created from Copilot chat.*
>

<!-- START COPILOT CODING AGENT TIPS -->
---

✨ Let Copilot coding agent [set things up for you](https://github.com/SBFRF/murgtools/issues/new?title=✨+Set+up+Copilot+instructions&body=Configure%20instructions%20for%20this%20repository%20as%20documented%20in%20%5BBest%20practices%20for%20Copilot%20coding%20agent%20in%20your%20repository%5D%28https://gh.io/copilot-coding-agent-tips%29%2E%0A%0A%3COnboard%20this%20repo%3E&assignees=copilot) — coding agent works faster and does higher quality work when set up for your repo.