From 0230d051a8c3208f75e8578b0ac27dff1e2fe3e5 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Mon, 27 Jul 2026 14:49:22 -0400 Subject: [PATCH 1/4] Stop the numba fast path claiming CRSs it cannot reproduce (#3697) The fast-path parameter extractors accepted CRS definitions whose PROJ parameters they do not model, then projected them with the wrong constants. Sweeping the EPSG database, 220 of the 2685 projected CRSs the fast path accepted came back more than 5 m from pyproj. Four causes: - _lcc_params and _aea_params defaulted lat_2 to a lat_1 that was already in radians, so math.radians ran over it twice. A 47 degree standard parallel became 0.82 degrees, the pair then looked distinct, and the two-parallel branch computed a cone constant for a parallel pair that does not exist. PROJ writes a single-parallel LCC with no lat_2 at all, so every 1SP definition hit it: 185 EPSG codes, worst 3819 m. - _aea_params, _cea_params, _laea_params, _stere_params, _sinu_params and _sterea_params never read units, so a CRS in feet was read as metres. They have no to_meter factor to hand back the way _lcc_params and _tmerc_params do, so they now decline a non-metre CRS instead. - Every extractor ignored +axis=. The kernels emit easting/northing, so a west/south grid came out negated and swapped (EPSG:2046-2055, 7085 km). - _is_wgs84_compatible_ellipsoid read a missing +datum= key as WGS84. PROJ only emits that key for datums it has an alias for, so GGRS87, Israel 1993 and VN-2000 took the fast path on their GRS80/WGS84 ellipsoid and their shift was dropped. The guard now measures the offset instead of matching names, with the threshold above the NAD83-as-WGS84 approximation this module already makes on purpose. After the fix the same EPSG-wide sweep reports 0 CRSs more than 5 m out, and the fast path still covers 2649 of them. Both dispatchers share these extractors, so this covers numpy, dask+numpy and the streaming path via try_numba_transform and cupy, dask+cupy via try_cuda_transform, plus output-bounds estimation in transform_points. --- .claude/sweep-accuracy-state.csv | 2 +- xrspatial/reproject/_projections.py | 108 +++++- .../test_reproject_fast_path_guard_3697.py | 310 ++++++++++++++++++ 3 files changed, 414 insertions(+), 6 deletions(-) create mode 100644 xrspatial/tests/test_reproject_fast_path_guard_3697.py diff --git a/.claude/sweep-accuracy-state.csv b/.claude/sweep-accuracy-state.csv index 6d2b39673..5b6f2a566 100644 --- a/.claude/sweep-accuracy-state.csv +++ b/.claude/sweep-accuracy-state.csv @@ -30,7 +30,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-07-22,3689,MEDIUM,5,"Cat5 backend divergence: allocation/direction argmin compared distances at different precision per backend: numpy brute force float32 (_distance rounds), CUDA kernel float64, cKDTree tie detection exact float64. On non-lattice coord grids (res 0.1/linspace/reprojected) geometric ties differ by ~1e-13 float64 noise but round equal at float32, so numpy applied the lowest-flat-index rule while cupy and dask kdtree each picked their own float64-closer target; 92/300 random 2-target layouts on a 40x40 res-0.1 grid had >=1 diverging pixel, and cupy vs cKDTree also disagreed with each other (squared vs rooted float64 orderings). proximity output unaffected (float32 either way). Fix #3689: round candidates to float32 in the CUDA argmin loop, detect cKDTree ties at float32, document float32 tie precision in both docstrings; 16-case cross-backend regression test on a res-0.1 grid (10 failed pre-fix). CUDA available; cupy + dask+cupy executed, 595 proximity tests pass. Cats 1-4,6 clean: matches scipy.ndimage.distance_transform_edt exactly and float64 brute force for all 3 metrics/modes/bounded+unbounded; target-pixel/translation invariants hold. Dedup honored: #3443 inclusive bound, #3392/#3555 f32-ulp bound widen, #3389 range-test precision, #3090 tie order. Test-coverage note: all prior tie-break tests used integer-lattice coords, which is why this class never fired." rasterize,2026-06-18,3384,HIGH,1;5,dask all_touched polygon-boundary supercover walk re-extracted per tile in tile-local float pixel coords; floor() tie on an on-grid boundary segment diverged from eager (Cat1 precision -> Cat5 backend split). Fixed by extracting boundary float segments in the global grid frame and shifting by integer tile offset (#3384). Verified numpy/cupy/dask+numpy/dask+cupy. -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)." +reproject,2026-07-27,3697,HIGH,4;6,"2026-07-27 sweep: filed #3697 (HIGH, cat 4;6) -- numba fast path accepted CRS definitions it does not model; 220 of the 2685 fast-path-accepted EPSG projected CRSs were >5 m from pyproj. Four causes, all fixed in PR for #3697: (a) _lcc_params/_aea_params defaulted lat_2 to an already-radians lat_1 so math.radians ran twice, breaking every 1SP LCC -- 185 codes, worst 3819 m (EPSG:8325); (b) _aea/_cea/_laea/_stere/_sinu/_sterea_params ignored units, reading feet as metres (EPSG:10481, 6594 km); (c) +axis= ignored, so wsu grids came out negated+swapped (EPSG:2046-2055, 7085 km); (d) _is_wgs84_compatible_ellipsoid read a missing +datum= key as WGS84, letting GGRS87/Israel 1993/VN-2000 skip their shift (24 codes, 325 m). After the fix the same EPSG-wide sweep reports 0 CRSs >5 m and the fast path still covers 2649. PRIOR FIXES RE-VERIFIED STILL HOLDING: #3274 (LAEA inverse rq, authalic series), #3275 (sphere/ellipsoid guard), #3276 (itrf ppm). VERIFIED CLEAN this run: merc/emerc/UTM/tmerc(incl. Zb lat_0!=0)/LCC 2SP/AEA/CEA/sinu/LAEA(obliq,equit,N/S pole)/polar stere(lat_ts, k_0, plain)/sterea forward+inverse vs pyproj <=1e-5 m, except the documented ~1.6 mm authalic truncation (cea/laea inverse) and the deliberate NAD83-as-WGS84 approximation (0.3-1.8 m on EPSG:2260/3978/5070, documented at _utm_params). Resampling kernels clean on all 4 backends: constant-raster partition of unity (<=4.4e-15), same-CRS nearest identity bitwise exact, no NaN leak into fully-valid stencils, numpy vs cupy vs dask+numpy vs dask+cupy agree to <=1.2e-13 for nearest/bilinear/cubic. KNOWN LOW, documented only (not fixed): _source_footprint_in_target probe array typo uses x-midpoint mx as a latitude in the last 3 ys entries (__init__.py ~2005) -- bbox stays a superset, correctness unaffected. NEW LOW, documented only: _place_same_crs (__init__.py ~2660) direct-copies a tile whose resolution differs from the output by up to 1%, so a wide tile drifts up to 0.01*width pixels; merge only. TEST-COVERAGE GAP for the test-coverage sweep: general_checks.py has no helper that exercises the projection param extractors, so every fast-path acceptance bug (#2651, #3275, #3697) reached main with green tests; a pyproj-differential test over a CRS list would have caught all three. gdal available but not used -- pyproj is the authoritative reference for the projection kernels and the resampling kernels were checked by invariant instead. cuda-available: cupy and dask+cupy 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-07-03,3626,MEDIUM,4,"Re-audit post-#1407 fix (ground distance OK). MEDIUM #3626: ray azimuths uniform in cell-INDEX angle not ground azimuth, so anisotropic cellsizes bias the azimuthal quadrature (45-deg ramp, true SVF 0.75: 0.728/0.773/0.693 at 1:1/1:2/2:1 aspect; +0.105 pure weighting error at 4:1); fixed by aiming rays in ground space (cos(phi)/cx, sin(phi)/cy normalized), square-cell results bit-identical. Cats 1-3,5 clean: float64 throughout, NaN center->NaN out, flat=1.0 exactly all 4 backends, flip/rot90/translate invariants at machine eps, numpy==dask==cupy==dask+cupy (max 2.2e-16) incl NaN holes straddling chunk boundaries, map_overlap depth==max_radius==stencil radius, depth>chunksize works (dask 2025.7). Cat6 skipped: richdem-unavailable rvt-unavailable grass-unavailable. LOWs (not fixed): NaN neighbor silently truncates ray (horizon beyond nodata ignored, undocumented but keeps numpy/dask edge semantics consistent); nearest-cell ray rounding gives ~-2.9% jitter bias on 45-deg slopes even square cells (inherent to integer-cell marching, reference tools interpolate along ray). For test-coverage sweep: cross-backend tests had no NaN-input and no anisotropic-cellsize cases (added aniso numpy/dask/cupy parity in #3626 PR); for style sweep: unused import _boundary_to_dask (F401)." diff --git a/xrspatial/reproject/_projections.py b/xrspatial/reproject/_projections.py index 44648cf39..e335b7082 100644 --- a/xrspatial/reproject/_projections.py +++ b/xrspatial/reproject/_projections.py @@ -24,6 +24,7 @@ import math import threading import warnings +from functools import lru_cache import numpy as np from numba import njit, prange @@ -156,6 +157,66 @@ def _crs_to_dict(crs): return crs.to_dict() +# Largest datum offset from WGS84 the fast-path kernels will absorb. The +# kernels project on the WGS84 ellipsoid and apply no datum shift, so a CRS +# whose datum sits further away than this has to go through pyproj instead. +# NAD83 and ETRS89 land 1-2 m from WGS84 and are deliberately treated as +# WGS84 here (see _is_supported_geographic), so the bar sits above them. +# Every datum this guard has to catch is far above it: GGRS87 is 325 m out, +# VN-2000 226 m, Israel 1993 79 m. See GH #3697. +_MAX_DATUM_OFFSET_M = 5.0 + + +@lru_cache(maxsize=256) +def _datum_offset_from_wgs84(crs): + """Metres between *crs*'s datum and WGS84, measured at its own origin. + + ``to_dict()`` only carries a ``+datum=`` key for datums PROJ has an + alias for. Everything else -- GGRS87, Israel 1993, VN-2000, + Hartebeesthoek94 -- arrives with the key missing, which the name test + in :func:`_is_wgs84_compatible_ellipsoid` used to read as "WGS84". + Measuring the offset is the only test that does not depend on PROJ + having a name for the datum. + + Returns ``None`` when the offset cannot be measured (pyproj absent, a + lite CRS with no geodetic base, a failed transform). Callers treat + that as "no evidence of a shift" and keep their existing verdict. + """ + geodetic = getattr(crs, 'geodetic_crs', None) + if geodetic is None: + return None + try: + import pyproj + except ImportError: + return None + try: + d = _crs_to_dict(crs) + lon = float(d.get('lon_0', 0.0)) + lat = float(d.get('lat_0', d.get('lat_1', 0.0))) + transformer = pyproj.Transformer.from_crs( + geodetic, pyproj.CRS.from_epsg(4326), always_xy=True, + ) + lon_w, lat_w = transformer.transform(lon, lat) + if not (math.isfinite(lon_w) and math.isfinite(lat_w)): + return None + _, _, dist = pyproj.Geod(ellps='WGS84').inv(lon, lat, lon_w, lat_w) + return abs(float(dist)) + except Exception: + return None + + +def _is_metric(d): + """True when a CRS's linear unit is the metre. + + ``_lcc_params`` and ``_tmerc_params`` hand back a ``to_meter`` factor + and their callers apply it. The other extractors build metre- + denominated constants with nowhere to put such a factor, so a CRS in + feet has to fall back to pyproj rather than be read as metres. See + GH #3697. + """ + return d.get('units', 'm') == 'm' + + def _get_datum_params(crs): """Return (dx, dy, dz, rx, ry, rz, ds, a_src, f_src) for a non-WGS84 datum. @@ -338,8 +399,16 @@ def _lcc_params(crs): if to_meter is None: return None - lat_1 = math.radians(d.get('lat_1', d.get('lat_0', 0.0))) - lat_2 = math.radians(d.get('lat_2', lat_1)) + # Take the degree values first, then convert. Defaulting lat_2 to the + # already-radians lat_1 ran math.radians over it twice, turning a 47 + # degree standard parallel into 0.82 degrees; the pair then looked + # distinct and the 2SP branch below computed a cone constant for a + # parallel pair that does not exist. PROJ writes a single-parallel LCC + # with no lat_2 at all, so every 1SP definition hit it (GH #3697). + lat_1_deg = d.get('lat_1', d.get('lat_0', 0.0)) + lat_2_deg = d.get('lat_2', lat_1_deg) + lat_1 = math.radians(lat_1_deg) + lat_2 = math.radians(lat_2_deg) lat_0 = math.radians(d.get('lat_0', 0.0)) lon_0 = math.radians(d.get('lon_0', 0.0)) k0_param = d.get('k_0', d.get('k', 1.0)) @@ -476,9 +545,15 @@ def _aea_params(crs): return None if not _is_wgs84_compatible_ellipsoid(crs): return None + if not _is_metric(d): + return None - lat_1 = math.radians(d.get('lat_1', 0.0)) - lat_2 = math.radians(d.get('lat_2', lat_1)) + # Degrees first, then convert -- see the matching note in _lcc_params + # for what the double conversion did to a 1SP definition (GH #3697). + lat_1_deg = d.get('lat_1', 0.0) + lat_2_deg = d.get('lat_2', lat_1_deg) + lat_1 = math.radians(lat_1_deg) + lat_2 = math.radians(lat_2_deg) lat_0 = math.radians(d.get('lat_0', 0.0)) lon_0 = math.radians(d.get('lon_0', 0.0)) @@ -582,6 +657,8 @@ def _cea_params(crs): return None if not _is_wgs84_compatible_ellipsoid(crs): return None + if not _is_metric(d): + return None lon_0 = math.radians(d.get('lon_0', 0.0)) lat_ts = math.radians(d.get('lat_ts', 0.0)) @@ -703,6 +780,8 @@ def _sinu_params(crs): return None if not _is_wgs84_compatible_ellipsoid(crs): return None + if not _is_metric(d): + return None lon_0 = math.radians(d.get('lon_0', 0.0)) fe = d.get('x_0', 0.0) fn = d.get('y_0', 0.0) @@ -769,6 +848,8 @@ def _laea_params(crs): return None if not _is_wgs84_compatible_ellipsoid(crs): return None + if not _is_metric(d): + return None lon_0 = math.radians(d.get('lon_0', 0.0)) lat_0 = math.radians(d.get('lat_0', 0.0)) @@ -963,6 +1044,8 @@ def _stere_params(crs): return None if not _is_wgs84_compatible_ellipsoid(crs): return None + if not _is_metric(d): + return None lat_0 = d.get('lat_0', 0.0) if abs(abs(lat_0) - 90.0) > 1e-6: @@ -1103,6 +1186,8 @@ def _sterea_params(crs): return None if not _is_wgs84_compatible_ellipsoid(crs): return None + if not _is_metric(d): + return None lat_0 = math.radians(d.get('lat_0', 0.0)) lon_0 = math.radians(d.get('lon_0', 0.0)) @@ -1792,6 +1877,13 @@ def _is_wgs84_compatible_ellipsoid(crs): d = _crs_to_dict(crs) except Exception: return False + # Axis order and direction. Every kernel here emits (easting, northing) + # in that order. A CRS that stores (west, south) -- the South African + # Lo grids, +axis=wsu -- or (north, east) needs the components negated + # and swapped, which none of them do, so hand it to pyproj (GH #3697). + axis = d.get('axis') + if axis is not None and axis != 'enu': + return False # Explicit sphere radius: never WGS84-compatible. if 'R' in d: return False @@ -1809,9 +1901,15 @@ def _is_wgs84_compatible_ellipsoid(crs): return False ellps = d.get('ellps', '') datum = d.get('datum', '') - # WGS84 and GRS80: no shift needed + # WGS84 and GRS80: no shift needed -- but only once we have checked + # that the datum really is WGS84-aligned. An unaliased datum reaches + # here with datum == '', so the name test alone let GGRS87, Israel + # 1993 and VN-2000 through on their GRS80/WGS84 ellipsoid (GH #3697). if (ellps in ('WGS84', 'GRS80', '') and datum in ('WGS84', 'NAD83', '')): + offset = _datum_offset_from_wgs84(crs) + if offset is not None and offset > _MAX_DATUM_OFFSET_M: + return False return True # Check if we have Helmert parameters for this datum key = datum if datum in _DATUM_PARAMS else ellps diff --git a/xrspatial/tests/test_reproject_fast_path_guard_3697.py b/xrspatial/tests/test_reproject_fast_path_guard_3697.py new file mode 100644 index 000000000..4e8a4c1c4 --- /dev/null +++ b/xrspatial/tests/test_reproject_fast_path_guard_3697.py @@ -0,0 +1,310 @@ +"""Regression tests for GH #3697. + +The numba fast path in ``xrspatial/reproject/_projections.py`` used to claim +CRS definitions whose PROJ parameters it does not model, and then project +them with the wrong constants. Four separate defects: + +1. ``_lcc_params`` / ``_aea_params`` defaulted ``lat_2`` to a ``lat_1`` that + had already been converted to radians, so ``math.radians`` ran twice and + a one-standard-parallel definition got a cone constant for a parallel + pair that does not exist. +2. ``_aea_params``, ``_cea_params``, ``_laea_params``, ``_stere_params``, + ``_sinu_params`` and ``_sterea_params`` ignored ``units``, reading a CRS + in feet as metres. +3. Every extractor ignored ``+axis=``, so a west/south oriented grid came + out negated and swapped. +4. ``_is_wgs84_compatible_ellipsoid`` read a missing ``+datum=`` key as + "WGS84", letting datums that only share the GRS80/WGS84 ellipsoid skip + their shift. +""" +import math + +import numpy as np +import pytest +import xarray as xr + +from xrspatial.reproject import _projections as P +from xrspatial.reproject import reproject + +pyproj = pytest.importorskip('pyproj') + + +WGS84 = 'EPSG:4326' + + +def _lon_field_error(crs_spec, lon_c, lat_c, span=0.5, n=401, out=32): + """Max georeferencing error, in metres, of a reprojected raster. + + The source stores its own longitude as the pixel value and the output + grid is pinned explicitly, so every output pixel should carry the + longitude pyproj assigns to that pixel's target coordinates. The source + is a smooth ramp sampled far finer than the output, so what is left is + georeferencing error rather than resampling error. + """ + tgt = pyproj.CRS.from_user_input(crs_spec) + src_crs = pyproj.CRS.from_epsg(4326) + fwd = pyproj.Transformer.from_crs(src_crs, tgt, always_xy=True) + inv = pyproj.Transformer.from_crs(tgt, src_crs, always_xy=True) + + lons = np.linspace(lon_c - span, lon_c + span, n) + lats = np.linspace(lat_c + span, lat_c - span, n) + src = xr.DataArray( + np.tile(lons, (n, 1)), dims=['y', 'x'], + coords={'y': lats, 'x': lons}, attrs={'crs': WGS84}, + ) + + cx, cy = fwd.transform(lon_c, lat_c) + ex, ey = fwd.transform( + np.array([lon_c - span / 2, lon_c + span / 2]), + np.array([lat_c - span / 2, lat_c + span / 2]), + ) + half = max(abs(ex[1] - ex[0]), abs(ey[1] - ey[0])) / 2 + result = reproject( + src, tgt, bounds=(cx - half, cy - half, cx + half, cy + half), + width=out, height=out, resampling='bilinear', + ) + + gx, gy = np.meshgrid(result.coords['x'].values, result.coords['y'].values) + true_lon, _ = inv.transform(gx.ravel(), gy.ravel()) + got_lon = np.asarray(result.data).ravel() + mask = np.isfinite(got_lon) & np.isfinite(true_lon) + assert mask.any(), f"no valid output pixels for {crs_spec}" + dlon = np.abs(got_lon[mask] - np.asarray(true_lon)[mask]) + return float(np.max(dlon) * 111320.0 * math.cos(math.radians(lat_c))) + + +# --------------------------------------------------------------------------- +# 1. lat_2 fallback +# --------------------------------------------------------------------------- + +class TestLccSingleStandardParallel: + """A 1SP Lambert Conformal Conic must use n = sin(lat_1).""" + + def test_cone_constant_is_sin_lat1(self): + # +proj=lcc +lat_1=45 +lat_0=45 with no lat_2 at all. + crs = pyproj.CRS.from_user_input( + '+proj=lcc +lat_1=45 +lat_0=45 +lon_0=-100 +k_0=1 ' + '+x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs' + ) + assert 'lat_2' not in P._crs_to_dict(crs) + params = P._lcc_params(crs) + assert params is not None + n = params[1] + assert n == pytest.approx(math.sin(math.radians(45)), abs=1e-12) + + @pytest.mark.parametrize('epsg,lon,lat', [ + ('EPSG:8325', -120.25, 44.75), # NAD83(2011) / Oregon Mitchell zone + ('EPSG:9549', 6.82, 45.18), # LTF2004(C) + ]) + def test_reproject_matches_pyproj(self, epsg, lon, lat): + assert _lon_field_error(epsg, lon, lat) < 1.0 + + def test_two_parallel_lcc_still_uses_both(self): + # EPSG:2154 carries lat_1 and lat_2, so the 2SP branch must run and + # produce a cone constant strictly between the two parallels. + crs = pyproj.CRS.from_epsg(2154) + d = P._crs_to_dict(crs) + assert d['lat_1'] != d['lat_2'] + n = P._lcc_params(crs)[1] + assert math.sin(math.radians(d['lat_2'])) < n + assert n < math.sin(math.radians(d['lat_1'])) + assert _lon_field_error('EPSG:2154', 3.0, 46.5) < 1.0 + + def test_aea_single_standard_parallel(self): + crs = pyproj.CRS.from_user_input( + '+proj=aea +lat_1=45 +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 ' + '+datum=WGS84 +units=m +no_defs' + ) + d = P._crs_to_dict(crs) + if 'lat_2' in d: + pytest.skip('this PROJ build emits an explicit lat_2') + n = P._aea_params(crs)[1] + assert n == pytest.approx(math.sin(math.radians(45)), abs=1e-12) + + +# --------------------------------------------------------------------------- +# 2. non-metre units +# --------------------------------------------------------------------------- + +class TestNonMetricUnitsFallBack: + """Extractors with no to_meter factor must decline a non-metre CRS.""" + + @pytest.mark.parametrize('proj,extractor', [ + ('+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96', + '_aea_params'), + ('+proj=cea +lon_0=0 +lat_ts=30', '_cea_params'), + ('+proj=laea +lat_0=45 +lon_0=-100', '_laea_params'), + ('+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45', '_stere_params'), + ('+proj=sinu +lon_0=0', '_sinu_params'), + ('+proj=sterea +lat_0=52 +lon_0=5 +k=0.9999', '_sterea_params'), + ]) + def test_feet_declined_metres_accepted(self, proj, extractor): + fn = getattr(P, extractor) + base = proj + ' +x_0=0 +y_0=0 +datum=WGS84 +no_defs' + metric = pyproj.CRS.from_user_input(base + ' +units=m') + feet = pyproj.CRS.from_user_input(base + ' +units=us-ft') + assert fn(metric) is not None, 'metre CRS must keep the fast path' + assert fn(feet) is None, 'foot CRS must fall back to pyproj' + + def test_reproject_us_survey_foot_albers(self): + # EPSG:10481 is +proj=aea +units=us-ft; it used to land 6594 km out. + assert _lon_field_error('EPSG:10481', -100.0, 31.25) < 1.0 + + def test_lcc_and_tmerc_keep_their_foot_support(self): + # These two do return a to_meter factor, so feet stay on the fast + # path. EPSG:2260 is tmerc in US survey feet, EPSG:2261 is lcc. + assert P._tmerc_params(pyproj.CRS.from_epsg(2260)) is not None + lcc_ft = pyproj.CRS.from_user_input( + '+proj=lcc +lat_1=41.03 +lat_2=40.66 +lat_0=40.16 +lon_0=-74 ' + '+x_0=300000 +y_0=0 +datum=WGS84 +units=us-ft +no_defs' + ) + assert P._lcc_params(lcc_ft) is not None + + +# --------------------------------------------------------------------------- +# 3. axis order / direction +# --------------------------------------------------------------------------- + +class TestAxisOrientation: + """A grid that is not east/north must not take the fast path.""" + + def test_wsu_declined(self): + # EPSG:2048 is the South African Lo19 grid: +axis=wsu. + crs = pyproj.CRS.from_epsg(2048) + assert P._crs_to_dict(crs).get('axis') == 'wsu' + assert P._tmerc_params(crs) is None + assert not P._is_wgs84_compatible_ellipsoid(crs) + + def test_enu_still_accepted(self): + crs = pyproj.CRS.from_epsg(32633) + assert P._is_wgs84_compatible_ellipsoid(crs) + + def test_reproject_south_oriented_grid(self): + assert _lon_field_error('EPSG:2048', 19.0, -29.0) < 1.0 + + +# --------------------------------------------------------------------------- +# 4. datum sharing the GRS80 / WGS84 ellipsoid +# --------------------------------------------------------------------------- + +class TestUnaliasedDatumGuard: + """A datum PROJ has no +datum= alias for must still be measured.""" + + # Lower bounds on the measured offset, well under the real values (the + # offset varies across each area of use) but far above + # _MAX_DATUM_OFFSET_M so the assertion pins the behaviour, not a + # PROJ-version-specific number. + @pytest.mark.parametrize('epsg,min_offset_m', [ + (2100, 200.0), # GGRS87 / Greek Grid + (2039, 50.0), # Israel 1993 / Israeli TM Grid + (9207, 100.0), # VN-2000 / TM-3 104-30 + ]) + def test_shifted_datum_declined(self, epsg, min_offset_m): + crs = pyproj.CRS.from_epsg(epsg) + d = P._crs_to_dict(crs) + # The dict carries no datum name, which is exactly what used to + # make the old name test read these as WGS84. + assert d.get('datum', '') == '' + assert d.get('ellps') in ('GRS80', 'WGS84') + offset = P._datum_offset_from_wgs84(crs) + assert offset is not None + assert offset > min_offset_m + assert not P._is_wgs84_compatible_ellipsoid(crs) + + @pytest.mark.parametrize('epsg', [ + 32633, # WGS 84 / UTM zone 33N + 2154, # RGF93 / Lambert-93 (ETRS89-realisation, ~0 m) + 3035, # ETRS89-extended / LAEA Europe + 5070, # NAD83 / Conus Albers (deliberately treated as WGS84) + 6350, # NAD83(2011) / Conus Albers + ]) + def test_aligned_datum_keeps_fast_path(self, epsg): + crs = pyproj.CRS.from_epsg(epsg) + assert P._is_wgs84_compatible_ellipsoid(crs) + offset = P._datum_offset_from_wgs84(crs) + if offset is not None: + assert offset <= P._MAX_DATUM_OFFSET_M + + @pytest.mark.parametrize('epsg,lon,lat', [ + (2100, 24.0, 38.0), + (2039, 35.2, 31.7), + ]) + def test_reproject_matches_pyproj(self, epsg, lon, lat): + assert _lon_field_error(f'EPSG:{epsg}', lon, lat) < 1.0 + + def test_lite_crs_has_no_geodetic_base(self): + # LiteCRS carries no geodetic_crs, so the measurement returns None + # and the curated built-in table keeps its verdict. + from xrspatial.reproject._lite_crs import CRS as LiteCRS + assert P._datum_offset_from_wgs84(LiteCRS(3857)) is None + + +# --------------------------------------------------------------------------- +# Backend parity: the guard runs on the CPU and CUDA dispatchers alike +# --------------------------------------------------------------------------- + +class TestBackendParity: + """Both dispatchers share the extractors, so both must decline.""" + + @pytest.mark.parametrize('epsg', [2100, 2048, 10481]) + def test_cpu_dispatcher_declines(self, epsg): + tgt = pyproj.CRS.from_epsg(epsg) + src = pyproj.CRS.from_epsg(4326) + fwd = pyproj.Transformer.from_crs(src, tgt, always_xy=True) + aou = tgt.area_of_use + cx, cy = fwd.transform((aou.west + aou.east) / 2, + (aou.south + aou.north) / 2) + bounds = (cx - 1000, cy - 1000, cx + 1000, cy + 1000) + assert P.try_numba_transform(src, tgt, bounds, (8, 8)) is None + + @pytest.mark.parametrize('epsg', [2100, 2048, 10481]) + def test_cuda_dispatcher_declines(self, epsg): + cuda = pytest.importorskip('numba.cuda') + if not cuda.is_available(): + pytest.skip('CUDA not available') + pytest.importorskip('cupy') + from xrspatial.reproject._projections_cuda import try_cuda_transform + + tgt = pyproj.CRS.from_epsg(epsg) + src = pyproj.CRS.from_epsg(4326) + fwd = pyproj.Transformer.from_crs(src, tgt, always_xy=True) + aou = tgt.area_of_use + cx, cy = fwd.transform((aou.west + aou.east) / 2, + (aou.south + aou.north) / 2) + bounds = (cx - 1000, cy - 1000, cx + 1000, cy + 1000) + assert try_cuda_transform(src, tgt, bounds, (8, 8)) is None + + @pytest.mark.parametrize('epsg,lon,lat', [ + (8325, -120.25, 44.75), + (2100, 24.0, 38.0), + ]) + def test_dask_matches_numpy(self, epsg, lon, lat): + da = pytest.importorskip('dask.array') + tgt = pyproj.CRS.from_epsg(epsg) + src_crs = pyproj.CRS.from_epsg(4326) + fwd = pyproj.Transformer.from_crs(src_crs, tgt, always_xy=True) + + n, span = 201, 0.4 + lons = np.linspace(lon - span, lon + span, n) + lats = np.linspace(lat + span, lat - span, n) + values = np.tile(lons, (n, 1)) + coords = {'y': lats, 'x': lons} + + cx, cy = fwd.transform(lon, lat) + bounds = (cx - 5000, cy - 5000, cx + 5000, cy + 5000) + + eager = reproject( + xr.DataArray(values, dims=['y', 'x'], coords=coords, + attrs={'crs': WGS84}), + tgt, bounds=bounds, width=24, height=24, resampling='bilinear', + ) + lazy = reproject( + xr.DataArray(da.from_array(values, chunks=(64, 64)), + dims=['y', 'x'], coords=coords, + attrs={'crs': WGS84}), + tgt, bounds=bounds, width=24, height=24, resampling='bilinear', + ) + np.testing.assert_allclose( + np.asarray(eager.data), np.asarray(lazy.data.compute()), + rtol=0, atol=1e-9, + ) From b7d5789a1bde82971e084d5fb6f53e78319a27a8 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Mon, 27 Jul 2026 14:55:24 -0400 Subject: [PATCH 2/4] Address review findings on the fast-path guard (#3697) - Fail closed when the datum offset cannot be measured. Returning None on an exception meant an unverifiable datum kept the fast path, which is the wrong direction for a correctness guard; it now returns inf so the caller sends the CRS to pyproj. None is still returned when there is genuinely nothing to measure (no pyproj, or a lite CRS with no geodetic base, where the curated built-in table already decides). - Write down why the axis check reads the PROJ dict rather than axis_info. Only axis direction matters, because always_xy=True already normalizes ordering; EPSG:2193 and EPSG:3346 are north/east in WKT and reproject correctly today. Widening the check would break them, so the comment says not to, and a test pins it. - Say in _is_wgs84_compatible_ellipsoid's docstring that it is the general fast-path precondition, not only an ellipsoid test. Kept the name: nine call sites and the #3275 tests reference it, and renaming buys nothing the docstring does not. - Note in _datum_offset_from_wgs84 that the lon_0/lat_0 probe point is deliberately approximate, since a Helmert shift varies slowly enough that an off-centre reading still lands far from the threshold. --- xrspatial/reproject/_projections.py | 48 +++++++++++++++---- .../test_reproject_fast_path_guard_3697.py | 29 +++++++++++ 2 files changed, 67 insertions(+), 10 deletions(-) diff --git a/xrspatial/reproject/_projections.py b/xrspatial/reproject/_projections.py index e335b7082..3eacadafb 100644 --- a/xrspatial/reproject/_projections.py +++ b/xrspatial/reproject/_projections.py @@ -178,9 +178,19 @@ def _datum_offset_from_wgs84(crs): Measuring the offset is the only test that does not depend on PROJ having a name for the datum. - Returns ``None`` when the offset cannot be measured (pyproj absent, a - lite CRS with no geodetic base, a failed transform). Callers treat - that as "no evidence of a shift" and keep their existing verdict. + The probe point is the CRS's own ``lon_0``/``lat_0``, which for a few + definitions sits outside the declared area of use. That is deliberate + and good enough: a Helmert shift varies slowly across a datum's extent, + so the reading is representative even off-centre, and the caller only + compares it against a threshold two orders of magnitude below the + shifts it has to catch. + + Returns ``None`` only when there is nothing to measure -- pyproj is + absent, or *crs* is a lite CRS with no geodetic base, in which case the + curated built-in table already decides. Returns ``inf`` when a + measurement was possible but failed, so the caller rejects the fast + path rather than assuming no shift: an unverifiable datum should go to + pyproj, not to kernels that silently skip its shift. """ geodetic = getattr(crs, 'geodetic_crs', None) if geodetic is None: @@ -198,11 +208,11 @@ def _datum_offset_from_wgs84(crs): ) lon_w, lat_w = transformer.transform(lon, lat) if not (math.isfinite(lon_w) and math.isfinite(lat_w)): - return None + return math.inf _, _, dist = pyproj.Geod(ellps='WGS84').inv(lon, lat, lon_w, lat_w) return abs(float(dist)) except Exception: - return None + return math.inf def _is_metric(d): @@ -1858,7 +1868,17 @@ def _is_supported_geographic(epsg): def _is_wgs84_compatible_ellipsoid(crs): - """True if *crs* uses WGS84/GRS80 OR a datum we can Helmert-shift. + """Whether the WGS84 kernels in this module can reproduce *crs*. + + Despite the name this is the general fast-path precondition, not just + an ellipsoid test. Every ``*_params`` extractor calls it, and it says + no on three separate grounds: an ellipsoid that is not WGS84/GRS80, an + axis direction that is not east/north, and a datum far enough from + WGS84 that skipping its shift would move pixels. Anything it rejects + falls back to pyproj. Keep new "can the kernels handle this?" checks + here rather than scattering them across the extractors (GH #3697). + + True if *crs* uses WGS84/GRS80 OR a datum we can Helmert-shift. Returns True for WGS84/NAD83 (no shift needed) and for datums with known Helmert parameters (NAD27, etc.) since the dispatch @@ -1877,10 +1897,18 @@ def _is_wgs84_compatible_ellipsoid(crs): d = _crs_to_dict(crs) except Exception: return False - # Axis order and direction. Every kernel here emits (easting, northing) - # in that order. A CRS that stores (west, south) -- the South African - # Lo grids, +axis=wsu -- or (north, east) needs the components negated - # and swapped, which none of them do, so hand it to pyproj (GH #3697). + # Axis direction. Every kernel here emits easting and northing, so a + # CRS that counts west or south positive -- the South African Lo grids, + # +axis=wsu -- comes out negated. Hand those to pyproj (GH #3697). + # + # Only direction matters, not axis *order*. Every transformer in this + # package is built with always_xy=True, which normalizes a northing- + # first CRS to easting-first for us; EPSG:2193 and EPSG:3346 are + # north/east in WKT and reproject correctly today. Reading `axis` from + # the PROJ dict is therefore the right test even though PROJ omits the + # key for those two: it omits it exactly when order is the only + # difference. Do NOT widen this to crs.axis_info -- that would reject + # both of them for a difference always_xy has already handled. axis = d.get('axis') if axis is not None and axis != 'enu': return False diff --git a/xrspatial/tests/test_reproject_fast_path_guard_3697.py b/xrspatial/tests/test_reproject_fast_path_guard_3697.py index 4e8a4c1c4..55892d653 100644 --- a/xrspatial/tests/test_reproject_fast_path_guard_3697.py +++ b/xrspatial/tests/test_reproject_fast_path_guard_3697.py @@ -238,6 +238,35 @@ def test_lite_crs_has_no_geodetic_base(self): from xrspatial.reproject._lite_crs import CRS as LiteCRS assert P._datum_offset_from_wgs84(LiteCRS(3857)) is None + def test_measurement_failure_fails_closed(self, monkeypatch): + # A CRS we could have measured but could not must be rejected, not + # waved through -- an unverifiable datum belongs on the pyproj path. + crs = pyproj.CRS.from_epsg(32633) + assert P._is_wgs84_compatible_ellipsoid(crs) + + def _boom(*args, **kwargs): + raise RuntimeError('transform unavailable') + + P._datum_offset_from_wgs84.cache_clear() + monkeypatch.setattr(pyproj.Transformer, 'from_crs', _boom) + try: + assert P._datum_offset_from_wgs84(crs) == math.inf + assert not P._is_wgs84_compatible_ellipsoid(crs) + finally: + P._datum_offset_from_wgs84.cache_clear() + + def test_axis_order_only_crs_keeps_fast_path(self): + # EPSG:2193 and EPSG:3346 are north/east in WKT but reproject + # correctly because always_xy=True normalizes the ordering. The + # axis guard must not reject them; see the note in + # _is_wgs84_compatible_ellipsoid. + for epsg in (2193, 3346): + crs = pyproj.CRS.from_epsg(epsg) + assert [a.direction for a in crs.axis_info] == ['north', 'east'] + assert P._crs_to_dict(crs).get('axis') is None + assert P._is_wgs84_compatible_ellipsoid(crs) + assert P._tmerc_params(crs) is not None + # --------------------------------------------------------------------------- # Backend parity: the guard runs on the CPU and CUDA dispatchers alike From 61901226ba8b6ca2309374db49b8f30185684e68 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Mon, 27 Jul 2026 14:57:31 -0400 Subject: [PATCH 3/4] Record the verified numbers for the reproject LOW findings (#3697) --- .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 5b6f2a566..47b4bc4f9 100644 --- a/.claude/sweep-accuracy-state.csv +++ b/.claude/sweep-accuracy-state.csv @@ -30,7 +30,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-07-22,3689,MEDIUM,5,"Cat5 backend divergence: allocation/direction argmin compared distances at different precision per backend: numpy brute force float32 (_distance rounds), CUDA kernel float64, cKDTree tie detection exact float64. On non-lattice coord grids (res 0.1/linspace/reprojected) geometric ties differ by ~1e-13 float64 noise but round equal at float32, so numpy applied the lowest-flat-index rule while cupy and dask kdtree each picked their own float64-closer target; 92/300 random 2-target layouts on a 40x40 res-0.1 grid had >=1 diverging pixel, and cupy vs cKDTree also disagreed with each other (squared vs rooted float64 orderings). proximity output unaffected (float32 either way). Fix #3689: round candidates to float32 in the CUDA argmin loop, detect cKDTree ties at float32, document float32 tie precision in both docstrings; 16-case cross-backend regression test on a res-0.1 grid (10 failed pre-fix). CUDA available; cupy + dask+cupy executed, 595 proximity tests pass. Cats 1-4,6 clean: matches scipy.ndimage.distance_transform_edt exactly and float64 brute force for all 3 metrics/modes/bounded+unbounded; target-pixel/translation invariants hold. Dedup honored: #3443 inclusive bound, #3392/#3555 f32-ulp bound widen, #3389 range-test precision, #3090 tie order. Test-coverage note: all prior tie-break tests used integer-lattice coords, which is why this class never fired." rasterize,2026-06-18,3384,HIGH,1;5,dask all_touched polygon-boundary supercover walk re-extracted per tile in tile-local float pixel coords; floor() tie on an on-grid boundary segment diverged from eager (Cat1 precision -> Cat5 backend split). Fixed by extracting boundary float segments in the global grid frame and shifting by integer tile offset (#3384). Verified numpy/cupy/dask+numpy/dask+cupy. -reproject,2026-07-27,3697,HIGH,4;6,"2026-07-27 sweep: filed #3697 (HIGH, cat 4;6) -- numba fast path accepted CRS definitions it does not model; 220 of the 2685 fast-path-accepted EPSG projected CRSs were >5 m from pyproj. Four causes, all fixed in PR for #3697: (a) _lcc_params/_aea_params defaulted lat_2 to an already-radians lat_1 so math.radians ran twice, breaking every 1SP LCC -- 185 codes, worst 3819 m (EPSG:8325); (b) _aea/_cea/_laea/_stere/_sinu/_sterea_params ignored units, reading feet as metres (EPSG:10481, 6594 km); (c) +axis= ignored, so wsu grids came out negated+swapped (EPSG:2046-2055, 7085 km); (d) _is_wgs84_compatible_ellipsoid read a missing +datum= key as WGS84, letting GGRS87/Israel 1993/VN-2000 skip their shift (24 codes, 325 m). After the fix the same EPSG-wide sweep reports 0 CRSs >5 m and the fast path still covers 2649. PRIOR FIXES RE-VERIFIED STILL HOLDING: #3274 (LAEA inverse rq, authalic series), #3275 (sphere/ellipsoid guard), #3276 (itrf ppm). VERIFIED CLEAN this run: merc/emerc/UTM/tmerc(incl. Zb lat_0!=0)/LCC 2SP/AEA/CEA/sinu/LAEA(obliq,equit,N/S pole)/polar stere(lat_ts, k_0, plain)/sterea forward+inverse vs pyproj <=1e-5 m, except the documented ~1.6 mm authalic truncation (cea/laea inverse) and the deliberate NAD83-as-WGS84 approximation (0.3-1.8 m on EPSG:2260/3978/5070, documented at _utm_params). Resampling kernels clean on all 4 backends: constant-raster partition of unity (<=4.4e-15), same-CRS nearest identity bitwise exact, no NaN leak into fully-valid stencils, numpy vs cupy vs dask+numpy vs dask+cupy agree to <=1.2e-13 for nearest/bilinear/cubic. KNOWN LOW, documented only (not fixed): _source_footprint_in_target probe array typo uses x-midpoint mx as a latitude in the last 3 ys entries (__init__.py ~2005) -- bbox stays a superset, correctness unaffected. NEW LOW, documented only: _place_same_crs (__init__.py ~2660) direct-copies a tile whose resolution differs from the output by up to 1%, so a wide tile drifts up to 0.01*width pixels; merge only. TEST-COVERAGE GAP for the test-coverage sweep: general_checks.py has no helper that exercises the projection param extractors, so every fast-path acceptance bug (#2651, #3275, #3697) reached main with green tests; a pyproj-differential test over a CRS list would have caught all three. gdal available but not used -- pyproj is the authoritative reference for the projection kernels and the resampling kernels were checked by invariant instead. cuda-available: cupy and dask+cupy paths executed." +reproject,2026-07-27,3697,HIGH,4;6,"2026-07-27 sweep: filed #3697 (HIGH, cat 4;6) -- numba fast path accepted CRS definitions it does not model; 220 of the 2685 fast-path-accepted EPSG projected CRSs were >5 m from pyproj. Four causes, all fixed in PR for #3697: (a) _lcc_params/_aea_params defaulted lat_2 to an already-radians lat_1 so math.radians ran twice, breaking every 1SP LCC -- 185 codes, worst 3819 m (EPSG:8325); (b) _aea/_cea/_laea/_stere/_sinu/_sterea_params ignored units, reading feet as metres (EPSG:10481, 6594 km); (c) +axis= ignored, so wsu grids came out negated+swapped (EPSG:2046-2055, 7085 km); (d) _is_wgs84_compatible_ellipsoid read a missing +datum= key as WGS84, letting GGRS87/Israel 1993/VN-2000 skip their shift (24 codes, 325 m). After the fix the same EPSG-wide sweep reports 0 CRSs >5 m and the fast path still covers 2649. PRIOR FIXES RE-VERIFIED STILL HOLDING: #3274 (LAEA inverse rq, authalic series), #3275 (sphere/ellipsoid guard), #3276 (itrf ppm). VERIFIED CLEAN this run: merc/emerc/UTM/tmerc(incl. Zb lat_0!=0)/LCC 2SP/AEA/CEA/sinu/LAEA(obliq,equit,N/S pole)/polar stere(lat_ts, k_0, plain)/sterea forward+inverse vs pyproj <=1e-5 m, except the documented ~1.6 mm authalic truncation (cea/laea inverse) and the deliberate NAD83-as-WGS84 approximation (0.3-1.8 m on EPSG:2260/3978/5070, documented at _utm_params). Resampling kernels clean on all 4 backends: constant-raster partition of unity (<=4.4e-15), same-CRS nearest identity bitwise exact, no NaN leak into fully-valid stencils, numpy vs cupy vs dask+numpy vs dask+cupy agree to <=1.2e-13 for nearest/bilinear/cubic. KNOWN LOW, documented only (not fixed): _source_footprint_in_target probe array typo uses x-midpoint mx as a latitude in the last 3 ys entries (__init__.py ~2005) -- bbox stays a superset, correctness unaffected. NEW LOW, documented only: _place_same_crs (__init__.py ~2660) direct-copies a tile whose resolution differs from the output by up to 1% with no resampling, so placement drifts linearly across the tile; verified a 0.5% mismatch on a 1000-px tile puts the last column 5 output pixels (556 m) from its true position. merge() only, and only for inputs whose resolutions disagree. TEST-COVERAGE GAP for the test-coverage sweep: general_checks.py has no helper that exercises the projection param extractors, so every fast-path acceptance bug (#2651, #3275, #3697) reached main with green tests; a pyproj-differential test over a CRS list would have caught all three. gdal available but not used -- pyproj is the authoritative reference for the projection kernels and the resampling kernels were checked by invariant instead. cuda-available: cupy and dask+cupy 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-07-03,3626,MEDIUM,4,"Re-audit post-#1407 fix (ground distance OK). MEDIUM #3626: ray azimuths uniform in cell-INDEX angle not ground azimuth, so anisotropic cellsizes bias the azimuthal quadrature (45-deg ramp, true SVF 0.75: 0.728/0.773/0.693 at 1:1/1:2/2:1 aspect; +0.105 pure weighting error at 4:1); fixed by aiming rays in ground space (cos(phi)/cx, sin(phi)/cy normalized), square-cell results bit-identical. Cats 1-3,5 clean: float64 throughout, NaN center->NaN out, flat=1.0 exactly all 4 backends, flip/rot90/translate invariants at machine eps, numpy==dask==cupy==dask+cupy (max 2.2e-16) incl NaN holes straddling chunk boundaries, map_overlap depth==max_radius==stencil radius, depth>chunksize works (dask 2025.7). Cat6 skipped: richdem-unavailable rvt-unavailable grass-unavailable. LOWs (not fixed): NaN neighbor silently truncates ray (horizon beyond nodata ignored, undocumented but keeps numpy/dask edge semantics consistent); nearest-cell ray rounding gives ~-2.9% jitter bias on 45-deg slopes even square cells (inherent to integer-cell marching, reference tools interpolate along ray). For test-coverage sweep: cross-backend tests had no NaN-input and no anisotropic-cellsize cases (added aniso numpy/dask/cupy parity in #3626 PR); for style sweep: unused import _boundary_to_dask (F401)." From 975cc49548345995faf9b2c772fd321ded040924 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Mon, 27 Jul 2026 14:59:55 -0400 Subject: [PATCH 4/4] Skip rather than error when PROJ lacks a test CRS (#3697) EPSG:10481, 9549 and 9207 are recent additions to the EPSG registry and PROJ ships a snapshot of it, so an older PROJ in CI would raise inside CRS construction instead of reaching the assertion. Route construction through a helper that skips in that case. --- .../test_reproject_fast_path_guard_3697.py | 44 ++++++++++++------- 1 file changed, 29 insertions(+), 15 deletions(-) diff --git a/xrspatial/tests/test_reproject_fast_path_guard_3697.py b/xrspatial/tests/test_reproject_fast_path_guard_3697.py index 55892d653..31833d07f 100644 --- a/xrspatial/tests/test_reproject_fast_path_guard_3697.py +++ b/xrspatial/tests/test_reproject_fast_path_guard_3697.py @@ -32,6 +32,20 @@ WGS84 = 'EPSG:4326' +def _crs_or_skip(spec): + """Build a CRS, skipping if this PROJ build's database lacks the code. + + Several CRSs used below (EPSG:10481, 9549, 9207) are recent additions + to the EPSG registry, and PROJ ships a snapshot of that registry. An + older PROJ should skip rather than error out of an unrelated + assertion. + """ + try: + return pyproj.CRS.from_user_input(spec) + except Exception: + pytest.skip(f'{spec} is not in this PROJ database') + + def _lon_field_error(crs_spec, lon_c, lat_c, span=0.5, n=401, out=32): """Max georeferencing error, in metres, of a reprojected raster. @@ -41,8 +55,8 @@ def _lon_field_error(crs_spec, lon_c, lat_c, span=0.5, n=401, out=32): is a smooth ramp sampled far finer than the output, so what is left is georeferencing error rather than resampling error. """ - tgt = pyproj.CRS.from_user_input(crs_spec) - src_crs = pyproj.CRS.from_epsg(4326) + tgt = _crs_or_skip(crs_spec) + src_crs = _crs_or_skip(4326) fwd = pyproj.Transformer.from_crs(src_crs, tgt, always_xy=True) inv = pyproj.Transformer.from_crs(tgt, src_crs, always_xy=True) @@ -82,7 +96,7 @@ class TestLccSingleStandardParallel: def test_cone_constant_is_sin_lat1(self): # +proj=lcc +lat_1=45 +lat_0=45 with no lat_2 at all. - crs = pyproj.CRS.from_user_input( + crs = _crs_or_skip( '+proj=lcc +lat_1=45 +lat_0=45 +lon_0=-100 +k_0=1 ' '+x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs' ) @@ -102,7 +116,7 @@ def test_reproject_matches_pyproj(self, epsg, lon, lat): def test_two_parallel_lcc_still_uses_both(self): # EPSG:2154 carries lat_1 and lat_2, so the 2SP branch must run and # produce a cone constant strictly between the two parallels. - crs = pyproj.CRS.from_epsg(2154) + crs = _crs_or_skip(2154) d = P._crs_to_dict(crs) assert d['lat_1'] != d['lat_2'] n = P._lcc_params(crs)[1] @@ -111,7 +125,7 @@ def test_two_parallel_lcc_still_uses_both(self): assert _lon_field_error('EPSG:2154', 3.0, 46.5) < 1.0 def test_aea_single_standard_parallel(self): - crs = pyproj.CRS.from_user_input( + crs = _crs_or_skip( '+proj=aea +lat_1=45 +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 ' '+datum=WGS84 +units=m +no_defs' ) @@ -170,13 +184,13 @@ class TestAxisOrientation: def test_wsu_declined(self): # EPSG:2048 is the South African Lo19 grid: +axis=wsu. - crs = pyproj.CRS.from_epsg(2048) + crs = _crs_or_skip(2048) assert P._crs_to_dict(crs).get('axis') == 'wsu' assert P._tmerc_params(crs) is None assert not P._is_wgs84_compatible_ellipsoid(crs) def test_enu_still_accepted(self): - crs = pyproj.CRS.from_epsg(32633) + crs = _crs_or_skip(32633) assert P._is_wgs84_compatible_ellipsoid(crs) def test_reproject_south_oriented_grid(self): @@ -200,7 +214,7 @@ class TestUnaliasedDatumGuard: (9207, 100.0), # VN-2000 / TM-3 104-30 ]) def test_shifted_datum_declined(self, epsg, min_offset_m): - crs = pyproj.CRS.from_epsg(epsg) + crs = _crs_or_skip(epsg) d = P._crs_to_dict(crs) # The dict carries no datum name, which is exactly what used to # make the old name test read these as WGS84. @@ -219,7 +233,7 @@ def test_shifted_datum_declined(self, epsg, min_offset_m): 6350, # NAD83(2011) / Conus Albers ]) def test_aligned_datum_keeps_fast_path(self, epsg): - crs = pyproj.CRS.from_epsg(epsg) + crs = _crs_or_skip(epsg) assert P._is_wgs84_compatible_ellipsoid(crs) offset = P._datum_offset_from_wgs84(crs) if offset is not None: @@ -241,7 +255,7 @@ def test_lite_crs_has_no_geodetic_base(self): def test_measurement_failure_fails_closed(self, monkeypatch): # A CRS we could have measured but could not must be rejected, not # waved through -- an unverifiable datum belongs on the pyproj path. - crs = pyproj.CRS.from_epsg(32633) + crs = _crs_or_skip(32633) assert P._is_wgs84_compatible_ellipsoid(crs) def _boom(*args, **kwargs): @@ -261,7 +275,7 @@ def test_axis_order_only_crs_keeps_fast_path(self): # axis guard must not reject them; see the note in # _is_wgs84_compatible_ellipsoid. for epsg in (2193, 3346): - crs = pyproj.CRS.from_epsg(epsg) + crs = _crs_or_skip(epsg) assert [a.direction for a in crs.axis_info] == ['north', 'east'] assert P._crs_to_dict(crs).get('axis') is None assert P._is_wgs84_compatible_ellipsoid(crs) @@ -277,7 +291,7 @@ class TestBackendParity: @pytest.mark.parametrize('epsg', [2100, 2048, 10481]) def test_cpu_dispatcher_declines(self, epsg): - tgt = pyproj.CRS.from_epsg(epsg) + tgt = _crs_or_skip(epsg) src = pyproj.CRS.from_epsg(4326) fwd = pyproj.Transformer.from_crs(src, tgt, always_xy=True) aou = tgt.area_of_use @@ -294,7 +308,7 @@ def test_cuda_dispatcher_declines(self, epsg): pytest.importorskip('cupy') from xrspatial.reproject._projections_cuda import try_cuda_transform - tgt = pyproj.CRS.from_epsg(epsg) + tgt = _crs_or_skip(epsg) src = pyproj.CRS.from_epsg(4326) fwd = pyproj.Transformer.from_crs(src, tgt, always_xy=True) aou = tgt.area_of_use @@ -309,8 +323,8 @@ def test_cuda_dispatcher_declines(self, epsg): ]) def test_dask_matches_numpy(self, epsg, lon, lat): da = pytest.importorskip('dask.array') - tgt = pyproj.CRS.from_epsg(epsg) - src_crs = pyproj.CRS.from_epsg(4326) + tgt = _crs_or_skip(epsg) + src_crs = _crs_or_skip(4326) fwd = pyproj.Transformer.from_crs(src_crs, tgt, always_xy=True) n, span = 201, 0.4