From 790c11298d089f41bbe757e08bee0d00caf2ddac Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Tue, 9 Jun 2026 16:46:07 -0700 Subject: [PATCH 1/2] Update accuracy sweep state for reproject (#3094, #3096) --- .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 974a9bebd..e13bddce2 100644 --- a/.claude/sweep-accuracy-state.csv +++ b/.claude/sweep-accuracy-state.csv @@ -27,7 +27,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-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-05-29,2721,MEDIUM,4;5,Bounded GREAT_CIRCLE on dask (both numpy+cupy) raised ValueError: map_overlap pad depth = max_distance/cellsize mixed metre distance with degree cellsize. numpy/cupy backends fine. Fixed by measuring per-pixel pitch with active metric (PR #2722). Cat1 float32 output is documented design choice; NaN/Inf masking via np.isfinite consistent; numpy GDAL-sweep matches exact nearest and cupy brute-force on tested grids. -reproject,2026-05-29,2620,HIGH,5,"Cat5 backend inconsistency: cupy _resample_cupy (cupyx map_coordinates) diverged from numpy/native on pyproj-fallback CRS pairs (projected->projected, e.g. EPSG:32633->3857). Edge-band cval=0.0 bleed (all modes, ~534/pixel) + cubic B-spline vs Catmull-Rom (~0.45 interior). Fixed PR for #2620: route eager+dask cupy through _resample_cupy_native. Other files clean: _merge numpy/cupy structurally identical; _datum_grids/_vertical/_itrf use -0.5 pixel-center interp and self-inequality NaN checks; WGS84/GRS80 constants correct; curvature correction n/a (no geodesic gradient here). LOW (not fixed): _transform._bilinear_interp_2ch docstring claims parallel but isn't." +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." 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 33f67e4f15d31bcd5e93e52b0452d0d89ce1669a Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Tue, 9 Jun 2026 16:48:05 -0700 Subject: [PATCH 2/2] Gate CUDA reproject fast path on WGS84-compatible datums (#3094) try_cuda_transform never got the non-WGS84 datum guard that #2651 added to the CPU fast paths. The projected CRS matchers accept any datum in the Helmert table, so cupy and dask+cupy reprojects between WGS84 and e.g. EPSG:27700 ran the WGS84 Krueger series with no datum shift and landed ~100 m off while numpy went through pyproj. Bail to the CPU path (which defers these pairs to pyproj) and add CUDA-side tests mirroring TestNonWgsDatumNumbaFastPath plus an end-to-end numpy-vs-cupy parity check. --- xrspatial/reproject/_projections_cuda.py | 18 +++++- xrspatial/tests/test_reproject.py | 73 ++++++++++++++++++++++++ 2 files changed, 88 insertions(+), 3 deletions(-) diff --git a/xrspatial/reproject/_projections_cuda.py b/xrspatial/reproject/_projections_cuda.py index e820f80ac..46d30af16 100644 --- a/xrspatial/reproject/_projections_cuda.py +++ b/xrspatial/reproject/_projections_cuda.py @@ -731,15 +731,27 @@ def try_cuda_transform(src_crs, tgt_crs, chunk_bounds, chunk_shape): import cupy as cp from ._projections import (_A_RECT, _ALPHA, _APA, _BETA, _CBG, _CGB, _MLFN_EN, _QP, - _WGS84_E2, _aea_params, _cea_params, _get_epsg, - _is_geographic_wgs84_or_nad83, _laea_params, _lcc_params, - _sinu_params, _stere_params, _tmerc_params, _utm_params) + _WGS84_E2, _aea_params, _cea_params, _get_datum_params, + _get_epsg, _is_geographic_wgs84_or_nad83, _laea_params, + _lcc_params, _sinu_params, _stere_params, _tmerc_params, + _utm_params) src_epsg = _get_epsg(src_crs) tgt_epsg = _get_epsg(tgt_crs) if src_epsg is None and tgt_epsg is None: return None + # The kernels below run in WGS84 with no datum shift. The projected + # CRS matchers (_tmerc_params, _lcc_params, ...) accept any datum in + # the Helmert table via _is_wgs84_compatible_ellipsoid, so without + # this guard a pair like EPSG:4326 <-> EPSG:27700 (OSGB36/Airy) is + # projected with WGS84 constants and lands ~100 m off. Mirror the + # CPU fast-path guard from GH #2651: bail out so the chunk worker + # falls back to the CPU path, which defers these pairs to pyproj. + if (_get_datum_params(src_crs) is not None + or _get_datum_params(tgt_crs) is not None): + return None + height, width = chunk_shape left, bottom, right, top = chunk_bounds res_x = (right - left) / width diff --git a/xrspatial/tests/test_reproject.py b/xrspatial/tests/test_reproject.py index 11d7a2829..5345def2c 100644 --- a/xrspatial/tests/test_reproject.py +++ b/xrspatial/tests/test_reproject.py @@ -6925,3 +6925,76 @@ def test_source_coords_match_pyproj_for_osgb36(self): # Guard against the old corruption: coords must be metres, not degrees. assert np.all(np.abs(src_x) > 1000.0) assert np.all(np.abs(src_y) > 1000.0) + + +@pytest.mark.skipif(not HAS_CUPY, reason="cupy required") +class TestNonWgsDatumCudaFastPath: + """The CUDA fast path must bail for non-WGS84 datums (GH #3094). + + GH #2651 gated the CPU fast paths so non-WGS84 datums fall back to + pyproj, but try_cuda_transform kept dispatching them. The projected + CRS matchers accept any datum in the Helmert table, so a pair like + EPSG:4326 <-> EPSG:27700 (OSGB36 / Airy) ran the WGS84 Krueger + series with no datum shift and returned coordinates ~100 m off, + making the cupy and dask+cupy backends diverge from numpy. + """ + + def test_cuda_fast_path_disabled_for_non_wgs_target(self): + from xrspatial.reproject._projections_cuda import try_cuda_transform + src = pyproj.CRS('EPSG:4326') + tgt = pyproj.CRS('EPSG:27700') + result = try_cuda_transform( + src, tgt, (400000.0, 200000.0, 410000.0, 210000.0), (4, 4), + ) + assert result is None + + def test_cuda_fast_path_disabled_for_projected_non_wgs_source(self): + from xrspatial.reproject._projections_cuda import try_cuda_transform + src = pyproj.CRS('EPSG:27700') + tgt = pyproj.CRS('EPSG:4326') + result = try_cuda_transform(src, tgt, (-2.0, 51.0, -1.0, 52.0), (4, 4)) + assert result is None + + def test_cuda_fast_path_disabled_for_geographic_non_wgs_source(self): + # NAD27 geographic (Clarke 1866 datum) -> Web Mercator. + from xrspatial.reproject._projections_cuda import try_cuda_transform + src = pyproj.CRS('EPSG:4267') + tgt = pyproj.CRS('EPSG:3857') + result = try_cuda_transform( + src, tgt, (-8000000.0, 4000000.0, -7900000.0, 4100000.0), (4, 4), + ) + assert result is None + + def test_cuda_wgs_fast_path_still_active(self): + # WGS84 UTM <-> WGS84 geographic must keep using the CUDA path. + from xrspatial.reproject._projections_cuda import try_cuda_transform + src = pyproj.CRS('EPSG:32617') + tgt = pyproj.CRS('EPSG:4326') + result = try_cuda_transform(src, tgt, (-84.0, 40.0, -83.0, 41.0), (4, 4)) + assert result is not None + + def test_cupy_reproject_matches_numpy_for_osgb36_target(self): + # End to end: cupy must agree with numpy for a non-WGS84 datum + # target. Before the fix the cupy backend sampled ~10 pixels away + # from the right source location (~100 m datum/ellipsoid error). + from xrspatial.reproject import reproject + rng = np.random.default_rng(3094) + data = rng.random((64, 64)) + coords = {'y': np.linspace(52.0, 51.0, 64), + 'x': np.linspace(-2.0, -1.0, 64)} + host = xr.DataArray( + data, dims=['y', 'x'], coords=coords, + attrs={'crs': 'EPSG:4326'}, + ) + eager = reproject(host, 'EPSG:27700').values + gpu = host.copy(data=cp.asarray(data)) + gpu_out = reproject(gpu, 'EPSG:27700') + gpu_arr = cp.asnumpy(gpu_out.data) + assert eager.shape == gpu_arr.shape + # NaN masks must agree cell for cell. + np.testing.assert_array_equal( + np.isfinite(eager), np.isfinite(gpu_arr), + ) + np.testing.assert_allclose( + eager, gpu_arr, rtol=1e-4, atol=1e-4, equal_nan=True, + )