From d3769f43cc03bff0fc2cb06e97d23d060a0d02c0 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Mon, 13 Jul 2026 13:23:26 +0000
Subject: [PATCH 01/13] Initial plan


From 57ca429036cf1a5ddeb2732a1ec2847d599a5d29 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Mon, 13 Jul 2026 13:27:57 +0000
Subject: [PATCH 02/13] Fix Argus overlay GeoTIFF extent conversion

---
 examples/test_wave_and_imagery.py | 17 ++-----
 murgtools/getdata/getDataFRF.py   | 45 ++++++++++++++----
 tests/test_getDataFRF.py          | 79 +++++++++++++++++++++++++++++++
 3 files changed, 120 insertions(+), 21 deletions(-)

diff --git a/examples/test_wave_and_imagery.py b/examples/test_wave_and_imagery.py
index 3ab0b86..151fbb1 100644
--- a/examples/test_wave_and_imagery.py
+++ b/examples/test_wave_and_imagery.py
@@ -127,20 +127,11 @@ 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)
-            sp_extent = get_geotiff_extent(tmp_path)
-            print(f"  State Plane extent: {sp_extent}")
+            native_extent = get_geotiff_extent(tmp_path)
+            print(f"  Native GeoTIFF extent: {native_extent}")
 
-            # 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
-
-            # 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']]
+            # Convert the Argus GeoTIFF bounds to lat/lon using the embedded CRS metadata.
+            argus_extent = get_geotiff_extent(tmp_path, to_latlon=True)
             print(f"  Lat/Lon extent: {argus_extent}")
     except Exception as e:
         print(f"  Error fetching Argus imagery: {e}")
diff --git a/murgtools/getdata/getDataFRF.py b/murgtools/getdata/getDataFRF.py
index 302616a..fa6f7b6 100755
--- a/murgtools/getdata/getDataFRF.py
+++ b/murgtools/getdata/getDataFRF.py
@@ -3343,7 +3343,7 @@ def getCSHOREOutput(self, prefix):
         return mod
 
 
-def get_geotiff_extent(filepath):
+def get_geotiff_extent(filepath, to_latlon=False):
     """Extract matplotlib extent from GeoTIFF using tifffile.
 
     Parses GeoTIFF tags (ModelTiepointTag and ModelPixelScaleTag) to compute
@@ -3351,33 +3351,62 @@ def get_geotiff_extent(filepath):
 
     Args:
         filepath (str): Path to GeoTIFF file.
+        to_latlon (bool, optional): When True, convert projected GeoTIFF bounds
+            to lon/lat using the GeoTIFF CRS metadata. Defaults to False.
 
     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 native
+            coordinates by default, or [lon_min, lon_max, lat_min, lat_max]
+            when ``to_latlon`` is True.
 
     Raises:
         KeyError: If required GeoTIFF tags are not present in the file.
+        ValueError: If ``to_latlon`` is True and the GeoTIFF CRS metadata
+            cannot be determined.
 
     Example:
         >>> extent = get_geotiff_extent('/path/to/image.tif')
         >>> plt.imshow(image, extent=extent)
+        >>> ll_extent = get_geotiff_extent('/path/to/image.tif', to_latlon=True)
 
     """
+    import pyproj
     import tifffile
 
     with tifffile.TiffFile(filepath) as tif:
-        tags = tif.pages[0].tags
+        page = tif.pages[0]
+        tags = page.tags
         # GeoTIFF tags: 33922=ModelTiepointTag, 33550=ModelPixelScaleTag
         tiepoint = tags[33922].value  # (i, j, k, x, y, z)
         scale = tags[33550].value     # (scaleX, scaleY, scaleZ)
-        height, width = tif.pages[0].shape[:2]
+        height, width = page.shape[:2]
         # Compute extent: [left, right, bottom, top]
         left = tiepoint[3]
-        top = tiepoint[4]
         right = left + width * scale[0]
-        bottom = top - height * scale[1]
-        return [left, right, bottom, top]
+        top = tiepoint[4]
+        edge_y = top - height * scale[1]
+        bottom, top = sorted((top, edge_y))
+
+        if not to_latlon:
+            return [min(left, right), max(left, right), bottom, top]
+
+        geotiff_tags = page.geotiff_tags or {}
+        epsg = geotiff_tags.get('ProjectedCSTypeGeoKey') or geotiff_tags.get('GeographicTypeGeoKey')
+        if hasattr(epsg, 'value'):
+            epsg = epsg.value
+        if epsg is None:
+            raise ValueError(f"GeoTIFF CRS metadata not found in {filepath}")
+
+        transformer = pyproj.Transformer.from_crs(f"EPSG:{epsg}", "EPSG:4326", always_xy=True)
+        corners = [
+            transformer.transform(left, bottom),
+            transformer.transform(left, top),
+            transformer.transform(right, bottom),
+            transformer.transform(right, top),
+        ]
+        lons = [corner[0] for corner in corners]
+        lats = [corner[1] for corner in corners]
+        return [min(lons), max(lons), min(lats), max(lats)]
 
 
 def getArgusImagery(dateOfInterest, filename=None, imageType="timex", verbose=True):
diff --git a/tests/test_getDataFRF.py b/tests/test_getDataFRF.py
index fb85eae..3bdbb5b 100644
--- a/tests/test_getDataFRF.py
+++ b/tests/test_getDataFRF.py
@@ -279,6 +279,85 @@ def test_returns_none_on_network_error(self):
             assert result is None
 
 
+class TestGetGeoTiffExtent:
+    """Tests for the get_geotiff_extent function."""
+
+    def test_returns_native_extent(self):
+        """Test native GeoTIFF extent extraction."""
+        import tempfile
+        import tifffile
+
+        from murgtools.getdata.getDataFRF import get_geotiff_extent
+
+        image = np.zeros((20, 10), dtype=np.uint8)
+        tiepoint = [0.0, 0.0, 0.0, 900500.0, 276000.0, 0.0]
+        scale = [10.0, 10.0, 0.0]
+
+        with tempfile.NamedTemporaryFile(suffix='.tif', delete=False) as tmp:
+            tmp_path = tmp.name
+
+        try:
+            tifffile.imwrite(
+                tmp_path,
+                image,
+                extratags=[
+                    (33922, 'd', 6, tiepoint, True),
+                    (33550, 'd', 3, scale, True),
+                ]
+            )
+
+            assert get_geotiff_extent(tmp_path) == [900500.0, 900600.0, 275800.0, 276000.0]
+        finally:
+            os.unlink(tmp_path)
+
+    def test_converts_projected_extent_to_latlon(self):
+        """Test projected GeoTIFF extent conversion to lon/lat."""
+        import tempfile
+        import tifffile
+        from pyproj import Transformer
+
+        from murgtools.getdata.getDataFRF import get_geotiff_extent
+
+        image = np.zeros((20, 10), dtype=np.uint8)
+        geokey_dir = (1, 1, 0, 3, 1024, 0, 1, 1, 1025, 0, 1, 1, 3072, 0, 1, 32618)
+        tiepoint = [0.0, 0.0, 0.0, 500000.0, 4000000.0, 0.0]
+        scale = [10.0, 10.0, 0.0]
+
+        with tempfile.NamedTemporaryFile(suffix='.tif', delete=False) as tmp:
+            tmp_path = tmp.name
+
+        try:
+            tifffile.imwrite(
+                tmp_path,
+                image,
+                extratags=[
+                    (33922, 'd', 6, tiepoint, True),
+                    (33550, 'd', 3, scale, True),
+                    (34735, 'H', len(geokey_dir), geokey_dir, True),
+                ]
+            )
+
+            extent = get_geotiff_extent(tmp_path, to_latlon=True)
+
+            transformer = Transformer.from_crs('EPSG:32618', 'EPSG:4326', always_xy=True)
+            corners = [
+                transformer.transform(500000.0, 3999800.0),
+                transformer.transform(500000.0, 4000000.0),
+                transformer.transform(500100.0, 3999800.0),
+                transformer.transform(500100.0, 4000000.0),
+            ]
+            expected = [
+                min(lon for lon, _ in corners),
+                max(lon for lon, _ in corners),
+                min(lat for _, lat in corners),
+                max(lat for _, lat in corners),
+            ]
+
+            assert extent == pytest.approx(expected)
+        finally:
+            os.unlink(tmp_path)
+
+
 class TestThreadGetArgusImagery:
     """Tests for the threadGetArgusImagery function."""
 

