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.

Author: brendancol

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
 

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,

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