Fix Argus overlay extent conversion in wave/imagery example by Copilot · Pull Request #57 · SBFRF/murgtools

Copilot · 2026-07-13T13:23:28Z

The Argus overlay in examples/test_wave_and_imagery.py could drift southeast of the satellite basemap because the overlay path assumed the GeoTIFF bounds were always NC State Plane. This change makes the overlay derive lon/lat bounds from the GeoTIFF’s embedded CRS instead of hard-coding a projection assumption.

CRS-aware Argus extent conversion
- Extend get_geotiff_extent() with an optional to_latlon mode.
- Read the GeoTIFF CRS from ProjectedCSTypeGeoKey / GeographicTypeGeoKey.
- Convert projected bounds to [lon_min, lon_max, lat_min, lat_max] for plotting when requested.
Example overlay path
- Update examples/test_wave_and_imagery.py to compute the Argus overlay extent from GeoTIFF metadata:
  - keep the native extent available for debugging/logging
  - use CRS-derived lon/lat bounds for imshow(...)
Regression coverage
- Add focused tests for:
  - native GeoTIFF extent extraction
  - projected GeoTIFF extent conversion to lon/lat

native_extent = get_geotiff_extent(tmp_path)
argus_extent = get_geotiff_extent(tmp_path, to_latlon=True)

ax2.imshow(argus_result["image"], extent=argus_extent, aspect="auto", alpha=0.5)

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>

Copilot

Pull request overview

This PR makes GeoTIFF extent handling CRS-aware so the Argus overlay in the wave/imagery example uses lon/lat bounds derived from the GeoTIFF’s embedded CRS instead of assuming NC State Plane, and adds regression tests for the new behavior.

Changes:

Extend get_geotiff_extent() to optionally convert native projected bounds to [lon_min, lon_max, lat_min, lat_max] using GeoTIFF CRS metadata.
Update examples/test_wave_and_imagery.py to use the CRS-derived lon/lat extent for the Argus overlay.
Add unit tests covering native extent extraction and projected->lat/lon conversion.

Reviewed changes

Copilot reviewed 4 out of 4 changed files in this pull request and generated 5 comments.

File	Description
`tests/test_getDataFRF.py`	Adds tests for native extent extraction and CRS-based lon/lat conversion.
`murgtools/getdata/getOutsideData.py`	Updates satellite imagery cropping to compute crop indices in projected (UTM) coordinates and derive a more accurate output extent.
`murgtools/getdata/getDataFRF.py`	Adds `to_latlon` mode to `get_geotiff_extent()` and implements CRS detection + coordinate conversion.
`examples/test_wave_and_imagery.py`	Switches Argus overlay extent calculation to the CRS-aware GeoTIFF conversion.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

    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


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>

SBFRF · 2026-07-14T13:33:28Z

Review Feedback Addressed in `cf68f3f`

Thanks for the thorough review! All 5 issues have been addressed:

1. Fallback path undefined variables

Added explicit initialization of tiepoint, scale, epsg to None before the try block, with a has_geotiff_metadata flag that tracks successful extraction. The cropping code now checks this flag and falls back to linear lat/lon mapping when GeoTIFF metadata is unavailable.

2. BBox UTM conversion with mixed coordinates

Now converts all 4 bbox corners to UTM and takes min/max eastings/northings, properly handling rotated/skewed AOIs.

3. int() truncation shrinking crop window

Changed to use math.floor() for x1/y1 and math.ceil() for x2/y2 to ensure the full requested bbox is included before clamping to image bounds.

4. resolution_m incorrect for non-UTM CRS

Added fallback path that computes resolution from deg_per_px when GeoTIFF metadata is unavailable, using the meters_per_deg conversion.

5. GeoKey parsing with integer IDs or string values

Updated get_geotiff_extent to support both string key names ('ProjectedCSTypeGeoKey') and integer key IDs (3072/2048), plus handling of string EPSG values like "EPSG:32618".

