reproject: geoid undulation lookup has half-pixel offset error (up to 2 m)

[brendancol]

Summary

xrspatial.reproject._vertical._interp_geoid_point (and the parallel _datum_grids._grid_interp_point) compute the bilinear interpolation index as if the geoid grid stored values at pixel-edge anchors, while the grid is actually pixel-center anchored. This produces a systematic half-pixel offset on every lookup.

For the vendored EGM96 15-arcmin grid this corresponds to a roughly 14 km horizontal shift in the geoid sampling location. The resulting error in the geoid undulation N depends on the local gradient; sampled at every 50th pixel center it is:

mean abs error : 0.24 m
p95  abs error : 0.69 m
max  abs error : 2.13 m

The same offset is present in _grid_interp_point, which is used to apply NADCON / OSTN / NTv2 horizontal datum shift grids; for those grids the offset is hundreds of metres of horizontal misregistration (the actual error in degrees depends on grid resolution and local gradient).

Reproducer

from xrspatial.reproject._vertical import _interp_geoid_point, _load_geoid

data, left, top, res_x, res_y, h, w = _load_geoid("EGM96")
# Pixel (0,0) center: lon=-180, lat=90.
# data[0,0] = 13.606 (the geoid undulation actually stored at that point)
center_lon = left + 0.5 * res_x   # -180.0
center_lat = top  - 0.5 * res_y   #  90.0

N = _interp_geoid_point(center_lon, center_lat, data, left, top,
                        res_x, res_y, h, w)
print(N, "vs stored", data[0, 0])
# 13.5497 vs stored 13.6062 -- a half-pixel blend with the next column

Cross-checked against pyproj's Transformer.from_crs('EPSG:4979', 'EPSG:5773'):

N at (-180, 90)  : pyproj 13.6062,  ours 13.5497  (diff 5.6 cm)
N at ( -74, 40.7): pyproj -32.7736, ours -32.8596 (diff 8.6 cm)

Root cause

_interp_geoid_point computes

col_f = (lon - left) / res_x          # position in pixel-edge units
row_f = (top - lat) / res_y

then does bilinear interpolation between data[r0, c0] and the next pixel. But the EGM96 GeoTIFF is pixel-center anchored: data[0, 0] is the value at (left + 0.5*res_x, top - 0.5*res_y), not at (left, top). To index into a pixel-center grid the correct continuous coord is (lon - left)/res_x - 0.5.

The same fix applies to _grid_interp_point in _datum_grids.py.

For comparison, the main reproject path already subtracts 0.5 when converting source CRS coords to source pixel indices (see _reproject_chunk_numpy: src_col_px = (src_x - src_left) / src_res_x - 0.5). So this is just a local inconsistency in the geoid / datum-grid path, not a deeper design issue.

Why existing tests miss it

tests/test_reproject.py::TestGeoidHeight::test_geoid_height_scalar checks reference values with tolerance < 3.0 m, which is looser than the 2.1 m worst-case offset error. The _REFERENCE_N dict was set with that loose tolerance precisely because the lookup itself was the source of the small disagreement against published values.

Categories

• Cat 4: Missing or wrong Earth curvature / projection corrections (half-pixel indexing inconsistent with how reproject itself handles pixel-center coords).

Severity

HIGH. The error is systematic, affects every geoid-shifted reprojection and every NADCON / OSTN / NTv2 datum shift, and falls in the 0.1-2 m band that matters for surveying / bathymetry / chart datum work which the vertical CRS API specifically targets.

Found by /sweep-accuracy against the reproject module on 2026-05-27.