Drop llh fill values before geocoding in geolocate_uavsar

[ZachHoppinen]

Addresses #34

What is going on

geolocate_uavsar writes the raw .llh lat and long grids to tifs and hands them to gdal.Warp as GEOLOCATION arrays. Outside the swath those grids are 0, so unless GDAL honors the nodata flag on the geolocation arrays, the suggested output extent runs from the scene all the way to lat/long 0.

Synthetic 400 x 400 scene near the Jemez site, fill left in, geolocation arrays with no nodata honored:

case                      secs       output shape   valid px
A zeros,   no bounds     12.20         (2, 70650)          0
B zeros,   bounds         0.04          (95, 200)          0
C nan,     no bounds      0.03             (1, 1)          0
D nan,     bounds         0.03          (95, 200)      19000
E nodata0, no bounds      0.03          (95, 200)      19000

Case A is the reported symptom: Too many points (529 out of 529) failed to transform, unable to compute output bounds, a 70,650 px wide grid, and no usable data. That run was at 0.001 deg; geolocate_uavsar uses 0.00005556 deg, 18x finer per axis, which is where 20 minutes comes from. Note case B, explicit bounds on their own are not enough, the fill has to come out of the arrays.

Caveat on verification

On current GDAL (tested 3.11.4 and 3.12.3) the existing code already runs fine, in 0.03 s on my synthetic scene, because rasterio writes nodata = 0 on the lat/long tifs, that survives into the VRT, and modern GDAL skips those pixels when computing the extent (case E). Handling of nodata in geolocation arrays has changed over GDAL's history and was being discussed as a problem in 2019. The issue was filed against GDAL 3.0.2, and I have no GDAL older than 3.11 available to reproduce the original 20 minute run directly, so I cannot claim this closes the issue in the reporter's environment. What it does do is make the result independent of that behavior.

Change

• Mask lat/long fill to NaN and set the tif nodata accordingly.
• Compute output bounds from the valid latitudes and longitudes and pass them to gdal.Warp, rather than letting GDAL infer the extent (case D above).
• yRes used spacing[0], so a non square spacing was silently ignored.
• The band loop read GetRasterBand(1) every iteration, so a multiband input took band 1's datatype for every band.
• Removed the printed "Ignore the error message: Unable to compute bounds" note, since that error is the failure being fixed.

Check

End to end run of geolocate_uavsar on a synthetic 300 x 300 scene whose first 30 rows are .llh fill, gdal 3.12.3:

runtime      : 0.03 s
output shape : 324 x 360 px
output extent: lon -106.5200 to -106.5000   lat 35.8800 to 35.8980
source extent: lon -106.5200 to -106.5000   lat 35.8800 to 35.8980
valid pixels : 116640 of 116640
value range  : 31.0 to 597.0  (source 1.0 to 599.0)

The extent matches the valid part of the scene and the low value is 31, the first row inside the swath, so the fill rows are dropped and nothing else is. Same output as pre-patch on this GDAL, one column narrower because the bounds now snap to the valid lat/long rather than to a nodata derived extent.

@jacktarricone if you still have the Jemez files, running this branch on them would say whether the 20 minute case is really gone on an older GDAL.

🤖 Generated with Claude Code

https://claude.ai/code/session_01AGcA5iStDRm4rsucaH5iDz