Add coordinate system detection for Argus GeoTIFF imagery

examples/diagnose_coordinates.py

@@ -11,7 +11,7 @@
 
 import datetime as DT
 import numpy as np
-from murgtools.getdata import getSatelliteImagery, findArgusImagery, get_geotiff_extent
+from murgtools.getdata import getSatelliteImagery, findArgusImagery, get_geotiff_extent, detect_geotiff_crs
 from murgtools.utils import geoprocess as gp
 import tempfile
 
@@ -98,12 +98,13 @@ def diagnose_coordinates():
             print(f"  Bottom: {raw_extent[2]:.6f}")
             print(f"  Top:    {raw_extent[3]:.6f}")
 
-            # Determine coordinate system
+            # Determine coordinate system using the new function
             print(f"\nCoordinate System Detection:")
             left, right, bottom, top = raw_extent
 
-            # Check if values look like State Plane
-            is_state_plane = (left > 800000 and bottom > 200000)
+            # Use the detect_geotiff_crs function
+            crs_type = detect_geotiff_crs(tmp_path)
+            is_state_plane = (crs_type == 'state_plane')
 
             if is_state_plane:
                 print(f"  ✓ DETECTED: North Carolina State Plane (NAD83)")

examples/test_wave_and_imagery.py

@@ -11,9 +11,13 @@
 import matplotlib.dates as mdates
 import numpy as np
 
-from murgtools.getdata import getObs, getSatelliteImagery, findArgusImagery, get_geotiff_extent
+from murgtools.getdata import getObs, getSatelliteImagery, findArgusImagery, get_geotiff_extent, detect_geotiff_crs
 from murgtools.utils import geoprocess as gp
 
+# Threshold for suspicious extent size (degrees)
+# More than 1 degree of lon/lat span is suspicious for FRF imagery
+SUSPICIOUS_EXTENT_THRESHOLD_DEGREES = 1.0
+
 
 def main():
     # Define time range: last 30 days with clean datetime boundaries
@@ -127,21 +131,43 @@ def main():
             offset_mins = argus_result.get('time_offset_minutes', 0)
             if offset_mins != 0:
                 print(f"  Time offset from requested: {offset_mins} minutes ({offset_mins/60:.1f} hours)")
-            # Extract extent from the GeoTIFF (in State Plane coordinates)
+            # Extract extent from the GeoTIFF
             sp_extent = get_geotiff_extent(tmp_path)
-            print(f"  State Plane extent: {sp_extent}")
+            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
 
-            # Convert State Plane corners to lat/lon
-            # sp_extent is [left, right, bottom, top] in State Plane Easting/Northing
-            sp_left, sp_right, sp_bottom, sp_top = sp_extent
+                # 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
 
-            # 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
+            # Validate extent is reasonable
+            if argus_extent:
+                lon_range = argus_extent[1] - argus_extent[0]
+                lat_range = argus_extent[3] - argus_extent[2]
 
-            # 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"  Lat/Lon extent: {argus_extent}")
+                if lon_range > SUSPICIOUS_EXTENT_THRESHOLD_DEGREES or lat_range > SUSPICIOUS_EXTENT_THRESHOLD_DEGREES:
+                    print(f"  ⚠ Warning: Extent spans {lon_range:.3f}° lon, {lat_range:.3f}° lat")
+                    print(f"  ⚠ This seems too large, coordinates may be incorrect")
     except Exception as e:
         print(f"  Error fetching Argus imagery: {e}")
 

murgtools/getdata/__init__.py

@@ -6,7 +6,7 @@
 """
 
 from .getDataFRF import (getObs, getDataTestBed, gettime, getnc, removeDuplicatesFromDictionary,
-                         get_geotiff_extent, getArgusImagery, threadGetArgusImagery, findArgusImagery,
+                         get_geotiff_extent, detect_geotiff_crs, getArgusImagery, threadGetArgusImagery, findArgusImagery,
                          getArgusPixelIntensity)
 from .getOutsideData import forecastData, getSatelliteImagery
 from .getPlotData import alt_PlotData
@@ -24,6 +24,7 @@
     "forecastData",
     "getSatelliteImagery",
     "get_geotiff_extent",
+    "detect_geotiff_crs",
     "getArgusImagery",
     "threadGetArgusImagery",
     "findArgusImagery",

murgtools/getdata/getDataFRF.py

@@ -3348,20 +3348,28 @@ def get_geotiff_extent(filepath):
 
     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 geographic coordinates
-            (typically lon/lat or projected coordinates depending on the GeoTIFF).
+        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)
 
     Raises:
         KeyError: If required GeoTIFF tags are not present in the file.
 
     Example:
-        >>> extent = get_geotiff_extent('/path/to/image.tif')
-        >>> plt.imshow(image, extent=extent)
+        >>> extent = get_geotiff_extent('argus_image.tif')
+        >>> # extent is in State Plane! Convert to lon/lat:
+        >>> 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']]
 
     """
     import tifffile
