Fix LAEA inverse rq normalization and authalic inverse series

.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

xrspatial/reproject/_projections.py

@@ -229,33 +229,37 @@ 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
 
 
 @njit(nogil=True, cache=True)
 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

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:

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