reproject: fix half-pixel offset in geoid and datum-shift grid lookups (#2508)#2516
Merged
brendancol merged 2 commits intoMay 28, 2026
Conversation
xarray-contrib#2508) `_interp_geoid_point` and `_grid_interp_point` indexed pixel-edge coordinates into pixel-center anchored grids. The bilinear stencil landed half a pixel off, producing up to 2.1 m error in the EGM96 geoid undulation (mean 24 cm, p95 69 cm) and an analogous offset on NADCON / OSTN / NTv2 horizontal datum shift grids. Subtract 0.5 from the row/col index so `data[r, c]` is queried at its pixel-center coordinate, matching how the main `reproject` pixel mapping already handles source pixel coords. Cross-checked against pyproj: New York geoid agrees to ~1e-10 m post-fix vs the 9 cm error before. Add `TestGeoidPixelCenterIndexing` with three regressions: a pixel-center identity test on the real EGM96 grid, a sub-cm pyproj cross-check at five locations including the antimeridian pole where the offset error was largest, and a synthetic-grid test for `_grid_interp_point` so the datum-grid fix doesn't depend on a downloaded NADCON file.
CI runners without local PROJ grid downloads and with network access disabled get a silent no-op transform from pyproj, returning 0 m at all locations instead of looking the geoid up. The cross-check test then sees our correct -32.77 m at New York vs pyproj's 0 m and fails. Probe at New York (~-32.8 m) before iterating sample points; skip if the magnitude is < 1 m (the no-grid fallback signature).
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
_interp_geoid_pointand_grid_interp_pointindexed pixel-edge coordinates into pixel-center anchored grids. The bilinear stencil landed half a pixel off, producing up to 2.1 m error on the EGM96 geoid (mean 24 cm, p95 69 cm at pixel centers) and an analogous offset on NADCON / OSTN / NTv2 horizontal datum shift grids. Cross-checked against pyproj: ~9 cm New York geoid error before the fix, ~1e-10 m after.Closes #2508.
What changed
_interp_geoid_point(inxrspatial/reproject/_vertical.py): subtract 0.5 from the row/col index sodata[r, c]is queried at its pixel-center coordinate. Edge handling extended so queries inside the outer pixel center's half-pixel border return the edge pixel value rather than NaN, matching PROJ'svgridshift. Negativecol_fis handled viamath.floorso the longitude wrap still picks the antipodal column correctly._grid_interp_point(inxrspatial/reproject/_datum_grids.py): same-0.5correction. Out-of-grid queries still return(0, 0)matching PROJ'shgridshift.TestGeoidPixelCenterIndexinginxrspatial/tests/test_reproject.py: three regressions_grid_interp_pointso it doesn't depend on a downloaded NADCON fileWhy existing tests missed it
TestGeoidHeightBehaviourchecked reference values with a 3 m tolerance, which is wider than the 2.1 m worst-case offset error.Test plan
pytest xrspatial/tests/test_reproject.py -k 'geoid or vertical or datum'-- 47 passed (44 prior + 3 new)pytest xrspatial/tests/test_reproject.py-- 348 passed (was 345)Found by
/sweep-accuracyagainstreprojecton 2026-05-27.