@@ -3380,6 +3388,48 @@ def get_geotiff_extent(filepath):
         return [left, right, bottom, top]
 
 
+def detect_geotiff_crs(filepath):
+    """Detect coordinate reference system of a GeoTIFF.
+
+    Analyzes the coordinate ranges in a GeoTIFF to determine if it uses
+    NC State Plane (NAD83) coordinates or geographic (lon/lat) coordinates.
+    This is useful for validating coordinate system assumptions before plotting.
+
+    Args:
+        filepath (str): Path to GeoTIFF file.
+
+    Returns:
+        str: 'state_plane', 'lonlat', or 'unknown'
+            - 'state_plane': NC State Plane NAD83 (typical at FRF)
+            - 'lonlat': Geographic coordinates (longitude/latitude)
+            - 'unknown': Could not determine coordinate system
+
+    Example:
+        >>> crs = detect_geotiff_crs('argus_image.tif')
+        >>> if crs == 'state_plane':
+        >>>     # Need to convert to lon/lat
+        >>>     extent = get_geotiff_extent(filepath)
+        >>>     # ... perform conversion ...
+        >>> elif crs == 'lonlat':
+        >>>     # Already in lon/lat, use directly
+        >>>     extent = get_geotiff_extent(filepath)
+
+    """
+    extent = get_geotiff_extent(filepath)
+    left, right, bottom, top = extent
+
+    # NC State Plane NAD83 typical ranges at FRF
+    # Easting: ~900,000 meters, Northing: ~270,000 meters
+    if (800000 < left < 1000000 and 200000 < bottom < 300000):
+        return 'state_plane'
+    # Geographic coordinates typical ranges
+    # Longitude: -180 to 180, Latitude: -90 to 90
+    elif (-180 <= left <= 180 and -90 <= bottom <= 90):
+        return 'lonlat'
+    else:
+        return 'unknown'
+
+
 def getArgusImagery(dateOfInterest, filename=None, imageType="timex", verbose=True):
     """Retrieve Argus orthophoto imagery from the FRF coastal imaging server.
 

tests/test_getDataFRF.py

@@ -1250,3 +1250,152 @@ def test_location_as_dict(self):
 
             assert result is not None
             np.testing.assert_array_equal(result['intensity'][0], [255, 128, 64])
+
+
+class TestGeotiffCRS:
+    """Tests for GeoTIFF coordinate system detection and extent extraction."""
+
+    def test_detect_geotiff_crs_state_plane(self):
+        """Test detection of NC State Plane coordinates."""
+        from murgtools.getdata.getDataFRF import detect_geotiff_crs
+        import tempfile
+        import tifffile
+        import numpy as np
+
+        # Create a temporary GeoTIFF with State Plane coordinates
+        tmp = tempfile.NamedTemporaryFile(suffix='.tif', delete=False)
+        tmp_path = tmp.name
+        tmp.close()
+
+        try:
+            mock_image = np.zeros((100, 200, 3), dtype=np.uint8)
+
+            # NC State Plane coordinates (typical FRF values)
+            # Easting: ~901,000 meters, Northing: ~274,000 meters
+            tiepoint = [0.0, 0.0, 0.0, 901000.0, 274000.0, 0.0]
+            scale = [10.0, 10.0, 0.0]  # 10 meters per pixel
+
+            extratags = [
+                (33922, 'd', 6, tiepoint, True),
+                (33550, 'd', 3, scale, True),
+            ]
+
+            tifffile.imwrite(tmp_path, mock_image, extratags=extratags)
+
+            crs = detect_geotiff_crs(tmp_path)
+            assert crs == 'state_plane'
+        finally:
+            import os
+            if os.path.exists(tmp_path):
+                os.unlink(tmp_path)
+
+    def test_detect_geotiff_crs_lonlat(self):
+        """Test detection of lon/lat geographic coordinates."""
+        from murgtools.getdata.getDataFRF import detect_geotiff_crs
+        import tempfile
+        import tifffile
+        import numpy as np
+
+        # Create a temporary GeoTIFF with lon/lat coordinates
+        tmp = tempfile.NamedTemporaryFile(suffix='.tif', delete=False)
+        tmp_path = tmp.name
+        tmp.close()
+
+        try:
+            mock_image = np.zeros((100, 200, 3), dtype=np.uint8)
+
+            # Geographic coordinates (typical FRF values)
+            # Longitude: ~-75.75, Latitude: ~36.18
+            tiepoint = [0.0, 0.0, 0.0, -75.75, 36.18, 0.0]
+            scale = [0.0001, 0.0001, 0.0]  # ~10 meters per pixel at this latitude
+
+            extratags = [
+                (33922, 'd', 6, tiepoint, True),
+                (33550, 'd', 3, scale, True),
+            ]
+
+            tifffile.imwrite(tmp_path, mock_image, extratags=extratags)
+
+            crs = detect_geotiff_crs(tmp_path)
+            assert crs == 'lonlat'
+        finally:
+            import os
+            if os.path.exists(tmp_path):
+                os.unlink(tmp_path)
+
+    def test_detect_geotiff_crs_unknown(self):
+        """Test detection of unknown coordinate system."""
+        from murgtools.getdata.getDataFRF import detect_geotiff_crs
+        import tempfile
+        import tifffile
+        import numpy as np
+
+        # Create a temporary GeoTIFF with unrecognized coordinates
+        tmp = tempfile.NamedTemporaryFile(suffix='.tif', delete=False)
+        tmp_path = tmp.name
+        tmp.close()
+
+        try:
+            mock_image = np.zeros((100, 200, 3), dtype=np.uint8)
+
+            # Random coordinates that don't match known systems
+            tiepoint = [0.0, 0.0, 0.0, 5000.0, 3000.0, 0.0]
+            scale = [10.0, 10.0, 0.0]
+
+            extratags = [
+                (33922, 'd', 6, tiepoint, True),
+                (33550, 'd', 3, scale, True),
+            ]
+
+            tifffile.imwrite(tmp_path, mock_image, extratags=extratags)
+
+            crs = detect_geotiff_crs(tmp_path)
+            assert crs == 'unknown'
+        finally:
+            import os
+            if os.path.exists(tmp_path):
+                os.unlink(tmp_path)
+
+    def test_get_geotiff_extent_format(self):
+        """Test that get_geotiff_extent returns correct format."""
+        from murgtools.getdata.getDataFRF import get_geotiff_extent
+        import tempfile
+        import tifffile
+        import numpy as np
+
+        tmp = tempfile.NamedTemporaryFile(suffix='.tif', delete=False)
+        tmp_path = tmp.name
+        tmp.close()
+
+        try:
+            mock_image = np.zeros((100, 200, 3), dtype=np.uint8)
+
+            # State Plane coordinates
+            tiepoint = [0.0, 0.0, 0.0, 901000.0, 274000.0, 0.0]
+            scale = [10.0, 10.0, 0.0]
+
+            extratags = [
+                (33922, 'd', 6, tiepoint, True),
+                (33550, 'd', 3, scale, True),
+            ]
+
+            tifffile.imwrite(tmp_path, mock_image, extratags=extratags)
+
+            extent = get_geotiff_extent(tmp_path)
+
+            # Verify extent format: [left, right, bottom, top]
+            assert len(extent) == 4
+            assert extent[0] < extent[1]  # left < right
+            assert extent[2] < extent[3]  # bottom < top
+
+            # Verify calculated values
+            # left = 901000, right = 901000 + 200*10 = 903000
+            # top = 274000, bottom = 274000 - 100*10 = 273000
+            assert extent[0] == 901000.0
+            assert extent[1] == 903000.0
+            assert extent[2] == 273000.0
+            assert extent[3] == 274000.0
+        finally:
+            import os
+            if os.path.exists(tmp_path):
+                os.unlink(tmp_path)