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/_itrf.py b/xrspatial/reproject/_itrf.py index 86de16d40..bcec14ea8 100644 --- a/xrspatial/reproject/_itrf.py +++ b/xrspatial/reproject/_itrf.py @@ -136,8 +136,14 @@ def _helmert14_point(X, Y, Z, t_epoch, t_obs, position_vector): """Apply 14-parameter Helmert transform to a single ECEF point. - Parameters are in metres (translations), ppb (scale), and - arcseconds (rotations). Rates are per year. + Parameters are in metres (translations), ppm (scale), and + arcseconds (rotations). Rates are per year. The units match + PROJ's ``+proj=helmert`` definition, which is what the ITRF data + files this module parses are written against: the published + ppb-magnitude ITRF scale terms appear in the files as ppm values + (e.g. ITRF2014->ITRF2000 carries ``+s=0.00212`` ppm = 2.12 ppb). + Treating them as ppb made the scale term 1000x too small, a ~2 cm + error at the Earth's surface (GH #3276). """ dt = t_obs - t_epoch @@ -145,7 +151,7 @@ def _helmert14_point(X, Y, Z, tx_e = tx + dtx * dt ty_e = ty + dty * dt tz_e = tz + dtz * dt - s_e = 1.0 + (s + ds * dt) * 1e-9 # ppb -> scale factor + s_e = 1.0 + (s + ds * dt) * 1e-6 # ppm -> scale factor # Rotations: arcsec -> radians AS2RAD = math.pi / (180.0 * 3600.0) rx_e = (rx + drx * dt) * AS2RAD diff --git a/xrspatial/tests/test_reproject_itrf_scale_3276.py b/xrspatial/tests/test_reproject_itrf_scale_3276.py new file mode 100644 index 000000000..fd462e567 --- /dev/null +++ b/xrspatial/tests/test_reproject_itrf_scale_3276.py @@ -0,0 +1,176 @@ +"""Regression tests for GH #3276. + +``_helmert14_point`` converted the Helmert scale with ``* 1e-9`` +("ppb"), but the parameters come from PROJ's ITRF data files and PROJ's +``+proj=helmert`` defines ``+s`` / ``+ds`` in ppm. The published +ppb-magnitude ITRF scale terms appear in the files as ppm values +(ITRF2014->ITRF2000 carries ``+s=0.00212`` ppm = 2.12 ppb), so the old +factor made the scale term 1000x too small: ~23 mm of height error for +that pair at epoch 2024, on a function whose stated purpose is mm-level +frame transforms. +""" +from __future__ import annotations + +import math + +import numpy as np +import pytest + +try: + import pyproj + HAS_PYPROJ = True +except ImportError: + HAS_PYPROJ = False + +pytestmark = pytest.mark.skipif( + not HAS_PYPROJ, reason="pyproj required for ITRF tests" +) + + +def _direct_params(src, tgt): + """Raw helmert params for a direct (non-reversed) src->tgt entry.""" + from xrspatial.reproject._itrf import _find_transform + + params, is_reverse = _find_transform(src, tgt) + if params is None or is_reverse: + pytest.skip(f"no direct {src}->{tgt} entry in the PROJ data files") + return params + + +def _pyproj_pipeline(params): + """Build a 4D cart->helmert->cart pipeline from a raw param dict.""" + keys = ('x', 'y', 'z', 's', 'rx', 'ry', 'rz', + 'dx', 'dy', 'dz', 'ds', 'drx', 'dry', 'drz', 't_epoch') + parts = [f"+{k}={params[k]}" for k in keys if k in params] + convention = params.get('convention', 'position_vector') + pipeline = ( + "+proj=pipeline " + "+step +proj=cart +ellps=WGS84 " + f"+step +proj=helmert {' '.join(parts)} +convention={convention} " + "+step +inv +proj=cart +ellps=WGS84" + ) + return pyproj.Transformer.from_pipeline(pipeline) + + +class TestItrfMatchesPyprojHelmert: + @pytest.mark.parametrize("src,tgt", [ + ("ITRF2014", "ITRF2000"), + ("ITRF2014", "ITRF2008"), + ("ITRF2020", "ITRF2014"), + ]) + def test_direct_pair_matches_pyproj_4d(self, src, tgt): + """itrf_transform vs a pyproj pipeline from the same file entry. + + pyproj only applies the rate terms when given 4D coordinates, + so the reference transform passes the epoch as the 4th value. + Sub-mm agreement requires the ppm scale interpretation; the old + ppb factor left a ~23 mm height gap on ITRF2014->ITRF2000. + """ + from xrspatial.reproject import itrf_transform + + params = _direct_params(src, tgt) + tr = _pyproj_pipeline(params) + + lon, lat, h = -74.0, 40.7, 100.0 + epoch = 2024.0 + rx, ry, rz, _ = tr.transform(lon, lat, h, epoch) + ox, oy, oz = itrf_transform(lon, lat, h, src=src, tgt=tgt, epoch=epoch) + + de = (ox - rx) * 111320.0 * math.cos(math.radians(lat)) + dn = (oy - ry) * 111320.0 + dh = oz - rz + err_mm = math.sqrt(de * de + dn * dn + dh * dh) * 1000.0 + assert err_mm < 0.5, f"{src}->{tgt} differs from pyproj by {err_mm:.3f} mm" + + def test_scale_term_not_dropped(self): + """The scale term must contribute at its ppm magnitude. + + ITRF2014->ITRF2000 has s + ds*(2024-2010) = 0.00366 ppm, which + is ~23 mm radially at the Earth's surface. With the broken ppb + factor the contribution collapsed to ~23 um, so zeroing the + scale barely changed the output. Compare the real transform + against one with s and ds removed and require the difference to + show up. + """ + from xrspatial.reproject import itrf_transform + + params = _direct_params("ITRF2014", "ITRF2000") + s = params.get('s', 0.0) + ds = params.get('ds', 0.0) + t_epoch = params.get('t_epoch', 2010.0) + epoch = 2024.0 + s_eff_ppm = s + ds * (epoch - t_epoch) + if abs(s_eff_ppm) < 1e-4: + pytest.skip("scale term too small in this PROJ data version") + + lon, lat, h = -74.0, 40.7, 100.0 + ox, oy, oz = itrf_transform(lon, lat, h, src='ITRF2014', + tgt='ITRF2000', epoch=epoch) + # Reference without scale: pyproj pipeline with s/ds zeroed. + params_noscale = dict(params) + params_noscale['s'] = 0.0 + params_noscale['ds'] = 0.0 + tr = _pyproj_pipeline(params_noscale) + nx, ny, nz, _ = tr.transform(lon, lat, h, epoch) + + de = (ox - nx) * 111320.0 * math.cos(math.radians(lat)) + dn = (oy - ny) * 111320.0 + dh = oz - nz + shift_mm = math.sqrt(de * de + dn * dn + dh * dh) * 1000.0 + # Expected radial effect: s_eff * R. R = 6.37e6 m approximates + # the geocentric radius at the test point (~6368 km); the 5% + # tolerance absorbs that approximation. With the broken ppb + # factor the measured shift is 1000x smaller and fails by three + # orders of magnitude. + expected_mm = abs(s_eff_ppm) * 1e-6 * 6.37e6 * 1000.0 + assert shift_mm == pytest.approx(expected_mm, rel=0.05) + + def test_reversed_lookup_matches_pyproj_inverse(self): + """Reversed lookups negate the file parameters, scale included. + + ITRF2000 -> ITRF2014 has no direct file entry; _find_transform + negates the ITRF2014 -> ITRF2000 parameters. The reference is + the same pipeline with an inverse helmert step. + """ + from xrspatial.reproject import itrf_transform + from xrspatial.reproject._itrf import _find_transform + + params, is_reverse = _find_transform("ITRF2000", "ITRF2014") + if params is None or not is_reverse: + pytest.skip("ITRF2000->ITRF2014 is not a reversed lookup here") + keys = ('x', 'y', 'z', 's', 'rx', 'ry', 'rz', + 'dx', 'dy', 'dz', 'ds', 'drx', 'dry', 'drz', 't_epoch') + parts = [f"+{k}={params[k]}" for k in keys if k in params] + convention = params.get('convention', 'position_vector') + tr = pyproj.Transformer.from_pipeline( + "+proj=pipeline " + "+step +proj=cart +ellps=WGS84 " + f"+step +inv +proj=helmert {' '.join(parts)} " + f"+convention={convention} " + "+step +inv +proj=cart +ellps=WGS84" + ) + lon, lat, h = -74.0, 40.7, 100.0 + epoch = 2024.0 + rx, ry, rz, _ = tr.transform(lon, lat, h, epoch) + ox, oy, oz = itrf_transform(lon, lat, h, src='ITRF2000', + tgt='ITRF2014', epoch=epoch) + de = (ox - rx) * 111320.0 * math.cos(math.radians(lat)) + dn = (oy - ry) * 111320.0 + dh = oz - rz + err_mm = math.sqrt(de * de + dn * dn + dh * dh) * 1000.0 + assert err_mm < 0.5 + + def test_array_and_scalar_agree(self): + from xrspatial.reproject import itrf_transform + + _direct_params("ITRF2014", "ITRF2000") + lons = np.array([-74.0, 10.0]) + lats = np.array([40.7, 50.0]) + hs = np.array([100.0, 0.0]) + alon, alat, ah = itrf_transform(lons, lats, hs, src='ITRF2014', + tgt='ITRF2000', epoch=2024.0) + slon, slat, sh = itrf_transform(-74.0, 40.7, 100.0, src='ITRF2014', + tgt='ITRF2000', epoch=2024.0) + assert alon[0] == pytest.approx(slon, abs=1e-12) + assert alat[0] == pytest.approx(slat, abs=1e-12) + assert ah[0] == pytest.approx(sh, abs=1e-9)