From 4c0eaace60de3096999cad443a7cb0dc719820ca Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Mon, 13 Jul 2026 13:28:52 +0000
Subject: [PATCH 03/13] Address review feedback on GeoTIFF fix

---
 murgtools/getdata/getDataFRF.py |  5 +++--
 tests/test_getDataFRF.py        | 16 ++++------------
 2 files changed, 7 insertions(+), 14 deletions(-)

diff --git a/murgtools/getdata/getDataFRF.py b/murgtools/getdata/getDataFRF.py
index fa6f7b6..9345ffb 100755
--- a/murgtools/getdata/getDataFRF.py
+++ b/murgtools/getdata/getDataFRF.py
@@ -3384,14 +3384,15 @@ def get_geotiff_extent(filepath, to_latlon=False):
         left = tiepoint[3]
         right = left + width * scale[0]
         top = tiepoint[4]
-        edge_y = top - height * scale[1]
-        bottom, top = sorted((top, edge_y))
+        y_edge = top - height * scale[1]
+        bottom, top = sorted((top, y_edge))
 
         if not to_latlon:
             return [min(left, right), max(left, right), bottom, top]
 
         geotiff_tags = page.geotiff_tags or {}
         epsg = geotiff_tags.get('ProjectedCSTypeGeoKey') or geotiff_tags.get('GeographicTypeGeoKey')
+        # tifffile may expose GeoTIFF keys as enums, so unwrap the numeric EPSG value when needed.
         if hasattr(epsg, 'value'):
             epsg = epsg.value
         if epsg is None:
diff --git a/tests/test_getDataFRF.py b/tests/test_getDataFRF.py
index 3bdbb5b..cef6307 100644
--- a/tests/test_getDataFRF.py
+++ b/tests/test_getDataFRF.py
@@ -1,12 +1,15 @@
 """Unit tests for getDataFRF module."""
 import datetime as DT
 import os
+import tempfile
 import numpy as np
+import tifffile
 import pytest
 import netCDF4 as nc
 from unittest.mock import MagicMock, patch, PropertyMock
+from pyproj import Transformer
 
-from murgtools.getdata.getDataFRF import gettime, removeDuplicatesFromDictionary
+from murgtools.getdata.getDataFRF import get_geotiff_extent, gettime, removeDuplicatesFromDictionary
 from murgtools.exceptions import InvalidGaugeError
 
 
@@ -284,11 +287,6 @@ class TestGetGeoTiffExtent:
 
     def test_returns_native_extent(self):
         """Test native GeoTIFF extent extraction."""
-        import tempfile
-        import tifffile
-
-        from murgtools.getdata.getDataFRF import get_geotiff_extent
-
         image = np.zeros((20, 10), dtype=np.uint8)
         tiepoint = [0.0, 0.0, 0.0, 900500.0, 276000.0, 0.0]
         scale = [10.0, 10.0, 0.0]
@@ -312,12 +310,6 @@ def test_returns_native_extent(self):
 
     def test_converts_projected_extent_to_latlon(self):
         """Test projected GeoTIFF extent conversion to lon/lat."""
-        import tempfile
-        import tifffile
-        from pyproj import Transformer
-
-        from murgtools.getdata.getDataFRF import get_geotiff_extent
-
         image = np.zeros((20, 10), dtype=np.uint8)
         geokey_dir = (1, 1, 0, 3, 1024, 0, 1, 1, 1025, 0, 1, 1, 3072, 0, 1, 32618)
         tiepoint = [0.0, 0.0, 0.0, 500000.0, 4000000.0, 0.0]

From d1c9fa252f2d40f35943cb5e24c9263b8609239f Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Mon, 13 Jul 2026 13:29:41 +0000
Subject: [PATCH 04/13] Polish GeoTIFF extent validation feedback

---
 murgtools/getdata/getDataFRF.py | 5 ++++-
 tests/test_getDataFRF.py        | 1 +
 2 files changed, 5 insertions(+), 1 deletion(-)

diff --git a/murgtools/getdata/getDataFRF.py b/murgtools/getdata/getDataFRF.py
index 9345ffb..022350e 100755
--- a/murgtools/getdata/getDataFRF.py
+++ b/murgtools/getdata/getDataFRF.py
@@ -3396,7 +3396,10 @@ def get_geotiff_extent(filepath, to_latlon=False):
         if hasattr(epsg, 'value'):
             epsg = epsg.value
         if epsg is None:
-            raise ValueError(f"GeoTIFF CRS metadata not found in {filepath}")
+            raise ValueError(
+                f"GeoTIFF CRS metadata (ProjectedCSTypeGeoKey or GeographicTypeGeoKey) "
+                f"not found in {filepath}"
+            )
 
         transformer = pyproj.Transformer.from_crs(f"EPSG:{epsg}", "EPSG:4326", always_xy=True)
         corners = [
diff --git a/tests/test_getDataFRF.py b/tests/test_getDataFRF.py
index cef6307..7536f1c 100644
--- a/tests/test_getDataFRF.py
+++ b/tests/test_getDataFRF.py
@@ -311,6 +311,7 @@ def test_returns_native_extent(self):
     def test_converts_projected_extent_to_latlon(self):
         """Test projected GeoTIFF extent conversion to lon/lat."""
         image = np.zeros((20, 10), dtype=np.uint8)
+        # GeoKey Directory: (version, revision, minor_rev, num_keys, key_id, location, count, value...)
         geokey_dir = (1, 1, 0, 3, 1024, 0, 1, 1, 1025, 0, 1, 1, 3072, 0, 1, 32618)
         tiepoint = [0.0, 0.0, 0.0, 500000.0, 4000000.0, 0.0]
         scale = [10.0, 10.0, 0.0]

From ff436b7e2656b971b0838654c18bc13393fe3e8c Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Mon, 13 Jul 2026 13:30:31 +0000
Subject: [PATCH 05/13] Finalize GeoTIFF extent helper cleanup

---
 murgtools/getdata/getDataFRF.py | 6 ++++--
 tests/test_getDataFRF.py        | 2 +-
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/murgtools/getdata/getDataFRF.py b/murgtools/getdata/getDataFRF.py
index 022350e..ebf233c 100755
--- a/murgtools/getdata/getDataFRF.py
+++ b/murgtools/getdata/getDataFRF.py
@@ -3388,10 +3388,12 @@ def get_geotiff_extent(filepath, to_latlon=False):
         bottom, top = sorted((top, y_edge))
 
         if not to_latlon:
-            return [min(left, right), max(left, right), bottom, top]
+            return [left, right, bottom, top]
 
         geotiff_tags = page.geotiff_tags or {}
-        epsg = geotiff_tags.get('ProjectedCSTypeGeoKey') or geotiff_tags.get('GeographicTypeGeoKey')
+        epsg = geotiff_tags.get('ProjectedCSTypeGeoKey')
+        if epsg is None:
+            epsg = geotiff_tags.get('GeographicTypeGeoKey')
         # tifffile may expose GeoTIFF keys as enums, so unwrap the numeric EPSG value when needed.
         if hasattr(epsg, 'value'):
             epsg = epsg.value
diff --git a/tests/test_getDataFRF.py b/tests/test_getDataFRF.py
index 7536f1c..45d9152 100644
--- a/tests/test_getDataFRF.py
+++ b/tests/test_getDataFRF.py
@@ -304,7 +304,7 @@ def test_returns_native_extent(self):
                 ]
             )
 
-            assert get_geotiff_extent(tmp_path) == [900500.0, 900600.0, 275800.0, 276000.0]
+            assert get_geotiff_extent(tmp_path) == pytest.approx([900500.0, 900600.0, 275800.0, 276000.0])
         finally:
             os.unlink(tmp_path)
 

From 0614dc80cab7adb246346aa027eb1d90f49217f9 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Mon, 13 Jul 2026 13:31:23 +0000
Subject: [PATCH 06/13] Clarify GeoTIFF CRS helper logic

---
 murgtools/getdata/getDataFRF.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/murgtools/getdata/getDataFRF.py b/murgtools/getdata/getDataFRF.py
index ebf233c..329ae83 100755
--- a/murgtools/getdata/getDataFRF.py
+++ b/murgtools/getdata/getDataFRF.py
@@ -3385,12 +3385,15 @@ def get_geotiff_extent(filepath, to_latlon=False):
         right = left + width * scale[0]
         top = tiepoint[4]
         y_edge = top - height * scale[1]
-        bottom, top = sorted((top, y_edge))
+        bottom = min(top, y_edge)
+        top = max(top, y_edge)
 
         if not to_latlon:
             return [left, right, bottom, top]
 
         geotiff_tags = page.geotiff_tags or {}
+        # Projected GeoTIFFs store their CRS in ProjectedCSTypeGeoKey, while
+        # geographic lon/lat GeoTIFFs use GeographicTypeGeoKey.
         epsg = geotiff_tags.get('ProjectedCSTypeGeoKey')
         if epsg is None:
             epsg = geotiff_tags.get('GeographicTypeGeoKey')

From b16e13ed77b70dd750642ea3947452ff17059f63 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Mon, 13 Jul 2026 13:32:10 +0000
Subject: [PATCH 07/13] Tidy GeoTIFF extent helper readability

---
 murgtools/getdata/getDataFRF.py |  8 ++++----
 tests/test_getDataFRF.py        | 10 ++++++----
 2 files changed, 10 insertions(+), 8 deletions(-)

diff --git a/murgtools/getdata/getDataFRF.py b/murgtools/getdata/getDataFRF.py
index 329ae83..f975bd4 100755
--- a/murgtools/getdata/getDataFRF.py
+++ b/murgtools/getdata/getDataFRF.py
@@ -3383,10 +3383,10 @@ def get_geotiff_extent(filepath, to_latlon=False):
         # Compute extent: [left, right, bottom, top]
         left = tiepoint[3]
         right = left + width * scale[0]
-        top = tiepoint[4]
-        y_edge = top - height * scale[1]
-        bottom = min(top, y_edge)
-        top = max(top, y_edge)
+        top_edge = tiepoint[4]
+        bottom_edge = top_edge - height * scale[1]
+        bottom = min(top_edge, bottom_edge)
+        top = max(top_edge, bottom_edge)
 
         if not to_latlon:
             return [left, right, bottom, top]
diff --git a/tests/test_getDataFRF.py b/tests/test_getDataFRF.py
index 45d9152..077bd3a 100644
--- a/tests/test_getDataFRF.py
+++ b/tests/test_getDataFRF.py
@@ -339,11 +339,13 @@ def test_converts_projected_extent_to_latlon(self):
                 transformer.transform(500100.0, 3999800.0),
                 transformer.transform(500100.0, 4000000.0),
             ]
+            lons = [lon for lon, _ in corners]
+            lats = [lat for _, lat in corners]
             expected = [
-                min(lon for lon, _ in corners),
-                max(lon for lon, _ in corners),
-                min(lat for _, lat in corners),
-                max(lat for _, lat in corners),
+                min(lons),
+                max(lons),
+                min(lats),
+                max(lats),
             ]
 
             assert extent == pytest.approx(expected)

From 793d97f59bd04b317efb395601c74bc11c57c1fb Mon Sep 17 00:00:00 2001
From: Spicer Bak <spicer.bak.frf@gmail.com>
Date: Tue, 14 Jul 2026 08:55:35 -0400
Subject: [PATCH 08/13] Fix satellite imagery coordinate alignment using UTM
 projection

The satellite imagery was offset ~784m east and ~182m north due to
using a linear lat/lon-to-pixel mapping. This assumed a constant
relationship between degrees and pixels across the scene, but UTM
projection doesn't have a linear relationship to lat/lon.

Fix by:
- Converting bbox corners from lat/lon to UTM coordinates
- Using GeoTIFF native tiepoint/scale for accurate pixel indices
- Converting actual cropped bounds back to lat/lon for extent

This ensures proper alignment between satellite imagery and other
geographic data like Argus orthophotos and gauge locations.

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
---
 murgtools/getdata/getOutsideData.py | 56 ++++++++++++++++++++++-------
 1 file changed, 44 insertions(+), 12 deletions(-)