All 160 tests pass and the visual output remains correctly aligned.

Copilot

Pull request overview

Copilot reviewed 4 out of 4 changed files in this pull request and generated 3 comments.

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>

SBFRF · 2026-07-14T13:52:52Z

Additional Review Feedback Addressed in `c82739e`

6. Negative scale_y flips vertical extent (getDataFRF.py:3388)

Fixed by computing both edges and normalizing with min/max for both X and Y axes:

x_edge1 = tiepoint[3]
x_edge2 = tiepoint[3] + width * scale[0]
left = min(x_edge1, x_edge2)
right = max(x_edge1, x_edge2)
# Same for Y

7. Test coverage for negative scale_y (tests/test_getDataFRF.py:310)

Added two new tests:

test_handles_negative_scale_y: Verifies extent with negative scale_y (row 0 at bottom)
test_handles_negative_scale_x: Verifies extent with negative scale_x (mirrored horizontally)

8. Empty slice causing division by zero (getOutsideData.py:523)

Added check after clamping crop indices in both UTM and fallback paths:

if x2 <= x1 or y2 <= y1:
    return None

All 162 tests pass (including 2 new tests).

Copilot

Pull request overview

Copilot reviewed 4 out of 4 changed files in this pull request and generated 3 comments.

+        # 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)
+
+        # 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
+        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


+    # Track whether we have GeoTIFF metadata for accurate UTM-based cropping
+    has_geotiff_metadata = False
+    tiepoint = None
+    scale = None
+    epsg = None


                        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)


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>

SBFRF · 2026-07-14T14:07:09Z

Review Feedback Addressed in `ef992b2`

9. Negative scale values invert pixel indices (getOutsideData.py:504)

Fixed by computing both pixel index edges and normalizing with min/max:

px_west = (bbox_utm_west - utm_origin_x) / scale_x
px_east = (bbox_utm_east - utm_origin_x) / scale_x
x1 = math.floor(min(px_west, px_east))
x2 = math.ceil(max(px_west, px_east))

Also normalize actual UTM bounds and use abs(scale_x) for resolution_m.

10. PR description should mention satellite behavior change (getOutsideData.py:398)

Updated PR description to explicitly call out the satellite imagery cropping behavior change, including before/after explanation.

11. EPSG detection defaults to 32618 before checking GeoKeys (getOutsideData.py:434)

Fixed by checking GeoKeys (ProjectedCSTypeGeoKey/3072) first, which are the authoritative source, before falling back to parsing the projection string. Only default to EPSG:32618 as absolute last resort.

All 162 tests pass.

Copilot

Pull request overview

Copilot reviewed 4 out of 4 changed files in this pull request and generated 1 comment.

Comments suppressed due to low confidence (1)

murgtools/getdata/getDataFRF.py:3350

The docstring says this computes “geographic bounds”, but the function now explicitly returns native GeoTIFF coordinates unless to_latlon=True. That wording is misleading for projected GeoTIFFs (and for any caller relying on the docstring).

    Parses GeoTIFF tags (ModelTiepointTag and ModelPixelScaleTag) to compute
    geographic bounds suitable for matplotlib imshow extent parameter.

+        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


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>

SBFRF · 2026-07-14T14:45:28Z

Review Feedback Addressed in `e192a2a`

12. resolution_m uses CRS units instead of meters (getOutsideData.py:535)

Fixed by checking the CRS unit using pyproj and converting appropriately:

crs = CRS.from_epsg(epsg)
crs_unit = crs.axis_info[0].unit_name
if crs_unit in ('metre', 'meter', 'm'):
    resolution_m = abs(scale_x)
elif crs_unit in ('degree', 'degrees'):
    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'):
    resolution_m = abs(scale_x) * 0.3048

This ensures resolution_m is always in meters regardless of the GeoTIFF CRS.

All 162 tests pass.

Copilot

Pull request overview

