From 66ffcd7ab414848518618e45c003a4226c451671 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Fri, 12 Jun 2026 07:51:50 -0700 Subject: [PATCH 1/3] Record accuracy sweep state for reproject (#3274, #3275, #3276) --- .claude/sweep-accuracy-state.csv | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.claude/sweep-accuracy-state.csv b/.claude/sweep-accuracy-state.csv index 9896a4df0..8c169280c 100644 --- a/.claude/sweep-accuracy-state.csv +++ b/.claude/sweep-accuracy-state.csv @@ -29,7 +29,7 @@ polygon_clip,2026-06-10,3186,HIGH,5,"Cat5 backend inconsistency: dask+cupy clip_ polygonize,2026-05-29,2606,HIGH,5,"Cat 5 HIGH: dask connectivity=8 cross-chunk merge filled diagonal notch where same-value regions meet only at a corner across a chunk boundary; total area exceeded raster. Hole ring was dropped because containment tested hole[0] (on exterior at pinch). Fixed via _ring_interior_point in PR for #2606. numpy, dask+numpy, dask+cupy area parity now holds; 4-conn was already correct. cupy + dask+cupy paths validated on GPU host. Other cats clean: NaN masked on numpy/cupy float paths (tested), _is_close handles +/-inf via exact-equality short-circuit, atol/rtol/simplify_tolerance reject NaN/inf, integer GPU CCL matches numpy." proximity,2026-06-09,3108,HIGH,4;5,"Cat5/Cat4: bounded GREAT_CIRCLE dask (numpy+cupy) missed targets across the +/-180 antimeridian seam: _halo_depth sized x-halo as linear parallel-arc sum, but haversine is periodic in lon and chords shorten near poles, so array-space adjacency is no lower bound on spherical distance; numpy/cupy (brute force) found the wrap target (~111 km), dask returned NaN. Fixed in #3108 via chord bound 2R asin(cos(lat_max)|sin(dlon/2)|) + x-axis fold when seam/180-deg chord within max_distance (covers over-pole too). CUDA host: cupy + dask+cupy executed, 417+ tests pass. Cat1-3 clean (float32 output documented; NaN via isfinite consistent; bounds guards correct; tie-break unified in #2881). LOW (not fixed): great_circle_distance uses WGS84 equatorial radius 6378137 as sphere radius (~0.1% vs mean-radius convention) but documented and exposed as param." rasterize,2026-06-09,3085,MEDIUM,2,"Cat2: non-finite burn values (NaN/inf, e.g. GeoDataFrame column with missing data) against integer output dtype silently cast to platform sentinel (NaN -> INT_MIN) on all 4 backends; bool dtype collapses NaN to True. numpy path suppresses the cast RuntimeWarning so it is fully silent. Fix #3085 mirrors NaN-fill guard #2504 and unsafe-int guard #3056; merge='count' exempt (never reads props). Cats 1/3/4/5 clean: CUDA available, ran 4-backend parity probe (mixed polygon/line/point, all 6 builtin merges, all_touched, chunked 17x23) -- bit-identical across numpy/cupy/dask+numpy/dask+cupy; 665 rasterize tests pass. GPU atomic min/max NaN handling (#2255) verified correct; non-atomic _merge_min_gpu/_merge_max_gpu device fns lack the NaN branch but are unreachable dead code (string merges always take the atomic ladder). Scanline ceil (GPU int-trunc emulation) matches CPU np.ceil for negative x. No Earth-curvature surface (planar vector rasterization by design)." -reproject,2026-06-09,3094,HIGH,4;5,"HIGH #3094: try_cuda_transform missing the #2651 non-WGS84 datum guard; cupy/dask+cupy 4326<->27700 (and DHDN/MGI/ED50/NAD27 entries) project with WGS84 Krueger + no datum shift, ~80-100 m coordinate error, end-to-end numpy-vs-cupy value diff 0.094 on [0,1] data with 5 NaN-mask flips (verified on GPU). MEDIUM #3096: empty-chunk returns in _reproject_chunk_numpy/_cupy and _reproject_block_adapter hardcode float64, so integer-source dask reproject computes float64 when any chunk misses the source footprint while meta advertises int (eager backend returns int). LOW (doc only): _place_same_crs accepts 1% res mismatch for direct placement, up to ~1% of tile width misregistration at the far edge on large same-CRS merge tiles. Cats 1-3 clean: f64 kernels, GDAL-style NaN renorm, correct bounds guards; projection series/constants match PROJ. CUDA available; GPU paths executed." +reproject,2026-06-12,3274,HIGH,1;4,"3 confirmed bugs, all kernel-vs-PROJ parity: #3274 HIGH LAEA inverse spurious /rq (2.6 km err for 3035) + _authalic_apa inverse series wrong (4.8 m in AEA/CEA inverses; PROJ 3-term = 1.6 mm), CPU+CUDA kernels both; #3275 HIGH _is_wgs84_compatible_ellipsoid passes R-defined spheres (MODIS sinusoidal 18.9 km err) and _aea_params/_cea_params lack the guard entirely (23.8 km on spherical aea/cea); #3276 MEDIUM itrf helmert scale 1e-9 but PROJ +s is ppm (1e-6), ~23 mm err. Verified clean: merc/emerc/UTM/tmerc/LCC/polar stere (incl lat_ts akm1) forward+inverse <=1e-5 m vs pyproj; resampling kernels NaN handling and GDAL renorm match across numpy/cupy/dask (CUDA run, gpu-vs-cpu 1.3e-7); dask footprint chunk-skip bbox is a superset in all probed cases (no holes). LOW (documented only): _source_footprint_in_target probe array typo uses x-midpoint mx as a latitude in last 3 ys entries (bbox superset, correctness unaffected)." resample,2026-05-29,2610,HIGH,3;5,"dask interp (nearest/bilinear) overlap depth=1 too small on downsample; block-centered source coord landed past chunk, map_coordinates clamped to edge -> wrong seam rows. Fixed PR #2627 via per-axis _downsample_radius. cupy+dask+cupy verified." 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 From df02e9df25d18387e47cd4497f3d31f19b878b1c Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Fri, 12 Jun 2026 08:00:39 -0700 Subject: [PATCH 2/3] Reject non-WGS84 ellipsoids in the projection fast-path dispatch (#3275) _is_wgs84_compatible_ellipsoid only looked at the ellps/datum names and treated empty values as WGS84. PROJ normalizes spherical definitions (+a=R +b=R) to a bare R key with no ellps, so the MODIS sinusoidal grid passed the check and ran the WGS84 meridional-arc kernel against a sphere: ~18.9 km error. Reject R-defined spheres and explicit semi-axes that do not match WGS84/GRS80 (0.5 m tolerance keeps GRS80 in). _aea_params and _cea_params had no ellipsoid check at all, letting sphere-based AEA/CEA definitions through with up to ~24 km error; add the same guard every other param extractor already calls. Rejected pairs return None from the dispatchers and fall back to pyproj, matching the non-WGS84 datum handling from #2651. The CUDA dispatch shares the helper and param extractors, so the GPU path is covered by the same change. --- xrspatial/reproject/_projections.py | 28 +++ ...t_reproject_sphere_ellipsoid_guard_3275.py | 183 ++++++++++++++++++ 2 files changed, 211 insertions(+) create mode 100644 xrspatial/tests/test_reproject_sphere_ellipsoid_guard_3275.py diff --git a/xrspatial/reproject/_projections.py b/xrspatial/reproject/_projections.py index 137be50eb..3f70f77fd 100644 --- a/xrspatial/reproject/_projections.py +++ b/xrspatial/reproject/_projections.py @@ -467,6 +467,8 @@ def _aea_params(crs): return None if d.get('proj') != 'aea': return None + if not _is_wgs84_compatible_ellipsoid(crs): + return None lat_1 = math.radians(d.get('lat_1', 0.0)) lat_2 = math.radians(d.get('lat_2', lat_1)) @@ -571,6 +573,8 @@ def _cea_params(crs): return None if d.get('proj') != 'cea': return None + if not _is_wgs84_compatible_ellipsoid(crs): + return None lon_0 = math.radians(d.get('lon_0', 0.0)) lat_ts = math.radians(d.get('lat_ts', 0.0)) @@ -1757,11 +1761,35 @@ def _is_wgs84_compatible_ellipsoid(crs): Returns True for WGS84/NAD83 (no shift needed) and for datums with known Helmert parameters (NAD27, etc.) since the dispatch will wrap the projection with a datum shift. + + CRSes that define the ellipsoid through explicit axes instead of an + ``ellps``/``datum`` name are only accepted when the axes match + WGS84/GRS80. PROJ normalizes a spherical definition (``+a=R +b=R``) + to a single ``R`` key with no ``ellps``, so without the axis check a + sphere-based CRS such as the MODIS sinusoidal grid + (``+proj=sinu +a=6371007.181 +b=6371007.181``) sailed through the + empty-string match and ran the WGS84 kernels, landing ~19 km off + (GH #3275). """ try: d = _crs_to_dict(crs) except Exception: return False + # Explicit sphere radius: never WGS84-compatible. + if 'R' in d: + return False + # Explicit semi-axes must match WGS84/GRS80. The two ellipsoids + # differ by ~1e-4 m in b, so a 0.5 m tolerance keeps GRS80 in while + # rejecting every other ellipsoid (Bessel, Clarke, intl, spheres). + a_axis = d.get('a') + if a_axis is not None and abs(float(a_axis) - _WGS84_A) > 0.5: + return False + b_axis = d.get('b') + if b_axis is not None and abs(float(b_axis) - _WGS84_B) > 0.5: + return False + rf = d.get('rf') + if rf is not None and abs(float(rf) - 1.0 / _WGS84_F) > 0.01: + return False ellps = d.get('ellps', '') datum = d.get('datum', '') # WGS84 and GRS80: no shift needed diff --git a/xrspatial/tests/test_reproject_sphere_ellipsoid_guard_3275.py b/xrspatial/tests/test_reproject_sphere_ellipsoid_guard_3275.py new file mode 100644 index 000000000..677441c15 --- /dev/null +++ b/xrspatial/tests/test_reproject_sphere_ellipsoid_guard_3275.py @@ -0,0 +1,183 @@ +"""Regression tests for GH #3275. + +The numba/CUDA fast-path dispatch applied WGS84 ellipsoid constants to +CRSes on other ellipsoids: + +1. ``_is_wgs84_compatible_ellipsoid`` only inspected the ``ellps`` / + ``datum`` keys of ``crs.to_dict()`` and treated empty values as + WGS84. PROJ normalizes spherical definitions (``+a=R +b=R``) to a + bare ``R`` key with no ``ellps``, so the standard MODIS sinusoidal + grid passed the check and the WGS84 meridional-arc kernel ran + against a sphere: ~18.9 km error at mid latitudes. +2. ``_aea_params`` and ``_cea_params`` had no ellipsoid check at all, + so sphere-based AEA/CEA definitions ran the WGS84 kernels with up to + ~24 km error. + +Non-representable pairs must return None from the dispatchers and fall +back to pyproj, like the non-WGS84 datum pairs from GH #2651 already do. +""" +from __future__ import annotations + +import warnings + +import numpy as np +import pytest + +from xrspatial.utils import has_cuda_and_cupy + +try: + import pyproj + HAS_PYPROJ = True +except ImportError: + HAS_PYPROJ = False + +HAS_CUPY = has_cuda_and_cupy() + +pytestmark = pytest.mark.skipif( + not HAS_PYPROJ, reason="pyproj required for reproject tests" +) + +MODIS_SINU = ( + "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 " + "+a=6371007.181 +b=6371007.181 +units=m +no_defs" +) +SPHERE_AEA = ( + "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 " + "+a=6370997 +b=6370997 +units=m +no_defs" +) +SPHERE_CEA = "+proj=cea +lat_ts=30 +a=6371228 +b=6371228 +units=m +no_defs" + + +class TestEllipsoidCompatibilityHelper: + def test_rejects_r_defined_sphere(self): + from xrspatial.reproject._projections import _is_wgs84_compatible_ellipsoid + + # PROJ normalizes a=b to a bare R key. + modis = pyproj.CRS(MODIS_SINU) + assert 'R' in modis.to_dict() + assert not _is_wgs84_compatible_ellipsoid(modis) + + def test_rejects_named_sphere(self): + from xrspatial.reproject._projections import _is_wgs84_compatible_ellipsoid + + crs = pyproj.CRS("+proj=laea +lat_0=45 +lon_0=-100 +ellps=sphere +units=m +no_defs") + assert not _is_wgs84_compatible_ellipsoid(crs) + + def test_accepts_wgs84_and_grs80(self): + # EPSG:3857 is deliberately absent: Web Mercator is defined on a + # sphere (a=b=6378137) and is dispatched by EPSG code without + # consulting this helper, so the helper correctly reports it as + # not WGS84-ellipsoid-compatible. + from xrspatial.reproject._projections import _is_wgs84_compatible_ellipsoid + + for code in (3395, 32633, 5070, 6933, 3035, 3031): + assert _is_wgs84_compatible_ellipsoid(pyproj.CRS(code)), code + + def test_accepts_explicit_wgs84_axes(self): + from xrspatial.reproject._projections import _is_wgs84_compatible_ellipsoid + + crs = pyproj.CRS( + "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 " + "+a=6378137 +b=6356752.314245179 +units=m +no_defs" + ) + assert _is_wgs84_compatible_ellipsoid(crs) + + def test_still_accepts_helmert_datums(self): + # Datums with Helmert parameters stay "compatible": the dispatch + # guards (try_numba_transform / transform_points) bail on them + # separately via _get_datum_params. + from xrspatial.reproject._projections import _is_wgs84_compatible_ellipsoid + + crs = pyproj.CRS(27700) # OSGB36 / Airy + assert _is_wgs84_compatible_ellipsoid(crs) + + +class TestSphereCrsesBailToPyproj: + def test_modis_sinusoidal_no_fast_path(self): + from xrspatial.reproject._projections import try_numba_transform + + geo = pyproj.CRS(4326) + modis = pyproj.CRS(MODIS_SINU) + cb = (0.0, 4.0e6, 1.0e5, 4.1e6) + assert try_numba_transform(geo, modis, cb, (4, 4)) is None + assert try_numba_transform(modis, geo, (0.0, 40.0, 1.0, 41.0), (4, 4)) is None + + def test_sphere_aea_no_fast_path(self): + from xrspatial.reproject._projections import transform_points + + geo = pyproj.CRS(4326) + crs = pyproj.CRS(SPHERE_AEA) + lons = np.array([-100.0, -90.0]) + lats = np.array([30.0, 40.0]) + assert transform_points(geo, crs, lons, lats) is None + + def test_sphere_cea_no_fast_path(self): + from xrspatial.reproject._projections import transform_points + + geo = pyproj.CRS(4326) + crs = pyproj.CRS(SPHERE_CEA) + lons = np.array([-100.0, 0.0]) + lats = np.array([30.0, 45.0]) + assert transform_points(geo, crs, lons, lats) is None + + def test_wgs84_pairs_keep_fast_path(self): + from xrspatial.reproject._projections import transform_points + + geo = pyproj.CRS(4326) + for code in (3857, 3395, 32633, 5070, 6933, 3035): + res = transform_points( + geo, pyproj.CRS(code), np.array([10.0]), np.array([45.0]) + ) + assert res is not None, code + + +class TestModisReprojectEndToEnd: + def test_fast_path_matches_exact_pyproj_path(self): + """4326 -> MODIS sinusoidal must agree with transform_precision=0. + + Before the guard, the default path ran the WGS84 sinusoidal + kernel against the MODIS sphere and sampled source pixels up to + ~19 km away from the exact pyproj path. + """ + import xarray as xr + + from xrspatial.reproject import reproject + + ny, nx = 120, 160 + y = np.linspace(60, 30, ny) + x = np.linspace(-20, 20, nx) + xx, yy = np.meshgrid(x, y) + data = np.sin(np.radians(xx * 4)) * np.cos(np.radians(yy * 4)) * 100 + raster = xr.DataArray( + data, dims=['y', 'x'], coords={'y': y, 'x': x}, + attrs={'crs': 'EPSG:4326'}, + ) + with warnings.catch_warnings(): + warnings.simplefilter('ignore') + exact = reproject(raster, MODIS_SINU, transform_precision=0) + fast = reproject(raster, MODIS_SINU) + mask = np.isfinite(exact.values) & np.isfinite(fast.values) + assert mask.sum() > 0.5 * mask.size + # Residual is the control-grid interpolation only (<0.1% of the + # +/-100 data range). Pre-fix the max difference was ~150. + diff = np.abs(exact.values[mask] - fast.values[mask]) + assert float(diff.max()) < 1.0 + + +@pytest.mark.skipif(not HAS_CUPY, reason="CUDA/cupy not available") +class TestCudaSphereGuard: + def test_cuda_modis_sinusoidal_no_fast_path(self): + from xrspatial.reproject._projections_cuda import try_cuda_transform + + geo = pyproj.CRS(4326) + modis = pyproj.CRS(MODIS_SINU) + cb = (0.0, 4.0e6, 1.0e5, 4.1e6) + assert try_cuda_transform(geo, modis, cb, (4, 4)) is None + + def test_cuda_wgs84_fast_path_still_active(self): + from xrspatial.reproject._projections_cuda import try_cuda_transform + + geo = pyproj.CRS(4326) + merc = pyproj.CRS(3857) + cb = (0.0, 0.0, 1.0e5, 1.0e5) + assert try_cuda_transform(geo, merc, cb, (4, 4)) is not None From 0a8e7913352f2a0bdd23664ef2fc8af8e95ad21d Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Fri, 12 Jun 2026 08:07:13 -0700 Subject: [PATCH 3/3] Address review: pin sphere-as-source dispatch branch (#3275) --- xrspatial/tests/test_reproject_sphere_ellipsoid_guard_3275.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/xrspatial/tests/test_reproject_sphere_ellipsoid_guard_3275.py b/xrspatial/tests/test_reproject_sphere_ellipsoid_guard_3275.py index 677441c15..d7776e123 100644 --- a/xrspatial/tests/test_reproject_sphere_ellipsoid_guard_3275.py +++ b/xrspatial/tests/test_reproject_sphere_ellipsoid_guard_3275.py @@ -110,6 +110,10 @@ def test_sphere_aea_no_fast_path(self): lons = np.array([-100.0, -90.0]) lats = np.array([30.0, 40.0]) assert transform_points(geo, crs, lons, lats) is None + # Sphere as the source goes through the other dispatch branch. + assert transform_points( + crs, geo, np.array([0.0]), np.array([0.0]) + ) is None def test_sphere_cea_no_fast_path(self): from xrspatial.reproject._projections import transform_points