diff --git a/murgtools/getdata/getOutsideData.py b/murgtools/getdata/getOutsideData.py
index d19c2ca..4c4daa2 100755
--- a/murgtools/getdata/getOutsideData.py
+++ b/murgtools/getdata/getOutsideData.py
@@ -441,19 +441,50 @@ def getSatelliteImagery(corners, filename=None, collection='sentinel-2-l2a',
     if collection == 'naip' and image.ndim == 3 and image.shape[-1] == 4:
         image = image[:, :, :3]  # Keep only RGB, drop NIR
 
-    # 7. Crop to bbox using actual image extent from GeoTIFF
+    # 7. Crop to bbox using UTM coordinates for accuracy
+    # The GeoTIFF stores UTM coordinates in tiepoint/scale, so we must convert
+    # the lat/lon bbox to UTM to get correct pixel indices.
     h, w = image.shape[:2]
-    px_per_deg_x = w / (scene_bbox[2] - scene_bbox[0])
-    px_per_deg_y = h / (scene_bbox[3] - scene_bbox[1])
 
-    x1 = int((bbox[0] - scene_bbox[0]) * px_per_deg_x)
-    x2 = int((bbox[2] - scene_bbox[0]) * px_per_deg_x)
-    y1 = int((scene_bbox[3] - bbox[3]) * px_per_deg_y)
-    y2 = int((scene_bbox[3] - bbox[1]) * px_per_deg_y)
+    # Create transformer from lat/lon to the image's UTM CRS
+    to_utm = Transformer.from_crs('EPSG:4326', f'EPSG:{epsg}', always_xy=True)
+
+    # Convert bbox corners to UTM
+    bbox_utm_west, _ = to_utm.transform(bbox[0], bbox[1])
+    bbox_utm_east, _ = to_utm.transform(bbox[2], bbox[1])
+    _, bbox_utm_south = to_utm.transform(bbox[0], bbox[1])
+    _, bbox_utm_north = to_utm.transform(bbox[0], bbox[3])
+
+    # Use GeoTIFF origin and scale for pixel calculation
+    # tiepoint[3] = UTM easting of pixel (0,0)
+    # tiepoint[4] = UTM northing of pixel (0,0) (top of image)
+    # scale[0] = meters per pixel in X
+    # scale[1] = meters per pixel in Y
+    utm_origin_x = tiepoint[3]
+    utm_origin_y = tiepoint[4]
+    scale_x = scale[0]
+    scale_y = scale[1]
+
+    # Compute pixel indices from UTM coordinates
+    x1 = int((bbox_utm_west - utm_origin_x) / scale_x)
+    x2 = int((bbox_utm_east - utm_origin_x) / scale_x)
+    y1 = int((utm_origin_y - bbox_utm_north) / scale_y)
+    y2 = int((utm_origin_y - bbox_utm_south) / scale_y)
 
     x1, x2 = max(0, x1), min(w, x2)
     y1, y2 = max(0, y1), min(h, y2)
 
+    # Compute actual UTM bounds of cropped region (for accurate extent)
+    actual_utm_west = utm_origin_x + x1 * scale_x
+    actual_utm_east = utm_origin_x + x2 * scale_x
+    actual_utm_north = utm_origin_y - y1 * scale_y
+    actual_utm_south = utm_origin_y - y2 * scale_y
+
+    # Convert actual UTM bounds back to lat/lon
+    from_utm = Transformer.from_crs(f'EPSG:{epsg}', 'EPSG:4326', always_xy=True)
+    actual_west, actual_south = from_utm.transform(actual_utm_west, actual_utm_south)
+    actual_east, actual_north = from_utm.transform(actual_utm_east, actual_utm_north)
+
     image = image[y1:y2, x1:x2]
 
     # 8. Ensure uint8 RGB format
@@ -461,11 +492,12 @@ def getSatelliteImagery(corners, filename=None, collection='sentinel-2-l2a',
         image = np.clip(image / image.max() * 255, 0, 255).astype(np.uint8)
 
     # 9. Compute georeferencing before rotation
-    deg_per_px_x = (bbox[2] - bbox[0]) / image.shape[1]
-    deg_per_px_y = (bbox[3] - bbox[1]) / image.shape[0]
+    # Use actual bounds for accurate deg_per_px
+    deg_per_px_x = (actual_east - actual_west) / image.shape[1]
+    deg_per_px_y = (actual_north - actual_south) / image.shape[0]
 
     meters_per_deg = 111320 * np.cos(np.radians(np.mean(lats)))
-    resolution_m = deg_per_px_x * meters_per_deg
+    resolution_m = scale_x  # Use actual GeoTIFF scale
 
     # 10. Rotate image to match AOI orientation
     if abs(rotation_angle) > 0.1:
@@ -474,8 +506,8 @@ def getSatelliteImagery(corners, filename=None, collection='sentinel-2-l2a',
     post_rot_shape = image.shape[:2]
 
     # 11. Compute georeferencing for rotated image
-    center_lon = (bbox[0] + bbox[2]) / 2
-    center_lat = (bbox[1] + bbox[3]) / 2
+    center_lon = (actual_west + actual_east) / 2
+    center_lat = (actual_south + actual_north) / 2
 
     rot_rad = np.radians(rotation_angle)
     cos_r, sin_r = np.cos(rot_rad), np.sin(rot_rad)

From cf68f3f5a41cd9fc71099ff8b2617fce4e31b562 Mon Sep 17 00:00:00 2001
From: Spicer Bak <spicer.bak.frf@gmail.com>
Date: Tue, 14 Jul 2026 09:31:09 -0400
Subject: [PATCH 09/13] Address PR review feedback on coordinate handling

Review feedback addressed:

1. Fallback path undefined variables: Added explicit initialization of
   tiepoint, scale, epsg to None with has_geotiff_metadata flag to track
   whether GeoTIFF metadata was successfully extracted.

2. BBox UTM conversion: Now converts all 4 bbox corners to UTM and takes
   min/max for proper bounds, handling rotated/skewed AOIs correctly.

3. int() truncation: Use math.floor for x1/y1 and math.ceil for x2/y2
   to include the full requested bbox before clamping.

4. resolution_m in non-UTM: Added fallback path that computes resolution
   from deg_per_px when GeoTIFF metadata is unavailable.

5. GeoKey parsing: Support both string key names and integer key IDs
   (3072/2048), and handle string EPSG values like "EPSG:32618".

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
---
 murgtools/getdata/getDataFRF.py     |  25 ++++--
 murgtools/getdata/getOutsideData.py | 121 +++++++++++++++++-----------
 2 files changed, 95 insertions(+), 51 deletions(-)

diff --git a/murgtools/getdata/getDataFRF.py b/murgtools/getdata/getDataFRF.py
index f975bd4..320d570 100755
--- a/murgtools/getdata/getDataFRF.py
+++ b/murgtools/getdata/getDataFRF.py
@@ -3392,14 +3392,27 @@ def get_geotiff_extent(filepath, to_latlon=False):
             return [left, right, bottom, top]
 
         geotiff_tags = page.geotiff_tags or {}
-        # Projected GeoTIFFs store their CRS in ProjectedCSTypeGeoKey, while
-        # geographic lon/lat GeoTIFFs use GeographicTypeGeoKey.
-        epsg = geotiff_tags.get('ProjectedCSTypeGeoKey')
-        if epsg is None:
-            epsg = geotiff_tags.get('GeographicTypeGeoKey')
-        # tifffile may expose GeoTIFF keys as enums, so unwrap the numeric EPSG value when needed.
+        # Projected GeoTIFFs store their CRS in ProjectedCSTypeGeoKey (3072), while
+        # geographic lon/lat GeoTIFFs use GeographicTypeGeoKey (2048).
+        # Support both string key names and integer key IDs.
+        epsg = (geotiff_tags.get('ProjectedCSTypeGeoKey') or
+                geotiff_tags.get(3072) or
+                geotiff_tags.get('GeographicTypeGeoKey') or
+                geotiff_tags.get(2048))
+
+        # tifffile may expose GeoTIFF keys as enums, so unwrap the numeric value.
         if hasattr(epsg, 'value'):
             epsg = epsg.value
+
+        # Handle string EPSG values like "EPSG:32618" or "32618"
+        if isinstance(epsg, str):
+            import re
+            match = re.search(r'(\d+)', epsg)
+            if match:
+                epsg = int(match.group(1))
+            else:
+                epsg = None
+
         if epsg is None:
             raise ValueError(
                 f"GeoTIFF CRS metadata (ProjectedCSTypeGeoKey or GeographicTypeGeoKey) "
diff --git a/murgtools/getdata/getOutsideData.py b/murgtools/getdata/getOutsideData.py
index 4c4daa2..5bd6667 100755
--- a/murgtools/getdata/getOutsideData.py
+++ b/murgtools/getdata/getOutsideData.py
@@ -391,6 +391,12 @@ def getSatelliteImagery(corners, filename=None, collection='sentinel-2-l2a',
         for chunk in resp.iter_content(chunk_size=8192):
             tmp.write(chunk)
 
+    # Track whether we have GeoTIFF metadata for accurate UTM-based cropping
+    has_geotiff_metadata = False
+    tiepoint = None
+    scale = None
+    epsg = None
+
     try:
         image = tifffile.imread(tmp_path)
 
@@ -421,6 +427,8 @@ def getSatelliteImagery(corners, filename=None, collection='sentinel-2-l2a',
                         epsg = 32600 + zone if hemisphere == 'N' else 32700 + zone
                     else:
                         epsg = 32618
+                else:
+                    epsg = 32618  # Default to UTM zone 18N if no projection string
 
                 # Convert UTM corners to lat/lon
                 transformer = Transformer.from_crs(f'EPSG:{epsg}', 'EPSG:4326', always_xy=True)
@@ -430,6 +438,7 @@ def getSatelliteImagery(corners, filename=None, collection='sentinel-2-l2a',
 
                 # scene_bbox in lat/lon: [west, south, east, north]
                 scene_bbox = [tl_lon, br_lat, br_lon, tl_lat]
+                has_geotiff_metadata = True
         except KeyError:
             # Missing GeoTIFF georeferencing tags; fall back to STAC bbox if available.
             if isinstance(item, dict) and 'bbox' in item:
@@ -441,49 +450,75 @@ def getSatelliteImagery(corners, filename=None, collection='sentinel-2-l2a',
     if collection == 'naip' and image.ndim == 3 and image.shape[-1] == 4:
         image = image[:, :, :3]  # Keep only RGB, drop NIR
 
-    # 7. Crop to bbox using UTM coordinates for accuracy
-    # The GeoTIFF stores UTM coordinates in tiepoint/scale, so we must convert
-    # the lat/lon bbox to UTM to get correct pixel indices.
+    # 7. Crop to bbox - use UTM coordinates if GeoTIFF metadata available
     h, w = image.shape[:2]
+    import math
+
+    if has_geotiff_metadata:
+        # Use UTM coordinates for accurate pixel calculation
+        to_utm = Transformer.from_crs('EPSG:4326', f'EPSG:{epsg}', always_xy=True)
+
+        # Convert all 4 bbox corners to UTM and take min/max for proper bounds
+        corners_ll = [
+            (bbox[0], bbox[1]),  # SW
+            (bbox[0], bbox[3]),  # NW
+            (bbox[2], bbox[1]),  # SE
+            (bbox[2], bbox[3]),  # NE
+        ]
+        corners_utm = [to_utm.transform(lon, lat) for lon, lat in corners_ll]
+        bbox_utm_west = min(c[0] for c in corners_utm)
+        bbox_utm_east = max(c[0] for c in corners_utm)
+        bbox_utm_south = min(c[1] for c in corners_utm)
+        bbox_utm_north = max(c[1] for c in corners_utm)
+
+        # Use GeoTIFF origin and scale for pixel calculation
+        utm_origin_x = tiepoint[3]
+        utm_origin_y = tiepoint[4]
+        scale_x = scale[0]
+        scale_y = scale[1]
+
+        # Compute pixel indices using floor/ceil to include full bbox
+        x1 = math.floor((bbox_utm_west - utm_origin_x) / scale_x)
+        x2 = math.ceil((bbox_utm_east - utm_origin_x) / scale_x)
+        y1 = math.floor((utm_origin_y - bbox_utm_north) / scale_y)
+        y2 = math.ceil((utm_origin_y - bbox_utm_south) / scale_y)
+
+        x1, x2 = max(0, x1), min(w, x2)
+        y1, y2 = max(0, y1), min(h, y2)
+
+        # Compute actual UTM bounds of cropped region
+        actual_utm_west = utm_origin_x + x1 * scale_x
+        actual_utm_east = utm_origin_x + x2 * scale_x
+        actual_utm_north = utm_origin_y - y1 * scale_y
+        actual_utm_south = utm_origin_y - y2 * scale_y
+
+        # Convert actual UTM bounds back to lat/lon
+        from_utm = Transformer.from_crs(f'EPSG:{epsg}', 'EPSG:4326', always_xy=True)
+        actual_west, actual_south = from_utm.transform(actual_utm_west, actual_utm_south)
+        actual_east, actual_north = from_utm.transform(actual_utm_east, actual_utm_north)
+
+        resolution_m = scale_x  # GeoTIFF scale is in meters for UTM
+    else:
+        # Fallback: use linear lat/lon mapping (less accurate but works without GeoTIFF metadata)
+        px_per_deg_x = w / (scene_bbox[2] - scene_bbox[0])
+        px_per_deg_y = h / (scene_bbox[3] - scene_bbox[1])
+
+        x1 = math.floor((bbox[0] - scene_bbox[0]) * px_per_deg_x)
+        x2 = math.ceil((bbox[2] - scene_bbox[0]) * px_per_deg_x)
+        y1 = math.floor((scene_bbox[3] - bbox[3]) * px_per_deg_y)
+        y2 = math.ceil((scene_bbox[3] - bbox[1]) * px_per_deg_y)
 
-    # Create transformer from lat/lon to the image's UTM CRS
-    to_utm = Transformer.from_crs('EPSG:4326', f'EPSG:{epsg}', always_xy=True)
-
-    # Convert bbox corners to UTM
-    bbox_utm_west, _ = to_utm.transform(bbox[0], bbox[1])
-    bbox_utm_east, _ = to_utm.transform(bbox[2], bbox[1])
-    _, bbox_utm_south = to_utm.transform(bbox[0], bbox[1])
-    _, bbox_utm_north = to_utm.transform(bbox[0], bbox[3])
-
-    # Use GeoTIFF origin and scale for pixel calculation
-    # tiepoint[3] = UTM easting of pixel (0,0)
-    # tiepoint[4] = UTM northing of pixel (0,0) (top of image)
-    # scale[0] = meters per pixel in X
-    # scale[1] = meters per pixel in Y
-    utm_origin_x = tiepoint[3]
-    utm_origin_y = tiepoint[4]
-    scale_x = scale[0]
-    scale_y = scale[1]
-
-    # Compute pixel indices from UTM coordinates
-    x1 = int((bbox_utm_west - utm_origin_x) / scale_x)
-    x2 = int((bbox_utm_east - utm_origin_x) / scale_x)
-    y1 = int((utm_origin_y - bbox_utm_north) / scale_y)
-    y2 = int((utm_origin_y - bbox_utm_south) / scale_y)
-
-    x1, x2 = max(0, x1), min(w, x2)
-    y1, y2 = max(0, y1), min(h, y2)
-
-    # Compute actual UTM bounds of cropped region (for accurate extent)
-    actual_utm_west = utm_origin_x + x1 * scale_x
-    actual_utm_east = utm_origin_x + x2 * scale_x
-    actual_utm_north = utm_origin_y - y1 * scale_y
-    actual_utm_south = utm_origin_y - y2 * scale_y
-
-    # Convert actual UTM bounds back to lat/lon
-    from_utm = Transformer.from_crs(f'EPSG:{epsg}', 'EPSG:4326', always_xy=True)
-    actual_west, actual_south = from_utm.transform(actual_utm_west, actual_utm_south)
-    actual_east, actual_north = from_utm.transform(actual_utm_east, actual_utm_north)
+        x1, x2 = max(0, x1), min(w, x2)
+        y1, y2 = max(0, y1), min(h, y2)
+
+        # Compute actual bounds from pixel indices
+        actual_west = scene_bbox[0] + x1 / px_per_deg_x
+        actual_east = scene_bbox[0] + x2 / px_per_deg_x
+        actual_north = scene_bbox[3] - y1 / px_per_deg_y
+        actual_south = scene_bbox[3] - y2 / px_per_deg_y
+
+        meters_per_deg = 111320 * np.cos(np.radians(np.mean(lats)))
+        resolution_m = (actual_east - actual_west) / (x2 - x1) * meters_per_deg
 
     image = image[y1:y2, x1:x2]
 
@@ -492,13 +527,9 @@ def getSatelliteImagery(corners, filename=None, collection='sentinel-2-l2a',
         image = np.clip(image / image.max() * 255, 0, 255).astype(np.uint8)
 
     # 9. Compute georeferencing before rotation
-    # Use actual bounds for accurate deg_per_px
     deg_per_px_x = (actual_east - actual_west) / image.shape[1]
     deg_per_px_y = (actual_north - actual_south) / image.shape[0]
 
-    meters_per_deg = 111320 * np.cos(np.radians(np.mean(lats)))
-    resolution_m = scale_x  # Use actual GeoTIFF scale
-
     # 10. Rotate image to match AOI orientation
     if abs(rotation_angle) > 0.1:
         image = scipy_rotate(image, -rotation_angle, reshape=True, order=1,

From c82739e9f2d94556c26d1d848ac9d3812185b7ee Mon Sep 17 00:00:00 2001
From: Spicer Bak <spicer.bak.frf@gmail.com>
Date: Tue, 14 Jul 2026 09:48:10 -0400
Subject: [PATCH 10/13] Address additional PR review feedback

1. Handle negative scale values in get_geotiff_extent:
   - Normalize both X and Y edges with min/max to handle negative
     scale_x or scale_y correctly

2. Add test coverage for negative scale values:
   - test_handles_negative_scale_y: Verify extent with negative scale_y
   - test_handles_negative_scale_x: Verify extent with negative scale_x

3. Prevent division by zero on empty crop:
   - Check for x2 <= x1 or y2 <= y1 after clamping crop indices
   - Return None if bbox is completely outside the scene
   - Applies to both UTM and fallback code paths

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
---
 murgtools/getdata/getDataFRF.py     | 16 +++++----
 murgtools/getdata/getOutsideData.py |  8 +++++
 tests/test_getDataFRF.py            | 55 +++++++++++++++++++++++++++++
 3 files changed, 73 insertions(+), 6 deletions(-)

diff --git a/murgtools/getdata/getDataFRF.py b/murgtools/getdata/getDataFRF.py
index 320d570..68ece6a 100755
--- a/murgtools/getdata/getDataFRF.py
+++ b/murgtools/getdata/getDataFRF.py
@@ -3381,12 +3381,16 @@ def get_geotiff_extent(filepath, to_latlon=False):
         scale = tags[33550].value     # (scaleX, scaleY, scaleZ)
         height, width = page.shape[:2]
         # Compute extent: [left, right, bottom, top]
-        left = tiepoint[3]
-        right = left + width * scale[0]
-        top_edge = tiepoint[4]
-        bottom_edge = top_edge - height * scale[1]
-        bottom = min(top_edge, bottom_edge)
-        top = max(top_edge, bottom_edge)
+        # Handle both positive and negative scale values by computing both
+        # edges and normalizing with min/max
+        x_edge1 = tiepoint[3]
+        x_edge2 = tiepoint[3] + width * scale[0]
+        y_edge1 = tiepoint[4]
+        y_edge2 = tiepoint[4] - height * scale[1]
+        left = min(x_edge1, x_edge2)
+        right = max(x_edge1, x_edge2)
+        bottom = min(y_edge1, y_edge2)
+        top = max(y_edge1, y_edge2)
 
         if not to_latlon:
             return [left, right, bottom, top]
diff --git a/murgtools/getdata/getOutsideData.py b/murgtools/getdata/getOutsideData.py
index 5bd6667..2ffcc2a 100755
--- a/murgtools/getdata/getOutsideData.py
+++ b/murgtools/getdata/getOutsideData.py
@@ -486,6 +486,10 @@ def getSatelliteImagery(corners, filename=None, collection='sentinel-2-l2a',
         x1, x2 = max(0, x1), min(w, x2)
         y1, y2 = max(0, y1), min(h, y2)
 
+        # Check for empty crop region (bbox outside scene)
+        if x2 <= x1 or y2 <= y1:
+            return None
+
         # Compute actual UTM bounds of cropped region
         actual_utm_west = utm_origin_x + x1 * scale_x
         actual_utm_east = utm_origin_x + x2 * scale_x
@@ -511,6 +515,10 @@ def getSatelliteImagery(corners, filename=None, collection='sentinel-2-l2a',
         x1, x2 = max(0, x1), min(w, x2)
         y1, y2 = max(0, y1), min(h, y2)
 
+        # Check for empty crop region (bbox outside scene)
+        if x2 <= x1 or y2 <= y1:
+            return None
+
         # Compute actual bounds from pixel indices
         actual_west = scene_bbox[0] + x1 / px_per_deg_x
         actual_east = scene_bbox[0] + x2 / px_per_deg_x
diff --git a/tests/test_getDataFRF.py b/tests/test_getDataFRF.py
index 077bd3a..a6cf5f9 100644
--- a/tests/test_getDataFRF.py
+++ b/tests/test_getDataFRF.py
@@ -352,6 +352,61 @@ def test_converts_projected_extent_to_latlon(self):
         finally:
             os.unlink(tmp_path)
 
+    def test_handles_negative_scale_y(self):
+        """Test GeoTIFF extent with negative scale_y (row 0 at bottom)."""
+        image = np.zeros((200, 100), dtype=np.uint8)
+        # Tiepoint at (0,0) pixel -> (900500, 275800) with NEGATIVE scale_y
+        # This means Y increases as row number increases (row 0 at bottom)
+        tiepoint = [0.0, 0.0, 0.0, 900500.0, 275800.0, 0.0]
+        scale = [1.0, -1.0, 0.0]  # Negative scale_y
+
+        with tempfile.NamedTemporaryFile(suffix='.tif', delete=False) as tmp:
+            tmp_path = tmp.name
+
+        try:
+            tifffile.imwrite(
+                tmp_path,
+                image,
+                extratags=[
+                    (33922, 'd', 6, tiepoint, True),
+                    (33550, 'd', 3, scale, True),
+                ]
+            )
+
+            extent = get_geotiff_extent(tmp_path)
+            # With negative scale_y, Y goes from 275800 to 275800 - 200*(-1) = 276000
+            # Extent should be normalized: [left, right, bottom, top]
+            assert extent == pytest.approx([900500.0, 900600.0, 275800.0, 276000.0])
+        finally:
+            os.unlink(tmp_path)
+
+    def test_handles_negative_scale_x(self):
+        """Test GeoTIFF extent with negative scale_x (mirrored horizontally)."""
+        image = np.zeros((200, 100), dtype=np.uint8)
+        # Tiepoint with NEGATIVE scale_x (X decreases as column increases)
+        tiepoint = [0.0, 0.0, 0.0, 900600.0, 276000.0, 0.0]
+        scale = [-1.0, 1.0, 0.0]  # Negative scale_x
+
+        with tempfile.NamedTemporaryFile(suffix='.tif', delete=False) as tmp:
+            tmp_path = tmp.name
+
+        try:
+            tifffile.imwrite(
+                tmp_path,
+                image,
+                extratags=[
+                    (33922, 'd', 6, tiepoint, True),
+                    (33550, 'd', 3, scale, True),
+                ]
+            )
+
+            extent = get_geotiff_extent(tmp_path)
+            # With negative scale_x, X goes from 900600 to 900600 + 100*(-1) = 900500
+            # Extent should be normalized: [left, right, bottom, top]
+            assert extent == pytest.approx([900500.0, 900600.0, 275800.0, 276000.0])
+        finally:
+            os.unlink(tmp_path)
+
 
 class TestThreadGetArgusImagery:
     """Tests for the threadGetArgusImagery function."""

From ef992b2ab656a63bb447f03771e9f766fcc2522c Mon Sep 17 00:00:00 2001
From: Spicer Bak <spicer.bak.frf@gmail.com>
Date: Tue, 14 Jul 2026 10:03:44 -0400
Subject: [PATCH 11/13] Handle negative scale and improve EPSG detection

1. Negative scale handling in getSatelliteImagery:
   - Normalize UTM bounds with min/max to handle negative scale
   - Normalize pixel indices with min/max before clamping
   - Normalize actual UTM bounds after cropping
   - Use abs(scale_x) for resolution_m

2. Improved EPSG detection:
   - Check GeoKeys (ProjectedCSTypeGeoKey/3072) first before parsing
     projection string, as GeoKeys are the authoritative source
   - Only default to EPSG:32618 as last resort when no CRS found

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
---
 murgtools/getdata/getOutsideData.py | 89 +++++++++++++++++++----------
 1 file changed, 60 insertions(+), 29 deletions(-)

diff --git a/murgtools/getdata/getOutsideData.py b/murgtools/getdata/getOutsideData.py
index 2ffcc2a..a528191 100755
--- a/murgtools/getdata/getOutsideData.py
+++ b/murgtools/getdata/getOutsideData.py
@@ -409,26 +409,44 @@ def getSatelliteImagery(corners, filename=None, collection='sentinel-2-l2a',
                 scale = tags[33550].value     # (scaleX, scaleY, scaleZ)
                 img_h, img_w = tif.pages[0].shape[:2]
 
-                # UTM bounds
-                utm_left = tiepoint[3]
-                utm_right = tiepoint[3] + img_w * scale[0]
-                utm_top = tiepoint[4]
-                utm_bottom = tiepoint[4] - img_h * scale[1]
-
-                # Determine UTM zone from GeoKey (tag 34737 contains projection info)
-                proj_str = tags.get(34737, None)
-                if proj_str:
-                    proj_str = proj_str.value
-                    # Parse UTM zone from string like "WGS 84 / UTM zone 18N"
-                    match = re.search(r'UTM zone (\d+)([NS])', str(proj_str))
-                    if match:
-                        zone = int(match.group(1))
-                        hemisphere = match.group(2)
-                        epsg = 32600 + zone if hemisphere == 'N' else 32700 + zone
-                    else:
-                        epsg = 32618
-                else:
-                    epsg = 32618  # Default to UTM zone 18N if no projection string
+                # UTM bounds - normalize with min/max to handle negative scale
+                x_edge1 = tiepoint[3]
+                x_edge2 = tiepoint[3] + img_w * scale[0]
+                y_edge1 = tiepoint[4]
+                y_edge2 = tiepoint[4] - img_h * scale[1]
+                utm_left = min(x_edge1, x_edge2)
+                utm_right = max(x_edge1, x_edge2)
+                utm_bottom = min(y_edge1, y_edge2)
+                utm_top = max(y_edge1, y_edge2)
+
+                # Determine EPSG from GeoKeys first, then fall back to parsing projection string
+                page = tif.pages[0]
+                geotiff_tags = page.geotiff_tags or {}
+                epsg = (geotiff_tags.get('ProjectedCSTypeGeoKey') or
+                        geotiff_tags.get(3072) or
+                        geotiff_tags.get('GeographicTypeGeoKey') or
+                        geotiff_tags.get(2048))
+                if hasattr(epsg, 'value'):
+                    epsg = epsg.value
+                if isinstance(epsg, str):
+                    match = re.search(r'(\d+)', epsg)
+                    epsg = int(match.group(1)) if match else None
+
+                # Fall back to parsing projection string if GeoKeys not found
+                if epsg is None:
+                    proj_str = tags.get(34737, None)
+                    if proj_str:
+                        proj_str = proj_str.value
+                        # Parse UTM zone from string like "WGS 84 / UTM zone 18N"
+                        match = re.search(r'UTM zone (\d+)([NS])', str(proj_str))
+                        if match:
+                            zone = int(match.group(1))
+                            hemisphere = match.group(2)
+                            epsg = 32600 + zone if hemisphere == 'N' else 32700 + zone
+
+                # Default to UTM zone 18N only as last resort
+                if epsg is None:
+                    epsg = 32618
 
                 # Convert UTM corners to lat/lon
                 transformer = Transformer.from_crs(f'EPSG:{epsg}', 'EPSG:4326', always_xy=True)
@@ -478,10 +496,17 @@ def getSatelliteImagery(corners, filename=None, collection='sentinel-2-l2a',
         scale_y = scale[1]
 
         # Compute pixel indices using floor/ceil to include full bbox
-        x1 = math.floor((bbox_utm_west - utm_origin_x) / scale_x)
-        x2 = math.ceil((bbox_utm_east - utm_origin_x) / scale_x)
-        y1 = math.floor((utm_origin_y - bbox_utm_north) / scale_y)
-        y2 = math.ceil((utm_origin_y - bbox_utm_south) / scale_y)
+        # Handle negative scale by computing both edges and normalizing
+        px_west = (bbox_utm_west - utm_origin_x) / scale_x
+        px_east = (bbox_utm_east - utm_origin_x) / scale_x
+        px_north = (utm_origin_y - bbox_utm_north) / scale_y
+        px_south = (utm_origin_y - bbox_utm_south) / scale_y
+
+        # Normalize pixel indices (handle negative scale)
+        x1 = math.floor(min(px_west, px_east))
+        x2 = math.ceil(max(px_west, px_east))
+        y1 = math.floor(min(px_north, px_south))
+        y2 = math.ceil(max(px_north, px_south))
 
         x1, x2 = max(0, x1), min(w, x2)
         y1, y2 = max(0, y1), min(h, y2)
@@ -491,17 +516,23 @@ def getSatelliteImagery(corners, filename=None, collection='sentinel-2-l2a',
             return None
 
         # Compute actual UTM bounds of cropped region
-        actual_utm_west = utm_origin_x + x1 * scale_x
-        actual_utm_east = utm_origin_x + x2 * scale_x
-        actual_utm_north = utm_origin_y - y1 * scale_y
-        actual_utm_south = utm_origin_y - y2 * scale_y
+        # Use the normalized image coordinates (utm_left/right/bottom/top from earlier)
+        # to ensure consistency regardless of scale sign
+        actual_utm_x1 = utm_origin_x + x1 * scale_x
+        actual_utm_x2 = utm_origin_x + x2 * scale_x
+        actual_utm_y1 = utm_origin_y - y1 * scale_y
+        actual_utm_y2 = utm_origin_y - y2 * scale_y
+        actual_utm_west = min(actual_utm_x1, actual_utm_x2)
+        actual_utm_east = max(actual_utm_x1, actual_utm_x2)
+        actual_utm_south = min(actual_utm_y1, actual_utm_y2)
+        actual_utm_north = max(actual_utm_y1, actual_utm_y2)
 
         # Convert actual UTM bounds back to lat/lon
         from_utm = Transformer.from_crs(f'EPSG:{epsg}', 'EPSG:4326', always_xy=True)
         actual_west, actual_south = from_utm.transform(actual_utm_west, actual_utm_south)
         actual_east, actual_north = from_utm.transform(actual_utm_east, actual_utm_north)
 
-        resolution_m = scale_x  # GeoTIFF scale is in meters for UTM
+        resolution_m = abs(scale_x)  # GeoTIFF scale magnitude in meters for UTM
     else:
         # Fallback: use linear lat/lon mapping (less accurate but works without GeoTIFF metadata)
         px_per_deg_x = w / (scene_bbox[2] - scene_bbox[0])

From e192a2a7ce4c3a75416332b67a671a58df6467fd Mon Sep 17 00:00:00 2001
From: Spicer Bak <spicer.bak.frf@gmail.com>
Date: Tue, 14 Jul 2026 10:43:43 -0400
Subject: [PATCH 12/13] Fix resolution_m to handle non-meter CRS units

Check the CRS unit using pyproj and convert appropriately:
- metre/meter/m: use scale directly
- degree/degrees: convert using meters_per_deg at image center
- foot/US survey foot/ft: convert using 0.3048 m/ft
- unknown: default to assuming meters

This ensures resolution_m is always in meters regardless of the
GeoTIFF CRS (UTM, geographic, or US survey feet).

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
---
 murgtools/getdata/getOutsideData.py | 19 ++++++++++++++++++-
 1 file changed, 18 insertions(+), 1 deletion(-)

diff --git a/murgtools/getdata/getOutsideData.py b/murgtools/getdata/getOutsideData.py
index a528191..6d54b57 100755
--- a/murgtools/getdata/getOutsideData.py
+++ b/murgtools/getdata/getOutsideData.py
@@ -532,7 +532,24 @@ def getSatelliteImagery(corners, filename=None, collection='sentinel-2-l2a',
         actual_west, actual_south = from_utm.transform(actual_utm_west, actual_utm_south)
         actual_east, actual_north = from_utm.transform(actual_utm_east, actual_utm_north)
 
-        resolution_m = abs(scale_x)  # GeoTIFF scale magnitude in meters for UTM
+        # Compute resolution in meters, accounting for CRS units
+        # UTM (326xx/327xx) uses meters, geographic CRS uses degrees
+        from pyproj import CRS
+        crs = CRS.from_epsg(epsg)
+        crs_unit = crs.axis_info[0].unit_name if crs.axis_info else 'metre'
+        if crs_unit in ('metre', 'meter', 'm'):
+            resolution_m = abs(scale_x)
+        elif crs_unit in ('degree', 'degrees'):
+            # Convert degrees to meters at image center latitude
+            center_lat = (actual_north + actual_south) / 2
+            meters_per_deg = 111320 * np.cos(np.radians(center_lat))
+            resolution_m = abs(scale_x) * meters_per_deg
+        elif crs_unit in ('foot', 'US survey foot', 'ft'):
+            # Convert feet to meters
+            resolution_m = abs(scale_x) * 0.3048
+        else:
+            # Default: assume meters for unknown units
+            resolution_m = abs(scale_x)
     else:
         # Fallback: use linear lat/lon mapping (less accurate but works without GeoTIFF metadata)
         px_per_deg_x = w / (scene_bbox[2] - scene_bbox[0])

From 8284fdfe66a2791653ea1a80419b222dec3ae600 Mon Sep 17 00:00:00 2001
From: Spicer Bak <spicer.bak.frf@gmail.com>
Date: Tue, 14 Jul 2026 11:42:13 -0400
Subject: [PATCH 13/13] Handle non-zero tiepoint pixel location in GeoTIFF

GeoTIFF ModelTiepointTag can reference any raster pixel (i,j), not just
(0,0). The code was assuming tiepoint[0]=tiepoint[1]=0, which caused
incorrect extent/bounds calculation for GeoTIFFs with non-zero tiepoint.

Fix by computing model-space origin for pixel (0,0):
  origin_x = tiepoint[3] - tiepoint[0] * scale[0]
  origin_y = tiepoint[4] + tiepoint[1] * scale[1]

Then use origin_x/origin_y instead of tiepoint[3]/tiepoint[4] for:
- Extent calculation in get_geotiff_extent()
- UTM bounds calculation in getSatelliteImagery()
- Pixel index calculation for cropping

Added test_handles_nonzero_tiepoint_pixel to verify correct handling.

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
---
 murgtools/getdata/getDataFRF.py     | 15 ++++++++---
 murgtools/getdata/getOutsideData.py | 39 ++++++++++++++++-------------
 tests/test_getDataFRF.py            | 30 ++++++++++++++++++++++
 3 files changed, 63 insertions(+), 21 deletions(-)

diff --git a/murgtools/getdata/getDataFRF.py b/murgtools/getdata/getDataFRF.py
index 68ece6a..b24b004 100755
--- a/murgtools/getdata/getDataFRF.py
+++ b/murgtools/getdata/getDataFRF.py
@@ -3380,13 +3380,20 @@ def get_geotiff_extent(filepath, to_latlon=False):
         tiepoint = tags[33922].value  # (i, j, k, x, y, z)
         scale = tags[33550].value     # (scaleX, scaleY, scaleZ)
         height, width = page.shape[:2]
+
+        # Compute model-space origin for pixel (0,0)
+        # Tiepoint specifies model coords (x,y) for raster pixel (i,j)
+        # origin_x = x - i * scaleX, origin_y = y + j * scaleY
+        origin_x = tiepoint[3] - tiepoint[0] * scale[0]
+        origin_y = tiepoint[4] + tiepoint[1] * scale[1]
+
         # Compute extent: [left, right, bottom, top]
         # Handle both positive and negative scale values by computing both
         # edges and normalizing with min/max
-        x_edge1 = tiepoint[3]
-        x_edge2 = tiepoint[3] + width * scale[0]
-        y_edge1 = tiepoint[4]
-        y_edge2 = tiepoint[4] - height * scale[1]
+        x_edge1 = origin_x
+        x_edge2 = origin_x + width * scale[0]
+        y_edge1 = origin_y
+        y_edge2 = origin_y - height * scale[1]
         left = min(x_edge1, x_edge2)
         right = max(x_edge1, x_edge2)
         bottom = min(y_edge1, y_edge2)
diff --git a/murgtools/getdata/getOutsideData.py b/murgtools/getdata/getOutsideData.py
index 6d54b57..799dbab 100755
--- a/murgtools/getdata/getOutsideData.py
+++ b/murgtools/getdata/getOutsideData.py
@@ -396,6 +396,8 @@ def getSatelliteImagery(corners, filename=None, collection='sentinel-2-l2a',
     tiepoint = None
     scale = None
     epsg = None
+    origin_x = None  # Model-space origin for pixel (0,0)
+    origin_y = None
 
     try:
         image = tifffile.imread(tmp_path)
@@ -409,11 +411,16 @@ def getSatelliteImagery(corners, filename=None, collection='sentinel-2-l2a',
                 scale = tags[33550].value     # (scaleX, scaleY, scaleZ)
                 img_h, img_w = tif.pages[0].shape[:2]
 
+                # Compute model-space origin for pixel (0,0)
+                # Tiepoint specifies model coords (x,y) for raster pixel (i,j)
+                origin_x = tiepoint[3] - tiepoint[0] * scale[0]
+                origin_y = tiepoint[4] + tiepoint[1] * scale[1]
+
                 # UTM bounds - normalize with min/max to handle negative scale
-                x_edge1 = tiepoint[3]
-                x_edge2 = tiepoint[3] + img_w * scale[0]
-                y_edge1 = tiepoint[4]
-                y_edge2 = tiepoint[4] - img_h * scale[1]
+                x_edge1 = origin_x
+                x_edge2 = origin_x + img_w * scale[0]
+                y_edge1 = origin_y
+                y_edge2 = origin_y - img_h * scale[1]
                 utm_left = min(x_edge1, x_edge2)
                 utm_right = max(x_edge1, x_edge2)
                 utm_bottom = min(y_edge1, y_edge2)
@@ -489,18 +496,17 @@ def getSatelliteImagery(corners, filename=None, collection='sentinel-2-l2a',
         bbox_utm_south = min(c[1] for c in corners_utm)
         bbox_utm_north = max(c[1] for c in corners_utm)
 
-        # Use GeoTIFF origin and scale for pixel calculation
-        utm_origin_x = tiepoint[3]
-        utm_origin_y = tiepoint[4]
+        # Use GeoTIFF origin (for pixel 0,0) and scale for pixel calculation
         scale_x = scale[0]
         scale_y = scale[1]
 
         # Compute pixel indices using floor/ceil to include full bbox
         # Handle negative scale by computing both edges and normalizing
-        px_west = (bbox_utm_west - utm_origin_x) / scale_x
-        px_east = (bbox_utm_east - utm_origin_x) / scale_x
-        px_north = (utm_origin_y - bbox_utm_north) / scale_y
-        px_south = (utm_origin_y - bbox_utm_south) / scale_y
+        # origin_x/origin_y are model coords for pixel (0,0), computed from tiepoint
+        px_west = (bbox_utm_west - origin_x) / scale_x
+        px_east = (bbox_utm_east - origin_x) / scale_x
+        px_north = (origin_y - bbox_utm_north) / scale_y
+        px_south = (origin_y - bbox_utm_south) / scale_y
 
         # Normalize pixel indices (handle negative scale)
         x1 = math.floor(min(px_west, px_east))
@@ -516,12 +522,11 @@ def getSatelliteImagery(corners, filename=None, collection='sentinel-2-l2a',
             return None
 
         # Compute actual UTM bounds of cropped region
-        # Use the normalized image coordinates (utm_left/right/bottom/top from earlier)
-        # to ensure consistency regardless of scale sign
-        actual_utm_x1 = utm_origin_x + x1 * scale_x
-        actual_utm_x2 = utm_origin_x + x2 * scale_x
-        actual_utm_y1 = utm_origin_y - y1 * scale_y
-        actual_utm_y2 = utm_origin_y - y2 * scale_y
+        # Use origin_x/origin_y (model coords for pixel 0,0) for consistency
+        actual_utm_x1 = origin_x + x1 * scale_x
+        actual_utm_x2 = origin_x + x2 * scale_x
+        actual_utm_y1 = origin_y - y1 * scale_y
+        actual_utm_y2 = origin_y - y2 * scale_y
         actual_utm_west = min(actual_utm_x1, actual_utm_x2)
         actual_utm_east = max(actual_utm_x1, actual_utm_x2)
         actual_utm_south = min(actual_utm_y1, actual_utm_y2)
diff --git a/tests/test_getDataFRF.py b/tests/test_getDataFRF.py
index a6cf5f9..5a9ce56 100644
--- a/tests/test_getDataFRF.py
+++ b/tests/test_getDataFRF.py
@@ -407,6 +407,36 @@ def test_handles_negative_scale_x(self):
         finally:
             os.unlink(tmp_path)
 
+    def test_handles_nonzero_tiepoint_pixel(self):
+        """Test GeoTIFF extent when tiepoint references non-(0,0) pixel."""
+        image = np.zeros((200, 100), dtype=np.uint8)
+        # Tiepoint at pixel (10, 20) -> model coords (900510, 275980)
+        # With scale (1, 1), pixel (0,0) should be at:
+        #   origin_x = 900510 - 10 * 1 = 900500
+        #   origin_y = 275980 + 20 * 1 = 276000
+        # So extent should be [900500, 900600, 275800, 276000]
+        tiepoint = [10.0, 20.0, 0.0, 900510.0, 275980.0, 0.0]
+        scale = [1.0, 1.0, 0.0]
+
+        with tempfile.NamedTemporaryFile(suffix='.tif', delete=False) as tmp:
+            tmp_path = tmp.name
+
+        try:
+            tifffile.imwrite(
+                tmp_path,
+                image,
+                extratags=[
+                    (33922, 'd', 6, tiepoint, True),
+                    (33550, 'd', 3, scale, True),
+                ]
+            )
+
+            extent = get_geotiff_extent(tmp_path)
+            # Verify extent is computed from pixel (0,0), not the tiepoint pixel
+            assert extent == pytest.approx([900500.0, 900600.0, 275800.0, 276000.0])
+        finally:
+            os.unlink(tmp_path)
+
 
 class TestThreadGetArgusImagery:
     """Tests for the threadGetArgusImagery function."""