Copilot reviewed 4 out of 4 changed files in this pull request and generated 5 comments.

+        utm_origin_x = tiepoint[3]
+        utm_origin_y = tiepoint[4]
+        scale_x = scale[0]
+        scale_y = scale[1]


+    # 7. Crop to bbox - use UTM coordinates if GeoTIFF metadata available
    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])
+    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)
+


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>

SBFRF · 2026-07-14T15:43:40Z

Review Feedback Addressed in `8284fdf`

13. Tiepoint assumes pixel (0,0) in get_geotiff_extent (getDataFRF.py:3389)

Fixed 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]

14. UTM bounds assumes pixel (0,0) in getSatelliteImagery (getOutsideData.py:416)

Same fix - compute origin for pixel (0,0) before deriving UTM edges.

15. Pixel index calculation assumes pixel (0,0) (getOutsideData.py:496)

Fixed by using origin_x/origin_y (stored from metadata extraction) instead of tiepoint[3]/tiepoint[4].

16. Add test for non-zero tiepoint pixel (tests/test_getDataFRF.py:307)

Added test_handles_nonzero_tiepoint_pixel which uses tiepoint at pixel (10, 20) and verifies the extent is correctly computed from pixel (0,0).

17. PR description should call out getSatelliteImagery changes (getOutsideData.py:478)

PR description was already updated in a previous commit to explicitly document the satellite imagery behavior changes.

All 163 tests pass (1 new test added).

SBFRF · 2026-07-14T15:48:59Z

Final Summary - Ready to Merge

This PR fixes the coordinate alignment between satellite imagery, Argus orthophotos, and wave gauge locations.

Root Cause

The satellite imagery was offset ~784m east and ~182m north due to using a linear lat/lon-to-pixel mapping instead of proper UTM coordinate transformation.

Key Fixes

Satellite Imagery (getSatelliteImagery)

Convert bbox corners to UTM using GeoTIFF's native CRS
Use tiepoint/scale for accurate pixel index calculation
Handle non-zero tiepoint pixel locations (i,j ≠ 0,0)
Handle negative scale values with proper normalization
Check GeoKeys for EPSG before falling back to projection string
Convert resolution_m correctly for non-meter CRS units
Return None for empty crop regions (prevents div/zero)

Argus Imagery (get_geotiff_extent)

Add to_latlon parameter for automatic CRS conversion
Handle non-zero tiepoint pixel locations
Handle negative scale values
Support both string and integer GeoKey IDs

Test Results

163 tests pass (5 new tests added)
Visual verification confirms proper alignment of satellite, Argus, and gauge locations

Review Feedback

All 17 review comments have been addressed across 6 commits.

Initial plan

d3769f4

Copilot AI assigned Copilot and SBFRF Jul 13, 2026

Copilot started work on behalf of SBFRF July 13, 2026 13:23 View session

Copilot AI linked an issue Jul 13, 2026 that may be closed by this pull request

satellite and argus overlap #56

Closed

Copilot AI added 6 commits July 13, 2026 13:27

Fix Argus overlay GeoTIFF extent conversion

57ca429

Address review feedback on GeoTIFF fix

4c0eaac

Polish GeoTIFF extent validation feedback

d1c9fa2

Finalize GeoTIFF extent helper cleanup

ff436b7

Clarify GeoTIFF CRS helper logic

0614dc8

Tidy GeoTIFF extent helper readability

b16e13e

Copilot AI changed the title ~~[WIP] Fix offset issue between satellite and Argus imagery~~ Fix Argus overlay extent conversion in wave/imagery example Jul 13, 2026

Copilot finished work on behalf of SBFRF July 13, 2026 13:33

Copilot AI requested a review from SBFRF July 13, 2026 13:33

SBFRF marked this pull request as ready for review July 14, 2026 13:00

Copilot AI review requested due to automatic review settings July 14, 2026 13:00

Copilot started reviewing on behalf of SBFRF July 14, 2026 13:00 View session