diff --git a/examples/diagnose_coordinates.py b/examples/diagnose_coordinates.py index 1433e65..904f81b 100644 --- a/examples/diagnose_coordinates.py +++ b/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)") diff --git a/examples/test_wave_and_imagery.py b/examples/test_wave_and_imagery.py index 3ab0b86..504ae37 100644 --- a/examples/test_wave_and_imagery.py +++ b/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}") diff --git a/murgtools/getdata/__init__.py b/murgtools/getdata/__init__.py index 1bd1428..3e3ca7b 100644 --- a/murgtools/getdata/__init__.py +++ b/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", diff --git a/murgtools/getdata/getDataFRF.py b/murgtools/getdata/getDataFRF.py index 302616a..8e5e011 100755 --- a/murgtools/getdata/getDataFRF.py +++ b/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. diff --git a/tests/test_getDataFRF.py b/tests/test_getDataFRF.py index fb85eae..7a41638 100644 --- a/tests/test_getDataFRF.py +++ b/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)