From 9a79981b2e40c27e16bf1a6d8839f243f7218adb Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Mon, 30 Mar 2026 11:22:26 -0700 Subject: [PATCH 1/2] Add longitude normalization to CUDA inverse projection kernels (#1088) Six CUDA inverse projection kernels returned raw `lam + lon0` without wrapping to [-pi, pi]. The CPU Numba kernels all pass through `_norm_lon_rad()` before converting to degrees. Without normalization, coordinates far from the central meridian (e.g., near the antimeridian) can produce longitude values outside [-180, 180]. Added `_d_norm_lon_rad` device function and applied it to: - LCC inverse (2 return paths) - AEA inverse - CEA inverse - Sinusoidal inverse - LAEA inverse - Polar Stereographic inverse Also removed dead code in `_tmerc_params()` where a loop computed a value that was immediately overwritten. --- xrspatial/reproject/_projections.py | 13 ++---- xrspatial/reproject/_projections_cuda.py | 23 +++++++--- xrspatial/tests/test_reproject.py | 54 ++++++++++++++++++++++++ 3 files changed, 74 insertions(+), 16 deletions(-) diff --git a/xrspatial/reproject/_projections.py b/xrspatial/reproject/_projections.py index 0b3739b6..2cd9883a 100644 --- a/xrspatial/reproject/_projections.py +++ b/xrspatial/reproject/_projections.py @@ -1739,17 +1739,12 @@ def _tmerc_params(crs): else: # Conformal latitude of origin Z = lat_0 + _clenshaw_sin_py(_CBG, lat_0) - # Forward Krueger correction at Ce=0 (central meridian) + # Forward Krueger correction at Ce=0 (central meridian): + # sin2=sin(2Z), cos2=cos(2Z), sinh2=0, cosh2=1 sin2Z = math.sin(2.0 * Z) cos2Z = math.cos(2.0 * Z) - dCn = 0.0 - for k in range(5, -1, -1): - dCn = cos2Z * dCn + _ALPHA[k] * sin2Z - # This is a simplified Clenshaw for Ce=0 (sinh=0, cosh=1) - # Actually, use the proper complex Clenshaw with Ce=0: - # sin2=sin(2Z), cos2=cos(2Z), sinh2=0, cosh2=1 - dCn_val = _clenshaw_complex_py(_ALPHA, sin2Z, cos2Z, 0.0, 1.0) - Zb = -Qn * (Z + dCn_val) + dCn = _clenshaw_complex_py(_ALPHA, sin2Z, cos2Z, 0.0, 1.0) + Zb = -Qn * (Z + dCn) return lon_0, k0, fe, fn, Zb, to_meter diff --git a/xrspatial/reproject/_projections_cuda.py b/xrspatial/reproject/_projections_cuda.py index 7dc95c8b..5c79223b 100644 --- a/xrspatial/reproject/_projections_cuda.py +++ b/xrspatial/reproject/_projections_cuda.py @@ -90,6 +90,15 @@ def _d_clenshaw_complex(a0, a1, a2, a3, a4, a5, dCe = sin2Cn * cosh2Ce * hi + cos2Cn * sinh2Ce * hr return dCn, dCe + @cuda.jit(device=True) + def _d_norm_lon_rad(lon): + """Normalize longitude to [-pi, pi].""" + while lon > math.pi: + lon -= 2.0 * math.pi + while lon < -math.pi: + lon += 2.0 * math.pi + return lon + # ----------------------------------------------------------------- # Web Mercator (EPSG:3857) -- spherical # ----------------------------------------------------------------- @@ -291,11 +300,11 @@ def _d_lcc_inv(x, y, lon0, n, c, rho0, k0, e, a): lam_n = math.atan2(x, rho0_y) if abs(rho) < 1e-30: lat = 90.0 if n > 0 else -90.0 - return math.degrees(lon0 + lam_n / n), lat + return math.degrees(_d_norm_lon_rad(lon0 + lam_n / n)), lat ts = math.pow(rho / (a * k0 * c), 1.0 / n) taup = math.sinh(math.log(1.0 / ts)) tau = _d_pj_sinhpsi2tanphi(taup, e) - return math.degrees(lam_n / n + lon0), math.degrees(math.atan(tau)) + return math.degrees(_d_norm_lon_rad(lam_n / n + lon0)), math.degrees(math.atan(tau)) @cuda.jit def _k_lcc_inverse(out_src_x, out_src_y, @@ -355,7 +364,7 @@ def _d_aea_inv(x, y, lon0, n, C, rho0, e, a, qp, ratio = -1.0 beta = math.asin(ratio) phi = _d_authalic_inv(beta, apa0, apa1, apa2, apa3, apa4) - return math.degrees(theta / n + lon0), math.degrees(phi) + return math.degrees(_d_norm_lon_rad(theta / n + lon0)), math.degrees(phi) @cuda.jit def _k_aea_inverse(out_src_x, out_src_y, @@ -404,7 +413,7 @@ def _d_cea_inv(x, y, lon0, k0, e, a, qp, apa0, apa1, apa2, apa3, apa4): ratio = -1.0 beta = math.asin(ratio) phi = _d_authalic_inv(beta, apa0, apa1, apa2, apa3, apa4) - return math.degrees(lam + lon0), math.degrees(phi) + return math.degrees(_d_norm_lon_rad(lam + lon0)), math.degrees(phi) @cuda.jit def _k_cea_inverse(out_src_x, out_src_y, @@ -476,7 +485,7 @@ def _d_sinu_inv(x, y, lon0, e2, a, en0, en1, en2, en3, en4): lam = 0.0 else: lam = x * math.sqrt(1.0 - e2 * s * s) / (a * c) - return math.degrees(lam + lon0), math.degrees(phi) + return math.degrees(_d_norm_lon_rad(lam + lon0)), math.degrees(phi) @cuda.jit def _k_sinu_inverse(out_src_x, out_src_y, @@ -594,7 +603,7 @@ def _d_laea_inv(x, y, lon0, sinb1, cosb1, ratio = -1.0 beta = math.asin(ratio) phi = _d_authalic_inv(beta, apa0, apa1, apa2, apa3, apa4) - return math.degrees(lam + lon0), math.degrees(phi) + return math.degrees(_d_norm_lon_rad(lam + lon0)), math.degrees(phi) @cuda.jit def _k_laea_inverse(out_src_x, out_src_y, @@ -670,7 +679,7 @@ def _d_stere_inv(x, y, lon0, akm1, e, is_south): phi = phi_new if is_south: phi = -phi - return math.degrees(lam + lon0), math.degrees(phi) + return math.degrees(_d_norm_lon_rad(lam + lon0)), math.degrees(phi) @cuda.jit def _k_stere_inverse(out_src_x, out_src_y, diff --git a/xrspatial/tests/test_reproject.py b/xrspatial/tests/test_reproject.py index d883d9ef..fd86351a 100644 --- a/xrspatial/tests/test_reproject.py +++ b/xrspatial/tests/test_reproject.py @@ -1315,6 +1315,60 @@ def test_bounds_overlap(self): assert not _bounds_overlap(a, (0, 11, 10, 20)) # no overlap y +class TestLongitudeNormalization: + """CPU projection round-trips should keep longitude in [-180, 180] (#1088).""" + + def test_sinusoidal_round_trip_stays_in_range(self): + """Sinusoidal inverse must normalize longitude near antimeridian.""" + from xrspatial.reproject._projections import ( + _sinu_fwd_point, _sinu_inv_point, _MLFN_EN, + ) + # Forward: WGS84 point near antimeridian + lon_in, lat_in = 179.5, 30.0 + lon0 = 0.0 # central meridian at 0 + x, y = _sinu_fwd_point(lon_in, lat_in, lon0, _WGS84_E2, _WGS84_A, _MLFN_EN) + # Inverse: should return longitude in [-180, 180] + lon_out, lat_out = _sinu_inv_point(x, y, lon0, _WGS84_E2, _WGS84_A, _MLFN_EN) + assert -180 <= lon_out <= 180, f"lon {lon_out} outside [-180, 180]" + assert abs(lon_out - lon_in) < 1e-6 + assert abs(lat_out - lat_in) < 1e-6 + + def test_lcc_round_trip_stays_in_range(self): + """LCC inverse must normalize longitude.""" + from xrspatial.reproject._projections import ( + _lcc_fwd_point, _lcc_inv_point, _WGS84_E, _WGS84_A, + ) + import math + # EPSG:2154 (France): lon0=3, lat1=44, lat2=49 + lon0 = math.radians(3.0) + lat1, lat2, lat0 = math.radians(44.0), math.radians(49.0), math.radians(46.5) + e = _WGS84_E + a = _WGS84_A + k0 = 1.0 + # Compute n, c, rho0 for LCC + from xrspatial.reproject._projections import _pj_tsfn + s1, s2 = math.sin(lat1), math.sin(lat2) + ts1 = _pj_tsfn(lat1, s1, e) + ts2 = _pj_tsfn(lat2, s2, e) + m1 = math.cos(lat1) / math.sqrt(1.0 - e * e * s1 * s1) + m2 = math.cos(lat2) / math.sqrt(1.0 - e * e * s2 * s2) + n = (math.log(m1) - math.log(m2)) / (math.log(ts1) - math.log(ts2)) + c = m1 / (n * math.pow(ts1, n)) + ts0 = _pj_tsfn(lat0, math.sin(lat0), e) + rho0 = a * k0 * c * math.pow(ts0, n) + # Forward + inverse round trip + lon_in, lat_in = 2.5, 47.0 + x, y = _lcc_fwd_point(lon_in, lat_in, lon0, n, c, rho0, k0, e, a) + lon_out, lat_out = _lcc_inv_point(x, y, lon0, n, c, rho0, k0, e, a) + assert -180 <= lon_out <= 180 + assert abs(lon_out - lon_in) < 1e-6 + assert abs(lat_out - lat_in) < 1e-6 + + +_WGS84_E2 = 2.0 * (1.0 / 298.257223563) - (1.0 / 298.257223563) ** 2 +_WGS84_A = 6378137.0 + + class TestReprojWithLiteCRS: def test_reproject_wgs84_to_utm_with_lite_crs(self): import xarray as xr From cbf2f90571e031ac2619e553e25e8aa2adc5c8da Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Mon, 30 Mar 2026 11:26:33 -0700 Subject: [PATCH 2/2] Address review: normalize CPU LCC early-return path, move constants (#1088) - Add _norm_lon_rad to CPU _lcc_inv_point early-return (rho < 1e-30) for consistency with the CUDA fix - Move _WGS84_E2/_WGS84_A constants to top of test file --- xrspatial/reproject/_projections.py | 2 +- xrspatial/tests/test_reproject.py | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/xrspatial/reproject/_projections.py b/xrspatial/reproject/_projections.py index 2cd9883a..b84253df 100644 --- a/xrspatial/reproject/_projections.py +++ b/xrspatial/reproject/_projections.py @@ -417,7 +417,7 @@ def _lcc_inv_point(x, y, lon0, n, c, rho0, k0, e, a): rho = math.hypot(x, rho0_y) lam_n = math.atan2(x, rho0_y) if abs(rho) < 1e-30: - return math.degrees(lon0 + lam_n / n), 90.0 if n > 0 else -90.0 + return math.degrees(_norm_lon_rad(lon0 + lam_n / n)), 90.0 if n > 0 else -90.0 ts = math.pow(rho / (a * k0 * c), 1.0 / n) # Recover phi from ts via Newton (pj_sinhpsi2tanphi) phi_approx = math.pi / 2.0 - 2.0 * math.atan(ts) diff --git a/xrspatial/tests/test_reproject.py b/xrspatial/tests/test_reproject.py index fd86351a..96c67423 100644 --- a/xrspatial/tests/test_reproject.py +++ b/xrspatial/tests/test_reproject.py @@ -27,6 +27,10 @@ not HAS_PYPROJ, reason="pyproj required for reproject tests" ) +# WGS84 constants for projection round-trip tests +_WGS84_E2 = 2.0 * (1.0 / 298.257223563) - (1.0 / 298.257223563) ** 2 +_WGS84_A = 6378137.0 + # --------------------------------------------------------------------------- # Helpers @@ -1365,10 +1369,6 @@ def test_lcc_round_trip_stays_in_range(self): assert abs(lat_out - lat_in) < 1e-6 -_WGS84_E2 = 2.0 * (1.0 / 298.257223563) - (1.0 / 298.257223563) ** 2 -_WGS84_A = 6378137.0 - - class TestReprojWithLiteCRS: def test_reproject_wgs84_to_utm_with_lite_crs(self): import xarray as xr