diff --git a/.claude/sweep-accuracy-state.csv b/.claude/sweep-accuracy-state.csv index 6d2b39673..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-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% 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)." diff --git a/xrspatial/reproject/_projections.py b/xrspatial/reproject/_projections.py index 44648cf39..3eacadafb 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,76 @@ 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. + + 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: + 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 math.inf + _, _, dist = pyproj.Geod(ellps='WGS84').inv(lon, lat, lon_w, lat_w) + return abs(float(dist)) + except Exception: + return math.inf + + +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 +409,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 +555,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 +667,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 +790,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 +858,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 +1054,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 +1196,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)) @@ -1773,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 @@ -1792,6 +1897,21 @@ def _is_wgs84_compatible_ellipsoid(crs): d = _crs_to_dict(crs) except Exception: return False + # 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 # Explicit sphere radius: never WGS84-compatible. if 'R' in d: return False @@ -1809,9 +1929,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..31833d07f --- /dev/null +++ b/xrspatial/tests/test_reproject_fast_path_guard_3697.py @@ -0,0 +1,353 @@ +"""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 _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. + + 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 = _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) + + 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 = _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' + ) + 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 = _crs_or_skip(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 = _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' + ) + 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 = _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 = _crs_or_skip(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 = _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. + 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 = _crs_or_skip(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 + + 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 = _crs_or_skip(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 = _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) + assert P._tmerc_params(crs) is not 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 = _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 + 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 = _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 + 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 = _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 + 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, + )