Inherit georef from level-0 IFD when reading overview IFDs (#1640)#1641
Merged
brendancol merged 2 commits intoMay 12, 2026
Conversation
GDAL-style COG writers, including this package's own to_geotiff, put the GeoKeyDirectory, ModelPixelScale, and ModelTiepoint only on the level-0 IFD. Before this fix, open_geotiff(path, overview_level=N) with N >= 1 returned a DataArray whose CRS, transform attr, and y/x coords were silently replaced with defaults (unit transform, integer pixel indices, no crs). Downstream spatial ops then computed pixel sizes from the wrong coords and produced silently-wrong answers. Add extract_geo_info_with_overview_inheritance in _geotags.py. When the selected IFD is a reduced-resolution overview (NewSubfileType bit 0) and its own extract_geo_info reports has_georef=False, the helper re-runs extract_geo_info on the first full-resolution IFD and rescales the pixel size by width_full / width_overview (and the same for height) before copying CRS-side metadata across. Overviews that already carry their own valid georef (some writers replicate it) are left untouched. Wire the helper into every read entry point: - read_to_array (eager numpy and dask chunk workers) - _read_cog_http (HTTP COGs) - _read_geo_info (dask metadata probe) - read_geotiff_gpu (cupy and dask+cupy paths) Add xrspatial/geotiff/tests/test_overview_geo_inheritance_1640.py covering all four backends (numpy, dask+numpy, cupy, dask+cupy) plus two edge cases: an overview that carries its own geokeys keeps them, and a file with no full-res sibling falls back without raising.
This was referenced May 12, 2026
brendancol
added a commit
that referenced
this pull request
May 12, 2026
PR #1641 (issue #1640) inherits level-0 georef on overview reads but keeps the level-0 origin_x / origin_y unchanged. That is correct for PixelIsArea -- origin is the upper-left corner of pixel (0, 0), which is also the upper-left corner of the overview's pixel (0, 0). It is wrong for PixelIsPoint (GeoKey 1025 = 2): the origin is the center of pixel (0, 0), and an overview pixel covering the first scale_x columns of level 0 has its center at the centroid of those level-0 pixels (origin + (scale - 1) * 0.5 * pixel_size_lvl0). Before this fix, open_geotiff on a 1024x1024 PixelIsPoint COG with 10 m pixels and origin (0, 0) returned x[:3] = [0, 20, 40] for overview_level=1 instead of [5, 25, 45]. Downstream sel / interp / reproject silently snaps to the wrong pixel for any DEM-style PixelIsPoint COG (USGS, OpenTopography, Copernicus DEM). Fix in extract_geo_info_with_overview_inheritance: choose the effective raster_type first (overview's own when it explicitly declared non-default, otherwise inherit from level 0), then apply origin_shift = (scale - 1) * 0.5 * pixel_size_lvl0 along each axis when that effective raster_type is PixelIsPoint. The PixelIsArea path is byte-equivalent to before. Add 13 regression tests in test_overview_pixel_is_point_1642.py: centroid identity across all four backends, transform-tuple values across all four backends, uniform grid step, unit-level helper tests for both raster_types via stubbed extract_geo_info, an own-geokeys- not-clobbered case on PixelIsPoint, and a PixelIsArea regression check so the #1640 contract still holds. All 1397 existing non-network geotiff tests still pass.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Closes #1640.
open_geotiff(path, overview_level=N)withN >= 1was silentlydropping CRS, the
transformattr, and georeferenced y/x coords.GDAL-style COG writers (including this package's
to_geotiff) putGeoKeys / ModelPixelScale / ModelTiepoint only on the level-0 IFD, so
the reader's per-IFD
extract_geo_infocall returned an emptyGeoTransformand no CRS for every overview read. Downstream spatialops then derived pixel sizes from the integer fallback coords and got
silently-wrong answers.
The fix lives in a new helper,
extract_geo_info_with_overview_inheritance. When the selected IFDis a reduced-resolution overview (NewSubfileType bit 0) and its own
georef is missing, the helper re-runs
extract_geo_infoon the firstfull-resolution IFD and rescales the pixel size by
width_full / width_overview(same for height) before copying theCRS-side metadata across. Overviews that already carry their own
geokeys (some writers replicate them) are returned unchanged.
The helper is wired into every read entry point so all four backends
benefit:
read_to_array(eager numpy and dask chunk workers),_read_cog_http(HTTP COGs),_read_geo_info(dask metadata probe),and
read_geotiff_gpu(cupy and dask+cupy).Test plan
tests/test_overview_geo_inheritance_1640.py:of the four backends (numpy, dask+numpy, cupy, dask+cupy)
deepcopy failures in
test_features.pyunrelated to this change)xrspatial/tests/suite passes (3528 passed, 15 skipped)