Pin surface_distance against scipy csgraph and document the allocation tie-break

.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)."

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,

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)