reproject: fix half-pixel offset in geoid and datum-shift grid lookups (#2508)

.claude/sweep-accuracy-state.csv

@@ -25,7 +25,7 @@ perlin,2026-04-10T12:00:00Z,,,,Improved Perlin noise implementation correct. Fad
 polygon_clip,2026-04-13T12:00:00Z,1197,,,crop=True + all_touched=True drops boundary pixels. Fix in PR #1200.
 polygonize,2026-05-19,2151,HIGH,5,cupy backend used exact-equality vs numba _is_close tolerance for float dtypes; fixed via tolerance-aware value grouping
 proximity,2026-03-30T15:00:00Z,,,,Direction >= boundary fragile but works due to truncated constant. Float32 truncation is design choice. No wrong-results bugs found.
-reproject,2026-05-17,2025,HIGH,5,"HIGH: _apply_vertical_shift crashed on dask (no fancy-index setitem), cupy (Numba JIT can't take cupy), and 3D multi-band (broadcast mismatch). Fixed by routing dask through map_blocks, round-tripping cupy via host, and looping per band. Previously documented MEDIUM (GPU dispatcher missing NAD27/Helmert wrapper) and MEDIUM (legacy _resample_cupy cubic NaN threshold 1e-6 diverges from CPU bilinear-fallback) still stand. Identity-CRS cupy reproject returns zeros - separate pre-existing bug, not in scope."
+reproject,2026-05-27,2508,HIGH,4;5,"HIGH (#2508, fixed): _interp_geoid_point and _grid_interp_point indexed pixel-edge coords on pixel-center anchored grids; up to 2.1 m geoid error (mean 24 cm, p95 69 cm), ~9 cm cross-check vs pyproj at NYC; same offset in NADCON/OSTN/NTv2 horizontal datum grids. Fixed by subtracting 0.5 from row/col indices; added TestGeoidPixelCenterIndexing with pyproj cross-check. MEDIUM (deferred): _resample_cupy cubic uses cardinal spline (cupyx.scipy.ndimage.map_coordinates order=3) while numpy / _resample_cupy_native use Catmull-Rom; ~0.007 absolute divergence on a smooth quadratic. Fires when CRS pair has no CUDA fast path. MEDIUM (deferred): _reproject_dask(is_cupy=True) returns dask with cupy chunks but numpy meta; _apply_vertical_shift_dask would call _apply_vertical_shift_numpy on cupy blocks and crash on np.asarray(cupy). Untriggered by tests (needs >GPU output). LOW: try_cuda_transform doesn't apply Helmert datum shifts (only WGS84/NAD83 supported) so non-WGS84 datums fall back to pyproj instead of running the GPU fast path; correctness is OK but performance asymmetric vs CPU path."
 resample,2026-04-14T12:00:00Z,1202,,,Interpolation used edge-aligned coords instead of block-centered. Fix in PR #1204.
 sieve,2026-04-13T12:00:00Z,,,,Union-find CCL correct. NaN excluded from labeling. All backends funnel through _sieve_numpy.
 sky_view_factor,2026-05-01,1407,HIGH,4,Horizon angle ignored cell size; fixed by passing cellsize_x/cellsize_y into CPU+GPU kernels and using ground distance

xrspatial/reproject/_datum_grids.py

@@ -226,11 +226,23 @@ def _grid_interp_point(lon, lat, dlat_grid, dlon_grid,
                        grid_h, grid_w):
     """Bilinear interpolation of a single point in the shift grid.
 
+    The shift grids are pixel-center anchored (matches rasterio's
+    ``bounds`` convention): ``dlat_grid[r, c]`` is the offset at
+    ``(grid_left + (c + 0.5) * grid_res_x, grid_top - (r + 0.5) * grid_res_y)``.
+    Indexing into the array therefore needs ``(coord - edge) / res - 0.5``,
+    not ``(coord - edge) / res``; without the ``-0.5`` shift the lookup
+    is biased by half a pixel, which produces a sub-decimetre to
+    decimetre horizontal residual depending on grid resolution
+    (NADCON5: ~5 mm typical, NTv2: 10+ cm). See GH #2508.
+
     Returns (dlat_arcsec, dlon_arcsec) or (0, 0) if outside the grid.
     """
-    col_f = (lon - grid_left) / grid_res_x
-    row_f = (grid_top - lat) / grid_res_y
+    col_f = (lon - grid_left) / grid_res_x - 0.5
+    row_f = (grid_top - lat) / grid_res_y - 0.5
 
+    # Outside the rectangle bounded by the outer pixel centers: no shift
+    # rather than an extrapolated one. This matches PROJ's hgridshift,
+    # which only applies inside the populated grid extent.
     if col_f < 0 or col_f > grid_w - 1 or row_f < 0 or row_f > grid_h - 1:
         return 0.0, 0.0
 

xrspatial/reproject/_vertical.py

@@ -121,29 +121,54 @@ def _load_geoid(model='EGM96'):
 
 @njit(nogil=True, cache=True)
 def _interp_geoid_point(lon, lat, data, left, top, res_x, res_y, h, w):
-    """Bilinear interpolation of geoid undulation at a single point."""
+    """Bilinear interpolation of geoid undulation at a single point.
+
+    The geoid GeoTIFF is pixel-center anchored: ``data[r, c]`` is the
+    value at ``(left + (c + 0.5) * res_x, top - (r + 0.5) * res_y)``.
+    The continuous interpolation index is therefore
+    ``(coord - edge) / res - 0.5`` rather than ``(coord - edge) / res``.
+    Without the ``-0.5`` correction the lookup at a pixel center
+    returned a blend of the target pixel and its neighbour, giving up to
+    ~2 m error on EGM96 at 15 arcmin. See GH #2508.
+    """
     # Wrap longitude to [-180, 180)
     lon_w = lon
     while lon_w < -180.0:
         lon_w += 360.0
     while lon_w >= 180.0:
         lon_w -= 360.0
 
-    col_f = (lon_w - left) / res_x
-    row_f = (top - lat) / res_y
-
-    if row_f < 0 or row_f > h - 1:
-        return math.nan  # outside latitude range
-
-    # Wrap column for global grids
-    c0 = int(col_f) % w
+    # Pixel-center index space: data[r, c] sits at fractional index (r, c).
+    col_f = (lon_w - left) / res_x - 0.5
+    row_f = (top - lat) / res_y - 0.5
+
+    # Latitude beyond the outermost pixel centers is treated as outside
+    # the defined region. Querying exactly on the southern edge (lat ==
+    # bottom + 0.5*res_y) gives row_f == h-1 which is still in-bounds.
+    if row_f < -0.5 or row_f > h - 0.5:
+        return math.nan
+
+    # Clamp half-pixel border to the edge pixel center so a query that
+    # falls between the outer pixel center and the raster bound returns
+    # that pixel's value rather than NaN. This matches PROJ's edge
+    # handling for vgridshift.
+    if row_f < 0.0:
+        row_f = 0.0
+    if row_f > h - 1:
+        row_f = h - 1.0
+
+    # Wrap column for global grids. ``int(...)`` rounds toward zero so
+    # negative fractions need ``math.floor`` to stay on the left
+    # neighbour (col_f = -0.5 -> c0 = -1 -> w-1, the antipodal column).
+    c0_raw = int(math.floor(col_f))
+    c0 = c0_raw % w
     c1 = (c0 + 1) % w
     r0 = int(row_f)
     if r0 >= h - 1:
         r0 = h - 2
     r1 = r0 + 1
 
-    dc = col_f - int(col_f)
+    dc = col_f - c0_raw
     dr = row_f - r0
 
     N = (data[r0, c0] * (1.0 - dr) * (1.0 - dc) +

xrspatial/tests/test_reproject.py

@@ -4163,6 +4163,123 @@ def test_geoid_height_raster_with_lat_lon_dims(self):
         assert np.isfinite(out.values).all()
 
 
+class TestGeoidPixelCenterIndexing:
+    """Regression coverage for the half-pixel offset bug (#2508).
+
+    The EGM96 GeoTIFF is pixel-center anchored: ``data[r, c]`` is the
+    value at ``(left + (c + 0.5) * res_x, top - (r + 0.5) * res_y)``.
+    Before #2508 the bilinear lookup indexed in pixel-edge space, which
+    produced up to ~2 m error at pixel centers and an 8-9 cm error at
+    representative locations like New York vs pyproj's geoid lookup.
+    """
+
+    def test_geoid_at_pixel_center_returns_stored_value(self):
+        """A query at the exact pixel center must return the stored cell
+        value (modulo float round-off), not a blend with the neighbour.
+        """
+        from xrspatial.reproject._vertical import (
+            _interp_geoid_point, _load_geoid,
+        )
+
+        data, left, top, res_x, res_y, h, w = _load_geoid('EGM96')
+
+        for (i, j) in [(0, 0), (10, 100), (h // 2, w // 2),
+                       (h - 1, w - 1)]:
+            lon_c = left + (j + 0.5) * res_x
+            lat_c = top - (i + 0.5) * res_y
+            N = _interp_geoid_point(
+                lon_c, lat_c, data, left, top, res_x, res_y, h, w,
+            )
+            assert abs(N - data[i, j]) < 1e-9, (
+                f"pixel ({i},{j}) center query expected "
+                f"data[{i},{j}]={data[i, j]!r}, got {N!r}; "
+                f"half-pixel offset bug from #2508?"
+            )
+
+    def test_geoid_height_matches_pyproj_within_cm(self):
+        """``geoid_height`` must agree with pyproj's EGM96 lookup to the
+        centimetre at well-sampled locations. The old half-pixel bias was
+        ~9 cm at New York; this test would fail by ~9 cm if reintroduced.
+        """
+        pyproj = pytest.importorskip('pyproj')
+        from xrspatial.reproject import geoid_height
+
+        src_crs = pyproj.CRS('EPSG:4979')
+        tgt_crs = pyproj.CRS('EPSG:5773')
+        transformer = pyproj.Transformer.from_crs(
+            src_crs, tgt_crs, always_xy=True,
+        )
+
+        # pyproj silently falls back to a no-op transform when the EGM96
+        # grid is not installed locally and PROJ network access is
+        # disabled (typical CI). Probe at New York: a real lookup gives
+        # ~-32.8 m, the no-grid fallback gives 0. Skip in the fallback
+        # case -- there's nothing to cross-check against.
+        _, _, h_probe = transformer.transform(-74.0, 40.7, 0.0)
+        if abs(h_probe) < 1.0:
+            pytest.skip(
+                "pyproj EGM96 grid unavailable on this runner "
+                "(transform returned ~0 at New York); cannot cross-check"
+            )
+
+        sample_points = [
+            (-74.0, 40.7),
+            (0.0, 0.0),
+            (139.7, 35.7),
+            (-150.0, 60.0),
+            (-180.0, 90.0),  # data[0,0]: the offset bug was largest here
+        ]
+        for lon, lat in sample_points:
+            _, _, h_ortho = transformer.transform(lon, lat, 0.0)
+            N_expected = -h_ortho  # h_ellip(=0) - h_ortho = -h_ortho
+            N_actual = geoid_height(lon, lat)
+            assert abs(N_actual - N_expected) < 1e-2, (
+                f"N({lon},{lat}) = {N_actual}, pyproj says "
+                f"{N_expected}; diff {N_actual - N_expected:.4f} m"
+            )
+
+    def test_grid_interp_point_pixel_center_returns_stored_value(self):
+        """``_grid_interp_point`` (datum shift grids) has the same
+        pixel-center anchoring as the geoid grid and was fixed in #2508.
+        Verify with a synthetic grid so the test doesn't depend on a
+        downloaded NADCON file.
+        """
+        from xrspatial.reproject._datum_grids import _grid_interp_point
+
+        # 4x5 synthetic grid with distinctive values; pixel-center anchored.
+        dlat_grid = np.array([
+            [10.0, 20.0, 30.0, 40.0, 50.0],
+            [11.0, 22.0, 33.0, 44.0, 55.0],
+            [12.0, 24.0, 36.0, 48.0, 60.0],
+            [13.0, 26.0, 39.0, 52.0, 65.0],
+        ], dtype=np.float64)
+        dlon_grid = dlat_grid * 2.0
+        grid_h, grid_w = dlat_grid.shape
+
+        grid_left = -110.0
+        grid_top = 45.0
+        grid_res_x = 1.0
+        grid_res_y = 1.0
+
+        for i in range(grid_h - 1):
+            for j in range(grid_w - 1):
+                lon_c = grid_left + (j + 0.5) * grid_res_x
+                lat_c = grid_top - (i + 0.5) * grid_res_y
+                dlat, dlon = _grid_interp_point(
+                    lon_c, lat_c, dlat_grid, dlon_grid,
+                    grid_left, grid_top, grid_res_x, grid_res_y,
+                    grid_h, grid_w,
+                )
+                assert abs(dlat - dlat_grid[i, j]) < 1e-12, (
+                    f"pixel ({i},{j}) center: expected dlat "
+                    f"{dlat_grid[i, j]}, got {dlat}"
+                )
+                assert abs(dlon - dlon_grid[i, j]) < 1e-12, (
+                    f"pixel ({i},{j}) center: expected dlon "
+                    f"{dlon_grid[i, j]}, got {dlon}"
+                )
+
+
 class TestVerticalHelperConversions:
     """Direct coverage for the four public vertical-conversion helpers.