diff --git a/.claude/sweep-accuracy-state.csv b/.claude/sweep-accuracy-state.csv index 9896a4df0..8c169280c 100644 --- a/.claude/sweep-accuracy-state.csv +++ b/.claude/sweep-accuracy-state.csv @@ -29,7 +29,7 @@ polygon_clip,2026-06-10,3186,HIGH,5,"Cat5 backend inconsistency: dask+cupy clip_ polygonize,2026-05-29,2606,HIGH,5,"Cat 5 HIGH: dask connectivity=8 cross-chunk merge filled diagonal notch where same-value regions meet only at a corner across a chunk boundary; total area exceeded raster. Hole ring was dropped because containment tested hole[0] (on exterior at pinch). Fixed via _ring_interior_point in PR for #2606. numpy, dask+numpy, dask+cupy area parity now holds; 4-conn was already correct. cupy + dask+cupy paths validated on GPU host. Other cats clean: NaN masked on numpy/cupy float paths (tested), _is_close handles +/-inf via exact-equality short-circuit, atol/rtol/simplify_tolerance reject NaN/inf, integer GPU CCL matches numpy." proximity,2026-06-09,3108,HIGH,4;5,"Cat5/Cat4: bounded GREAT_CIRCLE dask (numpy+cupy) missed targets across the +/-180 antimeridian seam: _halo_depth sized x-halo as linear parallel-arc sum, but haversine is periodic in lon and chords shorten near poles, so array-space adjacency is no lower bound on spherical distance; numpy/cupy (brute force) found the wrap target (~111 km), dask returned NaN. Fixed in #3108 via chord bound 2R asin(cos(lat_max)|sin(dlon/2)|) + x-axis fold when seam/180-deg chord within max_distance (covers over-pole too). CUDA host: cupy + dask+cupy executed, 417+ tests pass. Cat1-3 clean (float32 output documented; NaN via isfinite consistent; bounds guards correct; tie-break unified in #2881). LOW (not fixed): great_circle_distance uses WGS84 equatorial radius 6378137 as sphere radius (~0.1% vs mean-radius convention) but documented and exposed as param." rasterize,2026-06-09,3085,MEDIUM,2,"Cat2: non-finite burn values (NaN/inf, e.g. GeoDataFrame column with missing data) against integer output dtype silently cast to platform sentinel (NaN -> INT_MIN) on all 4 backends; bool dtype collapses NaN to True. numpy path suppresses the cast RuntimeWarning so it is fully silent. Fix #3085 mirrors NaN-fill guard #2504 and unsafe-int guard #3056; merge='count' exempt (never reads props). Cats 1/3/4/5 clean: CUDA available, ran 4-backend parity probe (mixed polygon/line/point, all 6 builtin merges, all_touched, chunked 17x23) -- bit-identical across numpy/cupy/dask+numpy/dask+cupy; 665 rasterize tests pass. GPU atomic min/max NaN handling (#2255) verified correct; non-atomic _merge_min_gpu/_merge_max_gpu device fns lack the NaN branch but are unreachable dead code (string merges always take the atomic ladder). Scanline ceil (GPU int-trunc emulation) matches CPU np.ceil for negative x. No Earth-curvature surface (planar vector rasterization by design)." -reproject,2026-06-09,3094,HIGH,4;5,"HIGH #3094: try_cuda_transform missing the #2651 non-WGS84 datum guard; cupy/dask+cupy 4326<->27700 (and DHDN/MGI/ED50/NAD27 entries) project with WGS84 Krueger + no datum shift, ~80-100 m coordinate error, end-to-end numpy-vs-cupy value diff 0.094 on [0,1] data with 5 NaN-mask flips (verified on GPU). MEDIUM #3096: empty-chunk returns in _reproject_chunk_numpy/_cupy and _reproject_block_adapter hardcode float64, so integer-source dask reproject computes float64 when any chunk misses the source footprint while meta advertises int (eager backend returns int). LOW (doc only): _place_same_crs accepts 1% res mismatch for direct placement, up to ~1% of tile width misregistration at the far edge on large same-CRS merge tiles. Cats 1-3 clean: f64 kernels, GDAL-style NaN renorm, correct bounds guards; projection series/constants match PROJ. CUDA available; GPU paths executed." +reproject,2026-06-12,3274,HIGH,1;4,"3 confirmed bugs, all kernel-vs-PROJ parity: #3274 HIGH LAEA inverse spurious /rq (2.6 km err for 3035) + _authalic_apa inverse series wrong (4.8 m in AEA/CEA inverses; PROJ 3-term = 1.6 mm), CPU+CUDA kernels both; #3275 HIGH _is_wgs84_compatible_ellipsoid passes R-defined spheres (MODIS sinusoidal 18.9 km err) and _aea_params/_cea_params lack the guard entirely (23.8 km on spherical aea/cea); #3276 MEDIUM itrf helmert scale 1e-9 but PROJ +s is ppm (1e-6), ~23 mm err. Verified clean: merc/emerc/UTM/tmerc/LCC/polar stere (incl lat_ts akm1) forward+inverse <=1e-5 m vs pyproj; resampling kernels NaN handling and GDAL renorm match across numpy/cupy/dask (CUDA run, gpu-vs-cpu 1.3e-7); dask footprint chunk-skip bbox is a superset in all probed cases (no holes). LOW (documented only): _source_footprint_in_target probe array typo uses x-midpoint mx as a latitude in last 3 ys entries (bbox superset, correctness unaffected)." resample,2026-05-29,2610,HIGH,3;5,"dask interp (nearest/bilinear) overlap depth=1 too small on downsample; block-centered source coord landed past chunk, map_coordinates clamped to edge -> wrong seam rows. Fixed PR #2627 via per-axis _downsample_radius. cupy+dask+cupy verified." sieve,2026-04-13T12:00:00Z,,,,Union-find CCL correct. NaN excluded from labeling. All backends funnel through _sieve_numpy. sky_view_factor,2026-05-01,1407,HIGH,4,Horizon angle ignored cell size; fixed by passing cellsize_x/cellsize_y into CPU+GPU kernels and using ground distance diff --git a/xrspatial/reproject/_projections.py b/xrspatial/reproject/_projections.py index 137be50eb..5fba1cca9 100644 --- a/xrspatial/reproject/_projections.py +++ b/xrspatial/reproject/_projections.py @@ -229,25 +229,29 @@ def _authalic_q(sinphi, e): def _authalic_apa(e): - """Precompute 6 coefficients for the authalic latitude inverse series. - - Returns array [APA0..APA5] used by _authalic_inv. - 6 terms give sub-centimetre accuracy (vs ~4m with 3 terms). - Coefficients from Snyder (1987) / Karney (2011). + """Precompute coefficients for the authalic latitude inverse series. + + Returns array [APA0..APA5] used by _authalic_inv (terms 4-6 are + zero, kept so existing kernel signatures stay unchanged). + + Coefficients match PROJ's pj_authset (Snyder 1987, eq. 3-18): + phi = beta + APA0*sin(2*beta) + APA1*sin(4*beta) + APA2*sin(6*beta). + Measured against exact numerical inversion of _authalic_q this + series is accurate to ~2.5e-10 rad (1.6 mm) for the WGS84 + eccentricity. The previous coefficients (e.g. 17/360 instead of + 23/360 at the e^4 order of APA1) did not invert _authalic_q and + were off by up to 7.5e-7 rad (~4.8 m); see GH #3274. """ e2 = e * e e4 = e2 * e2 e6 = e4 * e2 - e8 = e6 * e2 - e10 = e8 * e2 apa = np.empty(6, dtype=np.float64) - apa[0] = (e2 / 3.0 + 31.0 * e4 / 180.0 + 59.0 * e6 / 560.0 - + 17141.0 * e8 / 166320.0 + 28289.0 * e10 / 249480.0) - apa[1] = 17.0 * e4 / 360.0 + 61.0 * e6 / 1260.0 + 10217.0 * e8 / 120960.0 + 319.0 * e10 / 3024.0 - apa[2] = 383.0 * e6 / 45360.0 + 34729.0 * e8 / 1814400.0 + 192757.0 * e10 / 5765760.0 - apa[3] = 6007.0 * e8 / 272160.0 + 36941.0 * e10 / 1270080.0 - apa[4] = 33661.0 * e10 / 5765760.0 - apa[5] = 0.0 # 12th order term negligible for Earth + apa[0] = e2 / 3.0 + 31.0 * e4 / 180.0 + 517.0 * e6 / 5040.0 + apa[1] = 23.0 * e4 / 360.0 + 251.0 * e6 / 3780.0 + apa[2] = 761.0 * e6 / 45360.0 + apa[3] = 0.0 + apa[4] = 0.0 + apa[5] = 0.0 return apa @@ -255,7 +259,7 @@ def _authalic_apa(e): def _authalic_inv(beta, apa): """Inverse authalic latitude: beta (authalic, rad) -> phi (geodetic, rad). - 6-term Fourier series for sub-centimetre accuracy. + Snyder 3-term series (PROJ pj_authlat); ~1.6 mm max error on WGS84. """ t = 2.0 * beta return (beta @@ -871,13 +875,23 @@ def _laea_inv_point(x, y, lon0, sinb1, cosb1, else: lam = math.atan2(x_a, -y_a) else: # OBLIQ or EQUIT - # PROJ: x /= dd, y *= dd (undo the xmf/ymf scaling) - xn = x / (a * xmf) # = x / (a * rq * dd) - yn = y / (a * ymf) # = y / (a * rq / dd) = y * dd / (a * rq) + # PROJ: x /= dd, y *= dd (undo only the dd part of the xmf/ymf + # scaling). The rq factor must stay in rho because the angular + # distance below divides by rq exactly once. Dividing by + # (a * xmf) alone strips rq here and then divides by it again + # in sce, which inflated distances by ~0.11% (#3274). + xn = x / (a * xmf) * rq # = x / (a * dd) + yn = y / (a * ymf) * rq # = y * dd / a rho = math.hypot(xn, yn) if rho < 1e-30: return math.degrees(lon0), math.degrees(math.asin(sinb1)) - sce = 2.0 * math.asin(0.5 * rho / rq) + # Clamp the asin argument: points beyond the projection disc + # (rho > 2*rq) are invalid input; clamping matches the ratio + # clamps elsewhere in this kernel instead of emitting NaN. + half_rho_rq = 0.5 * rho / rq + if half_rho_rq > 1.0: + half_rho_rq = 1.0 + sce = 2.0 * math.asin(half_rho_rq) sinz = math.sin(sce) cosz = math.cos(sce) if mode == 0: # OBLIQ diff --git a/xrspatial/reproject/_projections_cuda.py b/xrspatial/reproject/_projections_cuda.py index 46d30af16..98f6183fa 100644 --- a/xrspatial/reproject/_projections_cuda.py +++ b/xrspatial/reproject/_projections_cuda.py @@ -584,12 +584,20 @@ def _d_laea_inv(x, y, lon0, sinb1, cosb1, else: lam = math.atan2(x_a, -y_a) else: - xn = x / (a * xmf) - yn = y / (a * ymf) + # Keep the rq factor in rho: sce below divides by rq exactly + # once (mirrors the CPU kernel fix for GH #3274). + xn = x / (a * xmf) * rq # = x / (a * dd) + yn = y / (a * ymf) * rq # = y * dd / a rho = math.hypot(xn, yn) if rho < 1e-30: return math.degrees(lon0), math.degrees(math.asin(sinb1)) - sce = 2.0 * math.asin(0.5 * rho / rq) + # Clamp the asin argument: points beyond the projection disc + # (rho > 2*rq) are invalid input; clamping matches the ratio + # clamps elsewhere in this kernel instead of emitting NaN. + half_rho_rq = 0.5 * rho / rq + if half_rho_rq > 1.0: + half_rho_rq = 1.0 + sce = 2.0 * math.asin(half_rho_rq) sinz = math.sin(sce) cosz = math.cos(sce) if mode == 0: diff --git a/xrspatial/tests/test_reproject_inverse_kernels_3274.py b/xrspatial/tests/test_reproject_inverse_kernels_3274.py new file mode 100644 index 000000000..5cafbac15 --- /dev/null +++ b/xrspatial/tests/test_reproject_inverse_kernels_3274.py @@ -0,0 +1,256 @@ +"""Regression tests for GH #3274. + +Two inverse-kernel bugs in ``xrspatial/reproject/_projections.py``: + +1. ``_laea_inv_point`` (oblique/equatorial) normalized the input by + ``x / (a * xmf)``, stripping the ``rq`` factor that the angular + distance ``sce = 2 * asin(0.5 * rho / rq)`` divides out. The double + division inflated distances by ~0.11% of the distance from the + projection origin (up to 2.6 km for EPSG:3035 over Europe). +2. ``_authalic_apa`` carried inverse-series coefficients that did not + invert ``_authalic_q`` (wrong at leading order in the second term: + 17/360 vs Snyder/PROJ's 23/360), putting every authalic-latitude + inverse (AEA, CEA, LAEA) off by up to ~4.8 m. + +The CUDA kernels in ``_projections_cuda.py`` shared both bugs. +""" +from __future__ import annotations + +import math + +import numpy as np +import pytest + +from xrspatial.utils import has_cuda_and_cupy + +try: + import pyproj + HAS_PYPROJ = True +except ImportError: + HAS_PYPROJ = False + +HAS_CUPY = has_cuda_and_cupy() + +pytestmark = pytest.mark.skipif( + not HAS_PYPROJ, reason="pyproj required for reproject tests" +) + + +def _inverse_parity_max_err_m(crs_code_or_obj, lons, lats): + """Max error (metres) of the numba inverse fast path vs pyproj. + + Forward-projects (lons, lats) with pyproj, runs the projected points + through ``transform_points`` back to EPSG:4326, and compares against + pyproj's own inverse. + """ + from xrspatial.reproject._projections import transform_points + + crs = pyproj.CRS(crs_code_or_obj) + geo = pyproj.CRS(4326) + fwd = pyproj.Transformer.from_crs(geo, crs, always_xy=True) + px, py = fwd.transform(lons, lats) + px = np.asarray(px) + py = np.asarray(py) + mask = np.isfinite(px) & np.isfinite(py) + px, py = px[mask], py[mask] + + result = transform_points(crs, geo, px, py) + assert result is not None, "expected a numba fast path for this pair" + tx, ty = result + + inv = pyproj.Transformer.from_crs(crs, geo, always_xy=True) + rx, ry = inv.transform(px, py) + rx = np.asarray(rx) + ry = np.asarray(ry) + err_deg_x = (tx - rx) * np.cos(np.radians(ry)) + err_deg_y = ty - ry + return float(np.hypot(err_deg_x, err_deg_y).max() * 111320.0) + + +class TestAuthalicInverseSeries: + def test_round_trips_authalic_q(self): + """phi -> q -> beta -> series must return phi to ~1e-9 rad. + + The old coefficients were off by up to 7.5e-7 rad (~4.8 m). + """ + from xrspatial.reproject._projections import (_authalic_apa, _authalic_q, _WGS84_E) + + apa = _authalic_apa(_WGS84_E) + qp = _authalic_q(1.0, _WGS84_E) + worst = 0.0 + for phi_deg in np.linspace(-89.5, 89.5, 359): + phi = math.radians(phi_deg) + beta = math.asin(_authalic_q(math.sin(phi), _WGS84_E) / qp) + t = 2.0 * beta + phi_back = beta + sum( + c * math.sin((k + 1) * t) for k, c in enumerate(apa[:5]) + ) + worst = max(worst, abs(phi_back - phi)) + # 5e-9 rad is ~3 cm on the ground; the correct series sits at + # ~2.5e-10 rad, the broken one at 7.5e-7 rad. + assert worst < 5e-9 + + def test_matches_proj_pj_authset_coefficients(self): + from xrspatial.reproject._projections import _authalic_apa, _WGS84_E2, _WGS84_E + + e2 = _WGS84_E2 + e4 = e2 * e2 + e6 = e4 * e2 + apa = _authalic_apa(_WGS84_E) + assert apa[0] == pytest.approx(e2 / 3 + 31 * e4 / 180 + 517 * e6 / 5040, rel=1e-15) + assert apa[1] == pytest.approx(23 * e4 / 360 + 251 * e6 / 3780, rel=1e-15) + assert apa[2] == pytest.approx(761 * e6 / 45360, rel=1e-15) + + +class TestInverseKernelParityVsPyproj: + """Inverse fast paths must match pyproj at the cm level. + + Tolerances are set an order of magnitude above the measured parity + (~1.6 mm) and three orders below the pre-fix errors (4.8 m / 2.6 km). + """ + + def test_laea_3035_inverse(self): + rng = np.random.default_rng(3274) + lons = rng.uniform(-10, 30, 500) + lats = rng.uniform(35, 70, 500) + assert _inverse_parity_max_err_m(3035, lons, lats) < 0.05 + + def test_cea_6933_inverse(self): + rng = np.random.default_rng(3274) + lons = rng.uniform(-179, 179, 500) + lats = rng.uniform(-85, 85, 500) + assert _inverse_parity_max_err_m(6933, lons, lats) < 0.05 + + def test_aea_wgs84_inverse(self): + # EPSG:5070 is NAD83; pyproj inserts a ~1-2 m WGS84->NAD83 step + # that the fast path deliberately skips, so the regression test + # uses an AEA definition on WGS84 directly. + aea = pyproj.CRS( + "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 " + "+x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs" + ) + rng = np.random.default_rng(3274) + lons = rng.uniform(-120, -70, 500) + lats = rng.uniform(25, 49, 500) + assert _inverse_parity_max_err_m(aea, lons, lats) < 0.05 + + def test_laea_polar_inverse_unchanged(self): + # The polar LAEA branch never had the rq bug; pin its parity so + # the fix does not disturb it. + laea_n = pyproj.CRS( + "+proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 " + "+ellps=WGS84 +units=m +no_defs" + ) + rng = np.random.default_rng(3274) + lons = rng.uniform(-180, 180, 500) + lats = rng.uniform(55, 89, 500) + assert _inverse_parity_max_err_m(laea_n, lons, lats) < 0.05 + + def test_laea_3035_forward_unchanged(self): + from xrspatial.reproject._projections import transform_points + + rng = np.random.default_rng(3274) + lons = rng.uniform(-10, 30, 500) + lats = rng.uniform(35, 70, 500) + geo = pyproj.CRS(4326) + crs = pyproj.CRS(3035) + tx, ty = transform_points(geo, crs, lons, lats) + fwd = pyproj.Transformer.from_crs(geo, crs, always_xy=True) + rx, ry = fwd.transform(lons, lats) + err = np.hypot(tx - np.asarray(rx), ty - np.asarray(ry)) + assert float(err.max()) < 1e-3 + + +class TestLaeaReprojectEndToEnd: + def test_fast_path_matches_exact_pyproj_path(self): + """reproject(4326 -> 3035) fast path vs transform_precision=0. + + transform_precision=0 is the documented exact-pyproj escape + hatch; before the fix the fast path sampled source pixels up to + 2.6 km away from where the exact path sampled them. + """ + import xarray as xr + + from xrspatial.reproject import reproject + + ny, nx = 200, 250 + y = np.linspace(70, 35, ny) + x = np.linspace(-10, 30, nx) + xx, yy = np.meshgrid(x, y) + data = np.sin(np.radians(xx * 3)) * np.cos(np.radians(yy * 3)) * 1000 + raster = xr.DataArray( + data, dims=['y', 'x'], coords={'y': y, 'x': x}, + attrs={'crs': 'EPSG:4326'}, + ) + + exact = reproject(raster, 'EPSG:3035', transform_precision=0) + fast = reproject(raster, 'EPSG:3035') + mask = np.isfinite(exact.values) & np.isfinite(fast.values) + assert mask.sum() > 0.5 * mask.size + # The remaining difference is the control-grid interpolation of + # the approximate transform, well under a thousandth of the data + # range. Pre-fix this was ~1.6 of a +/-1000 range at 7 km pixels. + diff = np.abs(exact.values[mask] - fast.values[mask]) + assert float(diff.max()) < 0.05 + + +@pytest.mark.skipif(not HAS_CUPY, reason="CUDA/cupy not available") +class TestCudaLaeaInverseParity: + def test_cuda_laea_inverse_matches_pyproj(self): + """try_cuda_transform for a 3035 chunk must match pyproj per pixel.""" + import cupy as cp + + from xrspatial.reproject._projections_cuda import try_cuda_transform + + geo = pyproj.CRS(4326) + crs = pyproj.CRS(3035) + # A chunk over central Europe, far enough from the projection + # origin that the old bug showed km-scale errors. + fwd = pyproj.Transformer.from_crs(geo, crs, always_xy=True) + x0, y0 = fwd.transform(25.0, 40.0) + h = w = 16 + cb = (x0 - 80000, y0 - 80000, x0 + 80000, y0 + 80000) + result = try_cuda_transform(geo, crs, cb, (h, w)) + assert result is not None + src_y, src_x = (cp.asnumpy(a) for a in result) + + res_x = (cb[2] - cb[0]) / w + res_y = (cb[3] - cb[1]) / h + xs = cb[0] + (np.arange(w) + 0.5) * res_x + ys = cb[3] - (np.arange(h) + 0.5) * res_y + xx, yy = np.meshgrid(xs, ys) + inv = pyproj.Transformer.from_crs(crs, geo, always_xy=True) + rlon, rlat = inv.transform(xx, yy) + err_m = np.hypot( + (src_x - rlon) * np.cos(np.radians(rlat)), src_y - rlat + ) * 111320.0 + assert float(err_m.max()) < 0.05 + + def test_cuda_matches_cpu_fast_path(self): + import cupy as cp + + from xrspatial.reproject._projections import try_numba_transform + from xrspatial.reproject._projections_cuda import try_cuda_transform + + geo = pyproj.CRS(4326) + crs = pyproj.CRS(3035) + cb = (4000000.0, 2000000.0, 4500000.0, 2500000.0) + shape = (32, 32) + cpu = try_numba_transform(geo, crs, cb, shape) + gpu = try_cuda_transform(geo, crs, cb, shape) + assert cpu is not None and gpu is not None + np.testing.assert_allclose(cp.asnumpy(gpu[0]), cpu[0], atol=1e-9) + np.testing.assert_allclose(cp.asnumpy(gpu[1]), cpu[1], atol=1e-9) + + +class TestLaeaEquatorialInverse: + def test_laea_equit_inverse(self): + """EQUIT-mode LAEA (lat_0=0) shares the fixed OBLIQ code path.""" + laea_eq = pyproj.CRS( + "+proj=laea +lat_0=0 +lon_0=20 +x_0=0 +y_0=0 " + "+ellps=WGS84 +units=m +no_defs" + ) + rng = np.random.default_rng(3274) + lons = rng.uniform(-20, 60, 500) + lats = rng.uniform(-40, 40, 500) + assert _inverse_parity_max_err_m(laea_eq, lons, lats) < 0.05