Two related problems in xrspatial/reproject around the vertical-datum-shift code path.
Hang on non-finite coords
_apply_vertical_shift in xrspatial/reproject/__init__.py (around lines 702-791) hands output-grid coordinates to _interp_geoid_2d directly. When the target CRS is projected, output coords are first inverse-projected to lon/lat at lines 773-776. Some projections (polar stereographic near the pole, the antimeridian on certain CRSes) emit inf or NaN from transformer.transform. Those values flow into _interp_geoid_point in xrspatial/reproject/_vertical.py (lines 127-130):
while lon_w < -180.0:
lon_w += 360.0
while lon_w >= 180.0:
lon_w -= 360.0
On inf, neither loop terminates and the worker hangs.
The public geoid_height API guards against non-finite input at _vertical.py:205-216, but _apply_vertical_shift calls the JIT batch directly and skips that guard.
Fix
After the optional inverse-projection step, replace any non-finite values in xx_strip and yy_strip with NaN before calling _interp_geoid_2d. The interpolator already returns math.nan for out-of-range input, so NaN-in NaN-out is fine. Update the per-strip is_valid mask to drop pixels with non-finite coordinates so NaN from the geoid does not contaminate the output.
The change is under ten lines, scoped to _apply_vertical_shift. _vertical.py is left alone.
No coverage for the vertical CRS path
The src_vertical_crs and tgt_vertical_crs parameters of reproject() (defined at __init__.py:507-516) have no end-to-end tests in xrspatial/tests/test_reproject.py. _apply_vertical_shift is roughly 90 lines and is currently exercised only through unit tests of geoid_height.
Tests to add
A TestVerticalShift class covering:
- EGM96 to ellipsoidal on a small lon/lat raster, verifying the shift matches the known undulation at a sample location.
- The reverse direction, confirming opposite sign.
- EGM96 to EGM2008, skipped if the EGM2008 grid is not available.
- Same-vertical-CRS no-op.
- Single-side vertical CRS no-op.
- A projected source/target (UTM 33N) to exercise the inverse-projection branch.
- A polar-stereographic case (EPSG:3413) on a 89-90 N raster to regression-test the hang. Assertion is that the call returns and produces finite output where the source is finite.
Two related problems in
xrspatial/reprojectaround the vertical-datum-shift code path.Hang on non-finite coords
_apply_vertical_shiftinxrspatial/reproject/__init__.py(around lines 702-791) hands output-grid coordinates to_interp_geoid_2ddirectly. When the target CRS is projected, output coords are first inverse-projected to lon/lat at lines 773-776. Some projections (polar stereographic near the pole, the antimeridian on certain CRSes) emitinfor NaN fromtransformer.transform. Those values flow into_interp_geoid_pointinxrspatial/reproject/_vertical.py(lines 127-130):On
inf, neither loop terminates and the worker hangs.The public
geoid_heightAPI guards against non-finite input at_vertical.py:205-216, but_apply_vertical_shiftcalls the JIT batch directly and skips that guard.Fix
After the optional inverse-projection step, replace any non-finite values in
xx_stripandyy_stripwith NaN before calling_interp_geoid_2d. The interpolator already returnsmath.nanfor out-of-range input, so NaN-in NaN-out is fine. Update the per-stripis_validmask to drop pixels with non-finite coordinates so NaN from the geoid does not contaminate the output.The change is under ten lines, scoped to
_apply_vertical_shift._vertical.pyis left alone.No coverage for the vertical CRS path
The
src_vertical_crsandtgt_vertical_crsparameters ofreproject()(defined at__init__.py:507-516) have no end-to-end tests inxrspatial/tests/test_reproject.py._apply_vertical_shiftis roughly 90 lines and is currently exercised only through unit tests ofgeoid_height.Tests to add
A
TestVerticalShiftclass covering: