From 66ffcd7ab414848518618e45c003a4226c451671 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Fri, 12 Jun 2026 07:51:50 -0700 Subject: [PATCH 1/3] Record accuracy sweep state for reproject (#3274, #3275, #3276) --- .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 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 From dfa0fd835b5074ae42c0c538ff52d95c41a24bcc Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Fri, 12 Jun 2026 07:55:46 -0700 Subject: [PATCH 2/3] Fix LAEA inverse rq normalization and authalic inverse series (#3274) The oblique/equatorial LAEA inverse divided out the rq factor twice (once in the input normalization, again in the angular distance), inflating distances by ~0.11% -- up to 2.6 km for EPSG:3035 over Europe. Keep rq in rho so sce divides by it exactly once, matching PROJ's laea.cpp. Same fix in the CUDA kernel. The authalic latitude inverse series coefficients did not invert _authalic_q (4.8 m max error). Replace them with PROJ's pj_authset / Snyder eq. 3-18 coefficients (1.6 mm max error), which feed the AEA, CEA, and LAEA inverses on both CPU and CUDA paths. Measured parity vs pyproj after the fix: 3035 inverse 1.0 mm (was 2.6 km), 6933 inverse 1.6 mm (was 4.8 m), WGS84 AEA inverse 1.5 mm (was ~4.9 m). Forward kernels unchanged. --- xrspatial/reproject/_projections.py | 44 ++-- xrspatial/reproject/_projections_cuda.py | 6 +- .../test_reproject_inverse_kernels_3274.py | 243 ++++++++++++++++++ 3 files changed, 273 insertions(+), 20 deletions(-) create mode 100644 xrspatial/tests/test_reproject_inverse_kernels_3274.py diff --git a/xrspatial/reproject/_projections.py b/xrspatial/reproject/_projections.py index 137be50eb..e8fb439a4 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,9 +875,13 @@ 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)) diff --git a/xrspatial/reproject/_projections_cuda.py b/xrspatial/reproject/_projections_cuda.py index 46d30af16..8175b2a06 100644 --- a/xrspatial/reproject/_projections_cuda.py +++ b/xrspatial/reproject/_projections_cuda.py @@ -584,8 +584,10 @@ 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)) 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..b256edd9d --- /dev/null +++ b/xrspatial/tests/test_reproject_inverse_kernels_3274.py @@ -0,0 +1,243 @@ +"""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) From cc2b42be8be1cea511c0b75e746a2f05bc8b526a Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Fri, 12 Jun 2026 08:05:29 -0700 Subject: [PATCH 3/3] Address review: clamp LAEA inverse asin argument, pin EQUIT mode (#3274) --- xrspatial/reproject/_projections.py | 8 +++++++- xrspatial/reproject/_projections_cuda.py | 8 +++++++- .../tests/test_reproject_inverse_kernels_3274.py | 13 +++++++++++++ 3 files changed, 27 insertions(+), 2 deletions(-) diff --git a/xrspatial/reproject/_projections.py b/xrspatial/reproject/_projections.py index e8fb439a4..5fba1cca9 100644 --- a/xrspatial/reproject/_projections.py +++ b/xrspatial/reproject/_projections.py @@ -885,7 +885,13 @@ def _laea_inv_point(x, y, lon0, sinb1, cosb1, 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 8175b2a06..98f6183fa 100644 --- a/xrspatial/reproject/_projections_cuda.py +++ b/xrspatial/reproject/_projections_cuda.py @@ -591,7 +591,13 @@ def _d_laea_inv(x, y, lon0, sinb1, cosb1, 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 index b256edd9d..5cafbac15 100644 --- a/xrspatial/tests/test_reproject_inverse_kernels_3274.py +++ b/xrspatial/tests/test_reproject_inverse_kernels_3274.py @@ -241,3 +241,16 @@ def test_cuda_matches_cpu_fast_path(self): 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