diff --git a/.claude/sweep-reference-validation-state.csv b/.claude/sweep-reference-validation-state.csv index 359f6f59f..8e67a3076 100644 --- a/.claude/sweep-reference-validation-state.csv +++ b/.claude/sweep-reference-validation-state.csv @@ -3,3 +3,4 @@ convolution,2026-07-02,3619,MEDIUM,CONVENTION-DIFF,0.0 exact (float64),scipy 1.1 edge_detection,2026-07-18,3677,,MATCHES,"1e-6 rel (observed ~1e-13 interior, bit-exact vs ndi.correlate)",scipy 1.16.1; sobel_x/sobel_y/prewitt_x/prewitt_y/laplacian interior match ndi.sobel/prewitt/laplace and ndi.correlate with same stencils; boundary nan/nearest/reflect/wrap == scipy modes nearest/reflect/wrap incl. edges; tilted plane analytic exact (sobel=8g prewitt=6g laplacian=0); coords/cellsize independent as documented; cupy+dask+numpy+dask+cupy bit-identical to numpy (CUDA run); golden-value test added via #3677 geotiff,2026-07-02,,,MATCHES,rtol=1e-6 float64; float32 codecs exact,"rasterio 1.4.4 (GDAL 3.10.3); gdal-cli 3.11.4; tifffile 2026.3.3; richdem-unavailable; osgeo-python-unavailable. Reader/writer round-trip vs rasterio/GDAL: codecs zstd/deflate/lzw/packbits/none + predictor 2&3 dmax=0; CRS EPSG:32611 & 4326 round-trip exact; GeoTransform north-up, non-square px, tiepoint & ModelTransformation all correct; nodata float-NaN->sentinel + int16/uint8/int32/uint16 sentinels exact; BigTIFF magic+data ok; windowed reads match rasterio.Window incl coord alignment; multiband round-trips (band-last is xrspatial's documented convention); COG valid per 'rio cogeo validate'; overview mean matches numpy block-mean to float32 ulp; tifffile tag check: Compression/Predictor/SampleFormat/GeoKeyDir all correct. Existing repo parity suite parity/test_reference.py 31 passed (incl dask+cupy). No divergence; golden_corpus already pins parity." pathfinding,2026-07-08,3655,,MATCHES,1e-12 rel (observed <=2e-15),"a_star_search vs scipy 1.16.1 sparse.csgraph.dijkstra + skimage 0.26.0 MCP_Geometric + analytic open-grid truth; 64x64 open/maze/friction/NaN-friction/anisotropic inputs, 4+8 conn; path chains adjacency+increment consistent; backend parity numpy/dask+numpy/cupy/dask+cupy confirmed (CUDA host); no prior golden test -> added golden-value tests (#3655); gdal/richdem/wbt not applicable to pathfinding" +surface_distance,2026-08-16,3726,MEDIUM,MATCHES,1e-6 rel (observed 0.0 exact vs scipy on all backends),"scipy 1.16.1 (sparse.csgraph.dijkstra over an explicit grid graph with sqrt(horizontal^2+dz^2) edge weights) used as the reference; gdal-unavailable richdem-unavailable whitebox-unavailable; skimage 0.26.0 MCP_Geometric not applicable (per-pixel cost model, not 3D edge length). cuda True. numpy, dask+numpy (bounded map_overlap and unbounded iterative), cupy and dask+cupy all reproduce the scipy reference bit for bit (max_abs=0.0) on a Gaussian hill, a tilted plane, a seeded rough surface and a NaN-barrier scene, at connectivity 4 and 8, with and without max_distance and target_values. Analytic tilted-plane ground truth matches to float32 eps (3.9e-8 rel). surface_direction bearings match the reference bearing convention on north-up rasters (N=360, E=90). Geodesic mode matches an independent haversine graph at R=6378137 exactly; that is the library-wide radius (proximity uses the same), and the spherical-vs-WGS84-ellipsoid difference is +-0.17% locally, +0.003% mean bias - convention, not a defect, though the docstring does not state the sphere. Two CRITICAL defects surfaced during validation were already filed by parallel sweep agents, so no duplicate issue was opened: #3722 (heap sized H*W overflows, glibc abort; reproduced here at 40x40 with a 50% dense source mask, peak occupancy 1.07x H*W) and #3721 (cupy max_iterations=H+W truncates relaxation). Added new evidence to #3721: with no barriers at all, only steep ridges, cupy returns finite distances up to 173.8% too large on 1314/1600 pixels - silent, not NaN. #3723 flagged as a duplicate of #3722. Filed #3726 (MEDIUM): surface_allocation/surface_direction resolve exact ties differently per backend (29/4096 px numpy vs cupy, 25/4096 numpy vs dask) with no documented rule; PR documents the tie-break as unspecified and adds golden-value reference-parity tests (no golden test existed before)." diff --git a/xrspatial/surface_distance.py b/xrspatial/surface_distance.py index c91aeb13e..94ebb0c70 100644 --- a/xrspatial/surface_distance.py +++ b/xrspatial/surface_distance.py @@ -1505,6 +1505,15 @@ def surface_allocation( ------- xr.DataArray or xr.Dataset 2-D array of allocation values (float32). + + Notes + ----- + When two or more targets are exactly equidistant from a pixel, which + one is reported is unspecified and varies by backend: the numpy and + dask paths break the tie on Dijkstra's pop order, the cupy path on + the order its relaxation passes converge. The distances themselves + agree; only the reported source differs. Do not depend on the + tie-break. """ result_data = _compute( raster, elevation, x, y, target_values, max_distance, @@ -1548,6 +1557,12 @@ def surface_direction( ------- xr.DataArray or xr.Dataset 2-D array of direction values (float32, degrees). + + Notes + ----- + When two or more targets are exactly equidistant from a pixel, which + one the bearing points at is unspecified and varies by backend. See + :func:`surface_allocation`. """ result_data = _compute( raster, elevation, x, y, target_values, max_distance, diff --git a/xrspatial/tests/test_surface_distance.py b/xrspatial/tests/test_surface_distance.py index 74ae62fbd..9e84b7bf9 100644 --- a/xrspatial/tests/test_surface_distance.py +++ b/xrspatial/tests/test_surface_distance.py @@ -792,3 +792,149 @@ def test_error_message_mentions_grid_size(self): surface_distance(raster, elevation) with pytest.raises(MemoryError, match="dask"): surface_distance(raster, elevation) + + +# --------------------------------------------------------------------------- +# Reference parity — golden values from scipy.sparse.csgraph.dijkstra +# --------------------------------------------------------------------------- +# +# The expected arrays below were produced by an independent reference: +# scipy 1.16.1 ``scipy.sparse.csgraph.dijkstra`` run over an explicitly +# built grid graph whose edge weights are ``sqrt(horizontal**2 + dz**2)``, +# the cost model surface_distance documents. xrspatial reproduced them +# bit for bit (max abs difference 0.0) across the numpy, dask+numpy, cupy +# and dask+cupy backends. These tests pin that parity so it cannot break +# silently; regenerate them only against the reference, never against +# xrspatial's own output. + +_REF_ELEV = np.array([ + [0., 10., 25., 40., 55., 40., 25., 10.], + [0., 12., 30., 60., 75., 60., 30., 12.], + [0., 15., 40., 90., 120., 90., 40., 15.], + [0., 15., 40., np.nan, np.nan, 90., 40., 15.], + [0., 12., 30., np.nan, np.nan, 60., 30., 12.], + [0., 10., 25., 40., 55., 40., 25., 10.], + [0., 8., 18., 28., 36., 28., 18., 8.], + [0., 5., 10., 15., 20., 15., 10., 5.], +]) + +_REF_SOURCE = np.zeros((8, 8), dtype=np.float64) +_REF_SOURCE[0, 0] = 1.0 +_REF_SOURCE[7, 7] = 2.0 + +_REF_CELLSIZE = 30.0 + +_REF_DISTANCE_8 = np.array([ + [0.0000, 31.6228, 65.1638, 98.7048, 132.2458, 165.7868, 199.3279, 210.6487], # noqa: E501 + [30.0000, 44.0908, 78.5269, 120.1638, 153.7048, 174.9659, 195.4324, 180.5821], # noqa: E501 + [60.0000, 74.2404, 94.9239, 152.0116, 193.6485, 198.7048, 167.2159, 150.4324], # noqa: E501 + [90.0000, 104.2404, 123.4847, np.nan, np.nan, 170.6465, 137.2159, 120.4324], + [120.0000, 134.0908, 149.2404, np.nan, np.nan, 128.2201, 105.5931, 90.2828], + [150.0000, 163.5890, 176.9478, 149.6591, 125.7057, 92.1646, 75.1793, 60.2162], # noqa: E501 + [180.0000, 193.1741, 164.1876, 134.4155, 107.0445, 75.9962, 44.3734, 30.1496], # noqa: E501 + [210.0000, 182.4829, 152.0691, 121.6553, 91.2414, 60.8276, 30.4138, 0.0000], +]) + +_REF_DISTANCE_4 = np.array([ + [0.0000, 31.6228, 65.1638, 98.7048, 132.2458, 165.7868, 199.3279, 210.6487], # noqa: E501 + [30.0000, 61.6894, 95.5776, 134.7603, 168.3013, 201.8424, 215.5678, 180.5821], # noqa: E501 + [60.0000, 91.8390, 127.2004, 177.1867, 219.6131, 233.8779, 184.3046, 150.4324], # noqa: E501 + [90.0000, 121.8390, 157.2004, np.nan, np.nan, 203.8779, 154.3046, 120.4324], + [120.0000, 151.9886, 186.9743, np.nan, np.nan, 161.4514, 122.6818, 90.2828], + [150.0000, 181.6228, 213.9233, 186.6618, 158.9370, 125.3959, 92.2680, 60.2162], # noqa: E501 + [180.0000, 211.0484, 183.1174, 154.3508, 124.1333, 93.0849, 61.4622, 30.1496], # noqa: E501 + [210.0000, 182.4829, 152.0691, 121.6553, 91.2414, 60.8276, 30.4138, 0.0000], +]) + +_REF_ALLOCATION_8 = np.array([ + [1., 1., 1., 1., 1., 1., 1., 2.], + [1., 1., 1., 1., 1., 1., 2., 2.], + [1., 1., 1., 1., 1., 1., 2., 2.], + [1., 1., 1., np.nan, np.nan, 2., 2., 2.], + [1., 1., 1., np.nan, np.nan, 2., 2., 2.], + [1., 1., 2., 2., 2., 2., 2., 2.], + [1., 1., 2., 2., 2., 2., 2., 2.], + [1., 2., 2., 2., 2., 2., 2., 2.], +]) + + +@pytest.mark.parametrize("connectivity,expected", + [(8, _REF_DISTANCE_8), (4, _REF_DISTANCE_4)]) +@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy']) +def test_reference_parity_scipy_csgraph(backend, connectivity, expected): + """Distances match scipy.sparse.csgraph.dijkstra on the same graph.""" + raster = _make_raster(_REF_SOURCE, backend=backend, chunks=(4, 4), + res=_REF_CELLSIZE) + elevation = _make_raster(_REF_ELEV, backend=backend, chunks=(4, 4), + res=_REF_CELLSIZE) + + result = _compute(surface_distance(raster, elevation, + connectivity=connectivity)) + + np.testing.assert_array_equal(np.isnan(result), np.isnan(expected)) + np.testing.assert_allclose(result, expected, rtol=1e-5, equal_nan=True) + + +@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy']) +def test_reference_parity_allocation(backend): + """Allocation matches the scipy reference's nearest-source labels.""" + raster = _make_raster(_REF_SOURCE, backend=backend, chunks=(4, 4), + res=_REF_CELLSIZE) + elevation = _make_raster(_REF_ELEV, backend=backend, chunks=(4, 4), + res=_REF_CELLSIZE) + + result = _compute(surface_allocation(raster, elevation)) + + np.testing.assert_array_equal(np.isnan(result), + np.isnan(_REF_ALLOCATION_8)) + np.testing.assert_allclose(result, _REF_ALLOCATION_8, equal_nan=True) + + +@pytest.mark.skipif(not has_cuda_and_cupy(), reason="cupy/cuda not available") +def test_reference_parity_cupy(): + """The cupy backend reproduces the same reference distances.""" + raster = _make_raster(_REF_SOURCE, backend='cupy', res=_REF_CELLSIZE) + elevation = _make_raster(_REF_ELEV, backend='cupy', res=_REF_CELLSIZE) + + result = _compute(surface_distance(raster, elevation)) + + np.testing.assert_allclose(result, _REF_DISTANCE_8, rtol=1e-5, + equal_nan=True) + + +@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy']) +def test_tilted_plane_matches_analytic_surface_distance(backend): + """A tilted plane has a closed-form surface distance. + + On a plane dipping at a constant gradient along +x, the true distance + along the surface from the origin is ``sqrt(dx**2 + dy**2 + dz**2)``. + The grid graph reproduces it exactly along the row, column and main + diagonal through the source, where the octile path is the straight + line, so no metrication error enters. + """ + n = 12 + cellsize = 30.0 + gradient = 0.5 + + _, xx = np.mgrid[0:n, 0:n].astype(np.float64) + elev = gradient * xx * cellsize + source = np.zeros((n, n), dtype=np.float64) + source[0, 0] = 1.0 + + raster = _make_raster(source, backend=backend, chunks=(6, 6), + res=cellsize) + elevation = _make_raster(elev, backend=backend, chunks=(6, 6), + res=cellsize) + result = _compute(surface_distance(raster, elevation)) + + k = np.arange(n, dtype=np.float64) + # Straight down-dip along row 0. + np.testing.assert_allclose( + result[0, :], k * cellsize * np.sqrt(1 + gradient ** 2), rtol=1e-5) + # Across-dip along column 0: the plane is level in y. + np.testing.assert_allclose(result[:, 0], k * cellsize, rtol=1e-5) + # Main diagonal. + np.testing.assert_allclose( + result[np.arange(n), np.arange(n)], + np.sqrt(2 * (k * cellsize) ** 2 + (gradient * k * cellsize) ** 2), + rtol=1e-5)