From 39f13a8495b3a60569771226ee620b490e7ce6e9 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Sun, 16 Aug 2026 12:17:47 -0400 Subject: [PATCH 1/9] Document surface_distance backends, geodesic limits, and add examples (#3714) The three public functions in surface_distance.py had no runnable example, no statement of which array backends they support, and no note that method='geodesic' works on NumPy input only. - State NumPy / CuPy / Dask+NumPy / Dask+CuPy support, matching the wording cost_distance() already uses. - Say on the method parameter that geodesic requires a NumPy-backed DataArray; the other three backends raise NotImplementedError. - Add a .. sourcecode:: python Examples block to each function, output pinned to a real run on this host. - Note on surface_direction() that the unbounded dask path currently returns wrong bearings, and link the fix to #3713. - Add docstring guard tests plus checks that the pinned example output is what the functions return. --- xrspatial/surface_distance.py | 141 ++++++++++++++++++++++- xrspatial/tests/test_surface_distance.py | 129 +++++++++++++++++++++ 2 files changed, 267 insertions(+), 3 deletions(-) diff --git a/xrspatial/surface_distance.py b/xrspatial/surface_distance.py index c91aeb13e..15f0796aa 100644 --- a/xrspatial/surface_distance.py +++ b/xrspatial/surface_distance.py @@ -1433,6 +1433,13 @@ def surface_distance( Edge cost accounts for both horizontal distance and elevation change: ``sqrt(horizontal_dist^2 + dz^2)``. + Surface distance supports NumPy, CuPy, Dask with NumPy, and Dask with + CuPy backed xarray DataArray. The return value is of the same type as + the input: a NumPy-backed input gives a NumPy-backed result, a + CuPy-backed input gives a CuPy-backed result, and a Dask-backed input + gives a Dask-backed result. ``method='geodesic'`` is the exception -- + it is implemented for NumPy-backed input only. + Parameters ---------- raster : xr.DataArray or xr.Dataset @@ -1456,13 +1463,53 @@ def surface_distance( method : str, default='planar' ``'planar'`` uses cell sizes in map units. ``'geodesic'`` computes great-circle horizontal distances - from lat/lon coordinates (elevation in meters). + from lat/lon coordinates (elevation in meters). Geodesic mode + requires a NumPy-backed DataArray; CuPy, Dask with NumPy, and + Dask with CuPy input raises ``NotImplementedError``. Returns ------- xr.DataArray or xr.Dataset 2-D array of surface distance values (float32). Source pixels have distance 0. Unreachable pixels are NaN. + + Examples + -------- + A single source in the top-left corner and a peak in the middle of an + otherwise flat 3x3 grid: + + .. sourcecode:: python + + >>> import numpy as np + >>> import xarray as xr + >>> from xrspatial import surface_distance + >>> source = np.array([ + ... [1., 0., 0.], + ... [0., 0., 0.], + ... [0., 0., 0.], + ... ]) + >>> elevation = np.array([ + ... [0., 0., 0.], + ... [0., 3., 0.], + ... [0., 0., 0.], + ... ]) + >>> n, m = source.shape + >>> raster = xr.DataArray(source, dims=['y', 'x'], name='raster') + >>> raster['y'] = np.arange(n)[::-1] + >>> raster['x'] = np.arange(m) + >>> elev = xr.DataArray( + ... elevation, dims=['y', 'x'], name='elevation') + >>> elev['y'] = np.arange(n)[::-1] + >>> elev['x'] = np.arange(m) + + >>> surface_distance(raster, elev).values + array([[0. , 1. , 2. ], + [1. , 3.3166249, 2.4142137], + [2. , 2.4142137, 3.4142137]], dtype=float32) + + Climbing to the summit costs ``sqrt(2 + 9) = 3.32``, while the far + corner is reached for 3.41 by walking around the peak along the flat + edges rather than over it. """ result_data = _compute( raster, elevation, x, y, target_values, max_distance, @@ -1492,6 +1539,11 @@ def surface_allocation( For each pixel, returns the value of the nearest target pixel by surface distance through the elevation model. + Surface allocation supports NumPy, CuPy, Dask with NumPy, and Dask + with CuPy backed xarray DataArray, and returns the same type it was + given. ``method='geodesic'`` is implemented for NumPy-backed input + only. + Parameters ---------- raster : xr.DataArray or xr.Dataset @@ -1504,7 +1556,41 @@ def surface_allocation( Returns ------- xr.DataArray or xr.Dataset - 2-D array of allocation values (float32). + 2-D array of allocation values (float32). Pixels that cannot + reach any target are NaN. + + Examples + -------- + Two sources on opposite corners, separated by a peak in the middle: + + .. sourcecode:: python + + >>> import numpy as np + >>> import xarray as xr + >>> from xrspatial import surface_allocation + >>> source = np.array([ + ... [1., 0., 0.], + ... [0., 0., 0.], + ... [0., 2., 0.], + ... ]) + >>> elevation = np.array([ + ... [0., 0., 0.], + ... [0., 3., 0.], + ... [0., 0., 0.], + ... ]) + >>> n, m = source.shape + >>> raster = xr.DataArray(source, dims=['y', 'x'], name='raster') + >>> raster['y'] = np.arange(n)[::-1] + >>> raster['x'] = np.arange(m) + >>> elev = xr.DataArray( + ... elevation, dims=['y', 'x'], name='elevation') + >>> elev['y'] = np.arange(n)[::-1] + >>> elev['x'] = np.arange(m) + + >>> surface_allocation(raster, elev).values + array([[1., 1., 1.], + [1., 2., 2.], + [2., 2., 2.]], dtype=float32) """ result_data = _compute( raster, elevation, x, y, target_values, max_distance, @@ -1535,6 +1621,11 @@ def surface_direction( nearest target pixel by surface distance. 0 = source pixel, 90 = east, 180 = south, 270 = west, 360 = north. + Surface direction supports NumPy, CuPy, Dask with NumPy, and Dask + with CuPy backed xarray DataArray, and returns the same type it was + given. ``method='geodesic'`` is implemented for NumPy-backed input + only. + Parameters ---------- raster : xr.DataArray or xr.Dataset @@ -1547,7 +1638,51 @@ def surface_direction( Returns ------- xr.DataArray or xr.Dataset - 2-D array of direction values (float32, degrees). + 2-D array of direction values (float32, degrees). Pixels that + cannot reach any target are NaN. + + Notes + ----- + On a Dask-backed input with an infinite ``max_distance``, or one large + enough that the search radius exceeds the chunk size, the bearings are + currently measured from the wrong origin in every chunk except the + first. Track the fix in + https://github.com/xarray-contrib/xarray-spatial/issues/3713. Setting + a finite ``max_distance`` smaller than the chunk size takes the + ``map_overlap`` path, which is unaffected. + + Examples + -------- + A single source at the centre of a flat 3x3 grid. Each pixel reports + the compass bearing back toward that source: + + .. sourcecode:: python + + >>> import numpy as np + >>> import xarray as xr + >>> from xrspatial import surface_direction + >>> source = np.array([ + ... [0., 0., 0.], + ... [0., 1., 0.], + ... [0., 0., 0.], + ... ]) + >>> elevation = np.zeros((3, 3)) + >>> n, m = source.shape + >>> raster = xr.DataArray(source, dims=['y', 'x'], name='raster') + >>> raster['y'] = np.arange(n)[::-1] + >>> raster['x'] = np.arange(m) + >>> elev = xr.DataArray( + ... elevation, dims=['y', 'x'], name='elevation') + >>> elev['y'] = np.arange(n)[::-1] + >>> elev['x'] = np.arange(m) + + >>> surface_direction(raster, elev).values + array([[135., 180., 225.], + [ 90., 0., 270.], + [ 45., 360., 315.]], dtype=float32) + + The pixel east of the source reads 270 (west) because the bearing + points back at the source, and the source itself reads 0. """ 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..24f27a062 100644 --- a/xrspatial/tests/test_surface_distance.py +++ b/xrspatial/tests/test_surface_distance.py @@ -1,5 +1,7 @@ """Tests for xrspatial.surface_distance.""" +import inspect + import numpy as np import pytest import xarray as xr @@ -792,3 +794,130 @@ def test_error_message_mentions_grid_size(self): surface_distance(raster, elevation) with pytest.raises(MemoryError, match="dask"): surface_distance(raster, elevation) + + +# --------------------------------------------------------------------------- +# Tests — docstring contract +# --------------------------------------------------------------------------- + + +_PUBLIC = [surface_distance, surface_allocation, surface_direction] + + +def _flat_doc(func): + """Docstring with line wrapping collapsed, so pinned phrases survive.""" + return " ".join(inspect.getdoc(func).split()) + + +@pytest.mark.parametrize("func", _PUBLIC) +def test_docstring_has_examples_section(func): + """Every public surface-distance function ships a runnable example.""" + doc = inspect.getdoc(func) + assert any(ln.strip() == "Examples" for ln in doc.splitlines()), ( + f"{func.__name__} docstring has no Examples section" + ) + assert ">>>" in doc, f"{func.__name__} Examples section has no code" + # Three dots renders as literal text instead of a code block. + assert "... sourcecode::" not in doc + + +@pytest.mark.parametrize("func", _PUBLIC) +def test_docstring_states_all_backends(func): + """All three functions dispatch to numpy, cupy, dask+numpy, dask+cupy. + + The docstrings named no backend at all, leaving users to guess. Pin the + phrases so the claim cannot silently disappear; reword the assertions, + not the docs, if the wording is revised while staying accurate. + """ + doc = _flat_doc(func) + assert "CuPy" in doc + assert "Dask with NumPy" in doc + assert "Dask with CuPy" in doc + + +@pytest.mark.parametrize("func", _PUBLIC) +def test_docstring_notes_geodesic_is_numpy_only(func): + """geodesic raises NotImplementedError on every non-numpy backend.""" + doc = _flat_doc(func) + assert "NumPy-backed input only" in doc or ( + "requires a NumPy-backed DataArray" in doc + ) + + +def test_geodesic_rejects_dask(): + """The documented geodesic limitation matches what the code does.""" + if da is None: + pytest.skip("dask not installed") + + source = np.zeros((4, 4), dtype=np.float64) + source[1, 1] = 1.0 + elev = np.zeros((4, 4), dtype=np.float64) + raster = _make_raster(source, backend='dask+numpy', chunks=(2, 2)) + elevation = _make_raster(elev, backend='dask+numpy', chunks=(2, 2)) + + with pytest.raises(NotImplementedError, match="geodesic"): + surface_distance(raster, elevation, method='geodesic') + + +def _docstring_example_rasters(source, elevation): + """Build the y-descending DataArrays used by the Examples blocks.""" + n, m = source.shape + raster = xr.DataArray(source, dims=['y', 'x'], name='raster') + raster['y'] = np.arange(n)[::-1] + raster['x'] = np.arange(m) + elev = xr.DataArray(elevation, dims=['y', 'x'], name='elevation') + elev['y'] = np.arange(n)[::-1] + elev['x'] = np.arange(m) + return raster, elev + + +_PEAK_ELEVATION = np.array([ + [0., 0., 0.], + [0., 3., 0.], + [0., 0., 0.], +]) + + +@pytest.mark.parametrize("func, source, expected", [ + ( + surface_distance, + np.array([[1., 0., 0.], + [0., 0., 0.], + [0., 0., 0.]]), + np.array([[0., 1., 2.], + [1., 3.3166249, 2.4142137], + [2., 2.4142137, 3.4142137]], dtype=np.float32), + ), + ( + surface_allocation, + np.array([[1., 0., 0.], + [0., 0., 0.], + [0., 2., 0.]]), + np.array([[1., 1., 1.], + [1., 2., 2.], + [2., 2., 2.]], dtype=np.float32), + ), +]) +def test_docstring_example_matches_output(func, source, expected): + """The pinned output in each Examples block is what the code returns.""" + raster, elev = _docstring_example_rasters(source, _PEAK_ELEVATION) + result = func(raster, elev) + assert result.dtype == np.float32 + np.testing.assert_allclose(result.values, expected, rtol=1e-6) + + +def test_direction_docstring_example_matches_output(): + """surface_direction()'s pinned Examples output, on the numpy backend.""" + source = np.array([ + [0., 0., 0.], + [0., 1., 0.], + [0., 0., 0.], + ]) + raster, elev = _docstring_example_rasters(source, np.zeros((3, 3))) + + expected = np.array([[135., 180., 225.], + [90., 0., 270.], + [45., 360., 315.]], dtype=np.float32) + result = surface_direction(raster, elev) + assert result.dtype == np.float32 + np.testing.assert_allclose(result.values, expected, rtol=1e-6) From b9278f17c81d93b57557686ad8317ab7f5989f9f Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Sun, 16 Aug 2026 12:17:51 -0400 Subject: [PATCH 2/9] Record the surface_distance documentation sweep in the state CSV --- .claude/sweep-documentation-state.csv | 1 + 1 file changed, 1 insertion(+) diff --git a/.claude/sweep-documentation-state.csv b/.claude/sweep-documentation-state.csv index b72ec9f31..2c7a7950d 100644 --- a/.claude/sweep-documentation-state.csv +++ b/.claude/sweep-documentation-state.csv @@ -10,4 +10,5 @@ interpolate-idw,2026-06-26,,MEDIUM,1;5,"Cat1 MEDIUM: idw (only public func in _i mahalanobis,2026-06-30,3579,MEDIUM,1;2;5,"Cat1 MEDIUM: mahalanobis (only public func) had Parameters/Returns but no Examples section (peers normalize.rescale/standardize do). Cat5 MEDIUM: NaN propagation (any non-finite band -> NaN pixel) and auto-stats N+1 all-finite requirement undocumented. Cat2 LOW: name param default not noted. Doc-only fix on deep-sweep-documentation-mahalanobis-2026-06-30: added Notes (NaN+backends) and a runnable Examples block (executed, output matches incl. NaN), set name default='mahalanobis'. Returns float64/shape and 4-backend claim verified against ArrayTypeFunctionMapping (accurate, left as-is). Cat3 n/a (no prior examples). Cat4 clean (listed in reference/utilities.rst). issue #3579. CUDA available: ran numpy example; cupy/dask covered by test_mahalanobis.py (24 pass).",1/1 pathfinding,2026-07-08,3651,MEDIUM,1;3;5,"Cat1 MEDIUM: multi_stop_search had Parameters/Returns/Raises but no Examples section (a_star_search has 13 example lines); added runnable example pinning waypoint_order/segment_costs/total_cost. Cat3 MEDIUM: a_star_search Examples used '... sourcecode:: python' (three dots, only occurrence in codebase) so the rendered page shows literal text instead of a code block; fixed to '..'; example itself executes clean. Cat5 MEDIUM: NaN surface cells are always impassable (_is_not_crossable, pinned by input_data_with_nans tests) but neither docstring said so (friction NaN rule was documented); documented in summary + barriers entries of both funcs. LOW fixed in passing: 'y: str' missing numpydoc space. Deferred to #3644 (api-consistency sibling, same day): 'values to bin' surface description + empty connectivity description. Clean: params match both signatures (grouped-entry aware check), both funcs in reference/pathfinding.rst, backend table verified by running all 4 backends (CUDA available; dask path returns equally-optimal alternate path, goal cost identical). Added 3 docstring-validation tests to test_pathfinding.py. Fix on deep-sweep-documentation-pathfinding-2026-07-08, issue #3651.",2/2 perlin,2026-06-23,,MEDIUM,2;5,"name param undocumented (Cat2) + float-dtype requirement/ValueError undocumented, no Raises section (Cat5); fixed in deep-sweep-documentation-perlin-2026-06-23; repo has issues disabled so no issue number; example runs and output matches; 1 public func (perlin) listed in reference/surface.rst",1/1 +surface_distance,2026-08-16,3714;3713,HIGH,1;2;5,"severity_max HIGH is the separately-filed code bug #3713; the documentation findings themselves top out at MEDIUM. 3 MEDIUM: zero >>> example blocks in the whole module (Cat 1); no backend-support statement although all 4 backends dispatch and run (Cat 5); method='geodesic' documented with no caveat but raises NotImplementedError on cupy, dask+numpy and dask+cupy (Cat 5). Fixed in PR for #3714: backend paragraph, geodesic caveat, runnable Examples on all 3 public funcs, plus 13 docstring guard tests. 3 LOW, not fixed: surface_allocation/surface_direction fold 6 signature params into one 'x, y, target_values, ... See surface_distance' entry (numpydoc parses it, nothing misleading); the UserWarning emitted on the unbounded dask path is undocumented; pre-existing flake8 F841 at test_surface_distance.py:236 left alone. Cat 4 clean: all 3 funcs appear in docs/source/reference/proximity.rst. Separately filed #3713 (HIGH, code bug not doc bug): surface_direction returns wrong bearings on the dask iterative path because _finalize_direction subtracts block-local pixel indices from the global src_row/src_col written by _run_tile; 105/144 pixels wrong on a 12x12 chunks=(6,6) case. CUDA available, all 4 backends executed.",3/3 templates,2026-06-26,3541,MEDIUM,3;5,"Cat3/Cat5 MEDIUM: from_template docstring example 'from_template(""FRA"", preserve=""shape"").attrs[""crs""]' showed 32631/UTM 31N but code returns 32630/UTM 30N (GADM FRA bbox includes overseas territories, centroid lon -2.98 -> UTM 30N; code correct per documented centroid-UTM contract). Doc-only fix on deep-sweep-documentation-templates-2026-06-26; issue #3541, PR pending. Both public funcs (from_template, list_templates) fully documented and in reference/templates.rst (no Cat1/Cat2/Cat4). CUDA available: ran all docstring examples incl. preserve paths; 63/63 test_templates.py pass.",2/2 From 80de4c993a8734ce8b8965ffe30e4bd7afb1dbf0 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Sun, 16 Aug 2026 12:19:45 -0400 Subject: [PATCH 3/9] Close the surface_distance backend test gaps and fix dask direction (#3725, #3713) The test file covered three public functions across four backends with 37 tests, but never ran anything on dask+cupy and only ever called surface_direction() on eager numpy. That second hole was hiding a real bug: _run_tile() records global source indices while _finalize_direction() built its pixel grid from block-local ones, so every tile except the top-left measured its bearings from the wrong origin. On a flat 6x6 raster chunked (3, 3) the iterative path returned the top-left tile's bearing field tiled across the whole output. _finalize_direction() and _extract_output() now take row/col offsets, and _assemble_sd() passes the tile's own offsets so both sides of the subtraction live in the same coordinate space. The eager numpy and bounded map_overlap paths keep the default of 0 because they seed local indices. New tests cover the holes the audit found: dask+cupy on all three public functions, surface_direction parity across every backend on both the bounded and iterative branches, Inf and all-NaN elevation, 1x1 and Nx1 rasters, connectivity=4 off the eager path, target_values on cupy, Dataset input through @supports_dataset, the geodesic NotImplementedError per backend, the dims rejection, and attrs/coords/dim preservation. Tests 37 -> 106, all passing with CUDA present. Branch coverage of xrspatial.surface_distance 80% -> 87% under NUMBA_DISABLE_JIT=1. --- xrspatial/surface_distance.py | 22 +- xrspatial/tests/test_surface_distance.py | 414 +++++++++++++++++++++++ 2 files changed, 431 insertions(+), 5 deletions(-) diff --git a/xrspatial/surface_distance.py b/xrspatial/surface_distance.py index c91aeb13e..b2b3c8fcc 100644 --- a/xrspatial/surface_distance.py +++ b/xrspatial/surface_distance.py @@ -357,13 +357,20 @@ def _finalize_alloc(alloc, dist, max_distance): def _finalize_direction(src_row, src_col, dist, cellsize_x, cellsize_y, - max_distance): + max_distance, row_offset=0, col_offset=0): """Compute compass bearing from each pixel to its allocated source. Uses pixel index differences scaled by cell size. + + ``src_row``/``src_col`` may hold indices in a wider coordinate space + than the block being finalized -- the dask iterative path records + global indices while each block sees only its own shape. The + offsets shift the pixel grid into that same space so the difference + is taken between comparable indices. """ H, W = dist.shape - row_idx, col_idx = np.meshgrid(np.arange(H), np.arange(W), indexing='ij') + row_idx, col_idx = np.meshgrid(np.arange(H) + row_offset, + np.arange(W) + col_offset, indexing='ij') # Coordinate differences (source - pixel) dx = (src_col.astype(np.float64) - col_idx) * cellsize_x @@ -381,7 +388,8 @@ def _finalize_direction(src_row, src_col, dist, cellsize_x, cellsize_y, def _extract_output(dist, alloc, src_row, src_col, - cellsize_x, cellsize_y, max_distance, mode): + cellsize_x, cellsize_y, max_distance, mode, + row_offset=0, col_offset=0): """Select and finalize the requested output from raw Dijkstra arrays.""" if mode == DISTANCE: return _finalize_dist(dist, max_distance) @@ -389,7 +397,8 @@ def _extract_output(dist, alloc, src_row, src_col, return _finalize_alloc(alloc, dist, max_distance) else: return _finalize_direction(src_row, src_col, dist, - cellsize_x, cellsize_y, max_distance) + cellsize_x, cellsize_y, max_distance, + row_offset, col_offset) # --------------------------------------------------------------------------- @@ -1174,8 +1183,11 @@ def _tile_fn(source_block, elev_block, block_info=None): max_distance, dy, dx, dd, row_offset, col_offset, ) + # _run_tile records global src_row/src_col, so the pixel grid has + # to be built in global coordinates too. return _extract_output(dist, alloc_arr, srow, scol, - cellsize_x, cellsize_y, max_distance, mode) + cellsize_x, cellsize_y, max_distance, mode, + row_offset, col_offset) return da.map_blocks( _tile_fn, diff --git a/xrspatial/tests/test_surface_distance.py b/xrspatial/tests/test_surface_distance.py index 74ae62fbd..c3b2db616 100644 --- a/xrspatial/tests/test_surface_distance.py +++ b/xrspatial/tests/test_surface_distance.py @@ -1,5 +1,7 @@ """Tests for xrspatial.surface_distance.""" +import warnings + import numpy as np import pytest import xarray as xr @@ -792,3 +794,415 @@ def test_error_message_mentions_grid_size(self): surface_distance(raster, elevation) with pytest.raises(MemoryError, match="dask"): surface_distance(raster, elevation) + + +# --------------------------------------------------------------------------- +# Backend coverage — every public function on every backend +# --------------------------------------------------------------------------- + +ALL_BACKENDS = ['numpy', 'dask+numpy', 'cupy', 'dask+cupy'] + +# Bounded vs iterative: with 6x6 chunks, max_distance=4 gives an overlap +# depth of 5, which fits inside a chunk, so the dask backends take the +# map_overlap branch. max_distance=inf forces the iterative tile sweep. +_BOUNDED = 4.0 + + +def _needs(backend): + if backend in ('dask+numpy', 'dask+cupy') and da is None: + pytest.skip("dask not installed") + if backend in ('cupy', 'dask+cupy') and not has_cuda_and_cupy(): + pytest.skip("cupy/cuda not available") + + +def _parity_scene(): + """A scene with relief and two distinguishable sources.""" + source = np.zeros((12, 12), dtype=np.float64) + source[2, 3] = 1.0 + source[9, 8] = 2.0 + elev = np.random.default_rng(42).uniform(0, 60, (12, 12)) + return source, elev + + +@pytest.mark.parametrize( + "func", [surface_distance, surface_allocation, surface_direction], + ids=['distance', 'allocation', 'direction'], +) +@pytest.mark.parametrize("backend", ALL_BACKENDS[1:]) +def test_bounded_matches_numpy(func, backend): + """Bounded (map_overlap) path matches the numpy baseline everywhere. + + surface_direction had no backend coverage at all before this test. + """ + _needs(backend) + source, elev = _parity_scene() + + expected = _compute(func(_make_raster(source), _make_raster(elev), + max_distance=_BOUNDED)) + actual = _compute(func(_make_raster(source, backend, chunks=(6, 6)), + _make_raster(elev, backend, chunks=(6, 6)), + max_distance=_BOUNDED)) + + np.testing.assert_allclose(actual, expected, rtol=1e-5, equal_nan=True) + + +@pytest.mark.parametrize( + "func", [surface_distance, surface_allocation, surface_direction], + ids=['distance', 'allocation', 'direction'], +) +@pytest.mark.parametrize("backend", ['dask+numpy', 'dask+cupy']) +def test_iterative_matches_numpy_all_modes(func, backend): + """Unbounded (iterative tile Dijkstra) path matches numpy. + + Regression test for the direction mode: ``_run_tile`` records global + source indices, so ``_finalize_direction`` has to build its pixel + grid in global coordinates too. Before the fix every tile but the + top-left one reported a bearing measured from its own origin. + """ + _needs(backend) + source, elev = _parity_scene() + + expected = _compute(func(_make_raster(source), _make_raster(elev))) + with pytest.warns(UserWarning, match="iterative"): + actual = _compute(func(_make_raster(source, backend, chunks=(6, 6)), + _make_raster(elev, backend, chunks=(6, 6)))) + + np.testing.assert_allclose(actual, expected, rtol=1e-4, equal_nan=True) + + +def test_iterative_direction_uses_global_pixel_indices(): + """Explicit check of the bearings the iterative path produces. + + A single source in the top-left chunk of a flat 6x6 raster gives an + unambiguous bearing for every pixel. The bug returned the top-left + chunk's bearing field tiled across the whole raster, so (3, 0) read + as 0 ("you are the source") instead of 360 (source due north). + """ + if da is None: + pytest.skip("dask not installed") + + source = np.zeros((6, 6), dtype=np.float64) + source[0, 0] = 1.0 + elev = np.zeros((6, 6), dtype=np.float64) + + with pytest.warns(UserWarning, match="iterative"): + result = _compute(surface_direction( + _make_raster(source, 'dask+numpy', chunks=(3, 3)), + _make_raster(elev, 'dask+numpy', chunks=(3, 3)), + )) + + # Only the source itself reports 0. + assert result[0, 0] == 0.0 + assert np.count_nonzero(result == 0.0) == 1 + # Column 0 below the source: the source is due north. + np.testing.assert_allclose(result[1:, 0], 360.0, atol=1e-4) + # Row 0 right of the source: the source is due west. + np.testing.assert_allclose(result[0, 1:], 270.0, atol=1e-4) + # Perfect diagonal: north-west. + for i in range(1, 6): + assert result[i, i] == pytest.approx(315.0, abs=1e-4) + + +@pytest.mark.skipif(not has_cuda_and_cupy(), reason="cupy/cuda not available") +@pytest.mark.skipif(da is None, reason="dask not installed") +def test_dask_cupy_returns_dask_cupy_array(): + """dask+cupy input keeps both wrappers on the way out.""" + source = np.zeros((12, 12), dtype=np.float64) + source[5, 5] = 1.0 + elev = np.zeros((12, 12), dtype=np.float64) + + result = surface_distance( + _make_raster(source, 'dask+cupy', chunks=(6, 6)), + _make_raster(elev, 'dask+cupy', chunks=(6, 6)), + max_distance=_BOUNDED, + ) + assert isinstance(result.data, da.Array) + assert is_cupy_array(result.data.compute()) + + +@pytest.mark.skipif(not has_cuda_and_cupy(), reason="cupy/cuda not available") +def test_cupy_direction_returns_cupy_array(): + """The cupy direction branch converts its numpy result back to cupy.""" + source = np.zeros((5, 5), dtype=np.float64) + source[2, 2] = 1.0 + elev = np.zeros((5, 5), dtype=np.float64) + + result = surface_direction(_make_raster(source, 'cupy'), + _make_raster(elev, 'cupy')) + assert is_cupy_array(result.data) + + +@pytest.mark.parametrize("backend", ALL_BACKENDS) +def test_target_values_all_backends(backend): + """target_values filtering behaves the same on every backend.""" + _needs(backend) + source = np.zeros((6, 6), dtype=np.float64) + source[0, 0] = 1.0 + source[5, 5] = 2.0 + elev = np.zeros((6, 6), dtype=np.float64) + + expected = _compute(surface_distance(_make_raster(source), + _make_raster(elev), + target_values=[2])) + raster = _make_raster(source, backend, chunks=(3, 3)) + elevation = _make_raster(elev, backend, chunks=(3, 3)) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + actual = _compute(surface_distance(raster, elevation, + target_values=[2])) + + # Only the pixel holding value 2 is a source. + assert actual[5, 5] == 0.0 + assert actual[0, 0] == pytest.approx(5 * np.sqrt(2), abs=1e-4) + np.testing.assert_allclose(actual, expected, rtol=1e-5, equal_nan=True) + + +@pytest.mark.parametrize("backend", ALL_BACKENDS) +def test_no_sources_all_nan_all_backends(backend): + """An all-zero source raster yields NaN everywhere on every backend.""" + _needs(backend) + zeros = np.zeros((4, 4), dtype=np.float64) + raster = _make_raster(zeros, backend, chunks=(2, 2)) + elevation = _make_raster(zeros, backend, chunks=(2, 2)) + + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + assert np.all(np.isnan(_compute(surface_distance(raster, elevation)))) + assert np.all(np.isnan(_compute(surface_allocation(raster, + elevation)))) + assert np.all(np.isnan(_compute(surface_direction(raster, + elevation)))) + + +# --------------------------------------------------------------------------- +# Edge cases — Inf and all-NaN elevation +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize("backend", ALL_BACKENDS) +def test_inf_elevation_is_a_barrier(backend): + """Non-finite elevation blocks paths whether it is NaN or Inf.""" + _needs(backend) + source = np.zeros((1, 5), dtype=np.float64) + source[0, 0] = 1.0 + elev = np.zeros((1, 5), dtype=np.float64) + elev[0, 2] = np.inf + + raster = _make_raster(source, backend, chunks=(1, 5)) + elevation = _make_raster(elev, backend, chunks=(1, 5)) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + result = _compute(surface_distance(raster, elevation)) + + assert result[0, 0] == 0.0 + assert result[0, 1] == pytest.approx(1.0, abs=1e-5) + # The Inf cell and everything behind it are unreachable. + assert np.all(np.isnan(result[0, 2:])) + + +@pytest.mark.parametrize("backend", ALL_BACKENDS) +def test_all_nan_elevation(backend): + """An entirely impassable elevation surface reaches nothing.""" + _needs(backend) + source = np.zeros((4, 4), dtype=np.float64) + source[1, 1] = 1.0 + elev = np.full((4, 4), np.nan, dtype=np.float64) + + raster = _make_raster(source, backend, chunks=(2, 2)) + elevation = _make_raster(elev, backend, chunks=(2, 2)) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + result = _compute(surface_distance(raster, elevation)) + + assert np.all(np.isnan(result)) + + +# --------------------------------------------------------------------------- +# Geometric degeneracies +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize("backend", ALL_BACKENDS) +def test_single_pixel_raster(backend): + """A 1x1 raster: the pixel is either its own source or unreachable.""" + _needs(backend) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + hit = _compute(surface_distance( + _make_raster(np.array([[1.0]]), backend, chunks=(1, 1)), + _make_raster(np.array([[0.0]]), backend, chunks=(1, 1)), + )) + miss = _compute(surface_distance( + _make_raster(np.array([[0.0]]), backend, chunks=(1, 1)), + _make_raster(np.array([[0.0]]), backend, chunks=(1, 1)), + )) + + assert hit.shape == (1, 1) + assert hit[0, 0] == 0.0 + assert np.isnan(miss[0, 0]) + + +@pytest.mark.parametrize("backend", ALL_BACKENDS) +def test_column_strip(backend): + """An Nx1 column exercises the kernel with no horizontal neighbours.""" + _needs(backend) + source = np.zeros((5, 1), dtype=np.float64) + source[0, 0] = 1.0 + elev = np.zeros((5, 1), dtype=np.float64) + + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + result = _compute(surface_distance( + _make_raster(source, backend, chunks=(2, 1)), + _make_raster(elev, backend, chunks=(2, 1)), + )) + + np.testing.assert_allclose(result.ravel(), [0.0, 1.0, 2.0, 3.0, 4.0], + rtol=1e-5) + + +# --------------------------------------------------------------------------- +# Parameter coverage +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize("backend", ALL_BACKENDS) +def test_connectivity_4_all_backends(backend): + """4-connectivity is honoured by every backend, not just eager numpy.""" + _needs(backend) + source = np.zeros((6, 6), dtype=np.float64) + source[0, 0] = 1.0 + elev = np.zeros((6, 6), dtype=np.float64) + + raster = _make_raster(source, backend, chunks=(3, 3)) + elevation = _make_raster(elev, backend, chunks=(3, 3)) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + conn4 = _compute(surface_distance(raster, elevation, connectivity=4)) + conn8 = _compute(surface_distance(raster, elevation, connectivity=8)) + + # Without diagonals the corner pixel costs two cardinal steps. + assert conn4[1, 1] == pytest.approx(2.0, abs=1e-5) + assert conn8[1, 1] == pytest.approx(np.sqrt(2), abs=1e-5) + # 4-connectivity can never beat 8-connectivity. + assert np.all(conn4 >= conn8 - 1e-5) + + +@pytest.mark.parametrize("backend", ['cupy', 'dask+numpy', 'dask+cupy']) +def test_geodesic_unsupported_backends_raise(backend): + """Geodesic mode is numpy-only and says so per backend.""" + _needs(backend) + raster = _make_raster(np.zeros((4, 4)), backend, chunks=(2, 2)) + elevation = _make_raster(np.zeros((4, 4)), backend, chunks=(2, 2)) + + with pytest.raises(NotImplementedError, match="geodesic"): + surface_distance(raster, elevation, method='geodesic') + + +def test_invalid_dims(): + """A raster whose dims do not match x/y is rejected by name.""" + source = _make_raster(np.zeros((3, 3))).rename({'y': 'x', 'x': 'y'}) + elevation = _make_raster(np.zeros((3, 3))) + with pytest.raises(ValueError, match=r"raster.dims should be"): + surface_distance(source, elevation) + + +# --------------------------------------------------------------------------- +# Dataset support +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize( + "func", [surface_distance, surface_allocation, surface_direction], + ids=['distance', 'allocation', 'direction'], +) +def test_dataset_input(func): + """@supports_dataset maps the function over every data variable.""" + source = np.zeros((1, 4), dtype=np.float64) + source[0, 0] = 1.0 + elevation = _make_raster(np.zeros((1, 4))) + + ds = xr.Dataset( + {'a': _make_raster(source), 'b': _make_raster(source * 3.0)}, + attrs={'source': 'test'}, + ) + result = func(ds, elevation) + + assert isinstance(result, xr.Dataset) + assert list(result.data_vars) == ['a', 'b'] + assert result.attrs == {'source': 'test'} + for name in ('a', 'b'): + expected = func(ds[name], elevation) + np.testing.assert_allclose(_compute(result[name]), + _compute(expected), equal_nan=True) + + +def test_dataset_allocation_keeps_per_variable_values(): + """Allocation over a Dataset reports each variable's own source values.""" + source = np.zeros((1, 4), dtype=np.float64) + source[0, 0] = 1.0 + elevation = _make_raster(np.zeros((1, 4))) + + ds = xr.Dataset({'a': _make_raster(source), + 'b': _make_raster(source * 3.0)}) + result = surface_allocation(ds, elevation) + + np.testing.assert_allclose(_compute(result['a']), np.full((1, 4), 1.0)) + np.testing.assert_allclose(_compute(result['b']), np.full((1, 4), 3.0)) + + +# --------------------------------------------------------------------------- +# Metadata preservation +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize( + "func", [surface_distance, surface_allocation, surface_direction], + ids=['distance', 'allocation', 'direction'], +) +@pytest.mark.parametrize("backend", ALL_BACKENDS) +def test_attrs_and_coords_preserved(func, backend): + """Input attrs, coords and dims survive the round trip. + + The module reads ``res`` off the input attrs to build its edge costs, + so losing them downstream would break any chained call. + """ + _needs(backend) + source = np.zeros((6, 6), dtype=np.float64) + source[1, 1] = 1.0 + elev = np.zeros((6, 6), dtype=np.float64) + + raster = _make_raster(source, backend, chunks=(3, 3)) + raster.attrs['crs'] = 4326 + raster.attrs['units'] = 'meters' + elevation = _make_raster(elev, backend, chunks=(3, 3)) + + result = func(raster, elevation, max_distance=_BOUNDED) + + assert result.attrs == raster.attrs + assert result.dims == raster.dims + assert result.shape == raster.shape + np.testing.assert_array_equal(result.y.values, raster.y.values) + np.testing.assert_array_equal(result.x.values, raster.x.values) + + +@pytest.mark.parametrize( + "func", [surface_distance, surface_allocation, surface_direction], + ids=['distance', 'allocation', 'direction'], +) +def test_custom_dim_names_preserved(func): + """lat/lon dims are not silently renamed to y/x.""" + source = np.zeros((4, 4), dtype=np.float64) + source[0, 0] = 1.0 + raster = xr.DataArray( + source, dims=['lat', 'lon'], + coords={'lat': np.arange(4, dtype=np.float64), + 'lon': np.arange(4, dtype=np.float64)}, + attrs={'res': (1.0, 1.0)}, + ) + elevation = raster.copy(data=np.zeros((4, 4))) + + result = func(raster, elevation, x='lon', y='lat') + + assert result.dims == ('lat', 'lon') + np.testing.assert_array_equal(result.lat.values, raster.lat.values) From e4b2da3bbdcb64d74688a020a1f482d04a0ca307 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Sun, 16 Aug 2026 12:20:29 -0400 Subject: [PATCH 4/9] Record the surface_distance test-coverage sweep (#3725) --- .claude/sweep-test-coverage-state.csv | 1 + 1 file changed, 1 insertion(+) diff --git a/.claude/sweep-test-coverage-state.csv b/.claude/sweep-test-coverage-state.csv index c57ee8e5c..cef6da301 100644 --- a/.claude/sweep-test-coverage-state.csv +++ b/.claude/sweep-test-coverage-state.csv @@ -29,6 +29,7 @@ rasterize,2026-06-18,2614;3102;3105;3296;3383,HIGH,4,,"Pass 7 (2026-06-18, deep- reproject,2026-06-09,2618;3050;3100;3101;3141,MEDIUM,1,,"CI follow-up same day: first CI run of the threaded streaming branch hard-crashed macos-arm64 py3.14 (SIGABRT in numba call_cfunc, two ThreadPoolExecutor threads concurrently inside try_numba_transform/tmerc_inverse) -- the projection kernels are @njit(parallel=True) and numba's workqueue threading layer aborts on concurrent entry; filed source bug #3141. Test fix: threaded parity test now uses transform_precision=0 (per-thread pyproj Transformer, no numba), the NaN multi-tile test and 3-D xfail forced serial (max_memory=1) so the numba fast path stays covered without concurrent entry. windows-3.14 failure was fail-fast collateral (its suite fully passed). Pass 2026-06-09 (deep-sweep test-coverage): delta re-sweep one day after the 2026-06-08 pass; module modified today by #3077 (datum-probe warning silencing) and #3081 (merge output-size guard backend-aware) -- both landed WITH their own tests (TestDatumProbeNoProjWarning; TestSecurityGuards merge-guard trio incl. the monkeypatched in-memory raise), so the delta added no gap; the guard branching is is_dask-only, so cupy eager shares the tested numpy branch (no per-backend guard test needed). Found one MEDIUM Cat 1 gap every prior pass missed: the 5th dispatch branch of reproject() -- the streaming fallback (_reproject_streaming / _process_tile_batch / _parse_max_memory, taken when source >512MB and dask is not importable) -- had zero coverage anywhere; _parse_max_memory only runs on that branch so the existing max_memory kwarg tests never reached it. Filed #3101, added test_reproject_streaming_3101.py (15 tests: parity vs in-memory numpy for threaded / serial(max_memory=1) / single-tile / nearest+NaN, plus 10 _parse_max_memory unit cases). Probe surfaced source bug #3100: streaming assembly allocates a 2-D output buffer but 3-D sources yield (h,w,b) tiles -> ValueError broadcast in both assembly loops; pinned with strict xfail, source fix left to #3100 (test-only PR, source untouched). CPU-only path so no GPU tests needed (CUDA host; file ran 14 passed + 1 xfailed). LOW carried (documented, not fixed): reproject(name=) / merge(name=) override values untested (only merge name fallback covered); non-square-cellsize successful anisotropic run; dask.bag distributed branch of _reproject_streaming still unexercised (needs a live distributed client). || PREVIOUS: Pass 2026-06-08 (deep-sweep test-coverage): #3050 closes the one live gap found this pass. reproject()'s dask+cupy backend was parity-tested only with resampling='cubic' (TestCupyPyprojFallbackParity::test_projected_to_projected_dask_cupy_match); nearest/bilinear were covered on numpy (end-to-end) and eager cupy (parametrized test_projected_to_projected_numpy_cupy_match) but never on the dask+cupy chunk-assembly path. Parametrized that test over ['nearest','bilinear','cubic']; all 3 RUN+PASS on a CUDA host. Cat 4 MEDIUM (resampling-mode parameter coverage on the dask+cupy backend). Test-only, source untouched. Re-confirmed _merge.merge() has NO genuine cupy/dask+cupy backend (_merge_inmemory/_merge_dask use _merge_arrays_numpy + raster.values; _merge_arrays_cupy is imported but never dispatched = dead code, not a test gap) matching the prior pass's observation. reproject() otherwise saturated across all 4 backends, NaN/Inf/all-NaN, degenerate shapes, metadata, vertical, bounds_policy, integer nodata. LOW (documented, not filed): dask+cupy resampling-mode parity is the only per-mode-per-backend cell that had been missing. || PREVIOUS: Pass 2026-05-29: reproject already has a deep suite (369 tests in test_reproject.py + coverage/gate files) covering all 4 backends, NaN/Inf/all-NaN/all-Inf, 1x1/2x2, metadata, vertical shift, bounds_policy x backends, integer nodata x backends. Gaps found: Cat 3 HIGH single-row (1xN) and single-col (Nx1) strip rasters never tested (hit size<2 branch of _validate_regular_axis + degenerate resampling axis); Cat 3 MEDIUM constant-value/zero-gradient raster never reprojected. Added TestDegenerateShapeReproject (12 tests): 1xN+Nx1 strips x numpy/dask/cupy/dask+cupy, constant raster numpy value-preservation + cross-backend parity. All 12 executed and passed on a CUDA host. Test-only, no source change (#2618). LOW (documented only): _merge._merge_arrays_cupy imported but never called by merge() (host-bounces via _merge_arrays_numpy) - dead-code source observation not a test gap; non-square cellsize reproject only covered via resolution-tuple validation errors not a successful anisotropic run." resample,2026-05-29,2547;2615,HIGH,1;2;3;5,,"Pass 2 (2026-05-29): added test_resample_cupy_agg_fallback_2615.py (6 tests, all passing on CUDA host). Closes Cat 1 MEDIUM backend-coverage gap: the cupy eager aggregate CPU fallback for average/min/max at a NON-integer downsample factor (_run_cupy fy==int(fy) branch in resample.py ~L957-973) was never exercised; existing TestCuPyParity used 12x12 scale 0.5 (integer factor 2 -> GPU reshape path) and only median/mode hit the host fallback. New tests use 10x10 scale 0.3 (factor 3.33) for average/min/max parity vs numpy plus a NaN-masked variant. Issue #2615. Module is otherwise very thoroughly covered (test_resample.py + 3 supplementary files); no remaining HIGH gaps found. Pass 1 (2026-05-27): added test_resample_coverage_2026_05_27.py with 70 tests (68 passing, 2 skipped). Closes Cat 3 HIGH Nx1 single-column gap across numpy/cupy/dask+numpy/dask+cupy x 8 methods (nearest/bilinear/cubic/average/min/max/median/mode) plus Nx1 upsample-nearest parity and Nx1 cross-backend aggregate parity. Closes Cat 2 MEDIUM NaN-parity gap on cupy and dask+cupy (existing TestCuPyParity/TestDaskCuPyParity used random data without NaN; the weight-mask gate and spline-prepad had no GPU NaN coverage). Closes Cat 3 MEDIUM all-equal-value raster across 8 methods (downsample) and 3 interp methods (upsample) plus a constant-with-NaN aggregate variant. Closes Cat 5 MEDIUM non-default dim-name propagation: lat/lon, latitude/longitude, and (channel, lat, lon) 3D round-trip without being renamed to y/x; per-dim attrs (units) preserved. Closes Cat 3 MEDIUM empty-raster behaviour pin: 0-row and 0-col rasters raise (currently IndexError) -- contract covered. Filed source-bug issue #2547: cubic on dask backends fails for Nx1 / arrays smaller than depth=16; the 2 skipped tests in this file gate on that fix landing. Source untouched." slope,2026-05-29,2697,MEDIUM,3,,"PR #2703: added degenerate-shape tests (1x1/1xN/Nx1) for all 4 planar backends + geodesic; no live bug, pins all-NaN+shape contract. CUDA host: cupy/dask+cupy ran. Backend/NaN/param/metadata coverage already complete." +surface_distance,2026-08-16,3725,HIGH,1;2;3;4;5,87,"Deep-sweep 2026-08-16 test-coverage on a CUDA host (cupy + dask+cupy RAN, not skipped). Before: 37 tests, branch_cov 80 (570/695). After: 106 tests, branch_cov 87. Cat 1 HIGH: _surface_distance_dask_cupy (:1232-1289) had ZERO coverage (neither the bounded map_overlap branch nor the unbounded convert-to-dask+numpy branch), and surface_direction was only ever called on eager numpy -- no cupy, no dask, no dask+cupy. That second hole was hiding a real SOURCE bug (#3713, filed by the sibling documentation sweep, no PR then): _run_tile records GLOBAL src_row/src_col while _finalize_direction built its pixel grid from BLOCK-LOCAL np.arange(H)/np.arange(W), so the dask iterative path returned the top-left tile's bearing field tiled across the whole raster (flat 6x6 chunked (3,3): 360 deg max error, four pixels reporting 0 'you are the source'). Fixed in the same PR by threading row_offset/col_offset through _extract_output -> _finalize_direction from _assemble_sd; eager numpy and bounded map_overlap keep offset 0 because they seed local indices. Fail-before/pass-after verified. Also covered: cupy target_values seeding (:560-561) and cupy no-source early return (:571-573); dask iterative + target_values (:698-701, :956-959); connectivity=4 off the eager path (:919-932). Cat 2 MEDIUM: Inf elevation and all-NaN elevation untested -- both behave correctly (non-finite = barrier) on all 4 backends, coverage gap not a bug. Cat 3 MEDIUM: no 1x1 and no Nx1 column-strip tests (only 1xN rows existed); both correct on all 4 backends. Cat 4 MEDIUM: all three public funcs carry @supports_dataset and none was ever called with a Dataset (same gap as the proximity trio) -- works correctly; error paths untested (raster.dims rejection :1306, geodesic NotImplementedError per backend :1371/:1379/:1401). Cat 5 HIGH: nothing asserted attrs/coords/dims preservation even though _compute reads res off attrs for its edge costs -- all 4 backends preserve them correctly. LOW (documented, not fixed): non-square cells (res_x != res_y) untested but correct (verified N=2.0/E=1.0 on a 2x1 grid); _available_memory_bytes psutil + 2GB-default fallbacks (:107-115), _available_gpu_memory_bytes failure path (:129-130), _check_gpu_memory GPU-branch raise (:156,:159), _precompute_dd_grid memory guard (:423-426) all unexercised; _dijkstra_geodesic visited/over-budget/non-finite branches (:295,:299, :313) unreached. Pre-existing flake8 F841 at test file line 236 (sd assigned but unused in test_allocation_consistency) left alone per surgical-change rule. Post-fix branch_cov of 87 is still measured with NUMBA_DISABLE_JIT=1, under which the cupy tests fail, so true with-GPU coverage is higher; full file is 106/106 GREEN with JIT on." templates,2026-06-30,3580,MEDIUM,1;4,,"Deep-sweep test-coverage re-run on a CUDA host (cuda available). Module is already heavily tested (281 tests): 4-backend matrix (numpy/dask+numpy/cupy/dask+cupy) + dask alias + bad-backend all green; preserve area/shape across backends; single-pixel + Nx1/1xN strips; cell-cap and chunk-count guards; padding/tiling helpers; country/region/city resolution + aliases; CF metadata + no-pyproj fallback. Cat 2 N/A (procedural generator, no raster input). Found two MEDIUM parameter-coverage gaps, no source bug. Cat 4: chunks=tuple only exercised via internal _estimate_n_chunks, never end-to-end through from_template (int/'auto' were); the dask chunk-count guard message on the explicit height/width path was untested (only the resolution-path message and the eager cell-cap message named the knob). Added test-only test_chunks_tuple_through_public_api and test_explicit_shape_chunk_count_message_names_height_width; both RAN and PASSED on the CUDA host (281 passed). LOW (documented, no test): non-NaN fill is only asserted on eager numpy (fill=0); probed live and works on all 4 backends but cross-backend fill value parity is not asserted. PR #3580 opened with the two tests." viewshed,2026-05-29,2693,HIGH,1;2;5,,"Pass 1 (2026-05-29): added 4 new test groups to test_viewshed.py (13 new tests + 1 xfail, all passing/xfailing on a CUDA+RTX host). Closes Cat 1 HIGH backend-coverage gap: the dask+cupy dispatch path in _viewshed_dask (Tier B) and _viewshed_windowed (max_distance) was registered but never invoked by any test -- added test_viewshed_dask_cupy_flat (analytical-angle parity, atol 0.03) and test_viewshed_dask_cupy_max_distance (windowed GPU run; observer cell 180, corners INVISIBLE). Both use non-zero flat terrain (1.3) because the RTX mesh builder rejects an all-zero raster (#1378). Closes Cat 5 HIGH metadata-preservation gap: only the numpy test_viewshed called general_output_checks; the cupy/dask/dask+cupy and max_distance paths never asserted attrs/coords/dims/array-type preservation. Added parametrised test_viewshed_metadata_preserved over {numpy,cupy,dask+numpy,dask+cupy} x {full, max_distance=2.0}: asserts attrs==, dims==, shape==, x/y coords allclose; runs general_output_checks (full type parity) for all backends except dask+cupy. Closes Cat 2 HIGH NaN-input gap and surfaced source bug #2693: viewshed on a numpy raster crashes with ValueError 'node not found' from _delete_from_tree when a NaN cell sits at certain positions (e.g. (2,4) in a 5x5 with observer at (2,2)), while NaN at (1,1)/(0,0)/(4,4) runs fine. Added test_viewshed_nan_input_supported_positions (parametrised working positions, asserts observer=180 and NaN cell is INVISIBLE/NaN) plus test_viewshed_nan_input_crashing_position (xfail strict, raises, links #2693). Noted but NOT fixed (source change out of scope for test sweep): the dask+cupy backend does not preserve the cupy backing -- _viewshed_dask computes then rewraps via da.from_array(result_np), so the output computes to numpy not cupy; general_output_checks is skipped for dask+cupy for that reason (candidate for the metadata/backend-parity sweep). LOW (documented only): non-square cell sizes; 1x1 and 1xN geometry covered behaviourally by probing (run without error). Test-only PR; viewshed.py untouched." visibility,2026-06-10,3192,HIGH,1;2;4,,"cupy cumulative_viewshed/visibility_frequency broken (numpy count + cupy viewshed) -> issue #3192 (dup #3193), fix in flight in #3205 with its own cupy parity tests, xfail pins dropped to avoid an XPASS race; added cupy _extract_transect+line_of_sight parity, NaN LOS, Fresnel-blocked branch; dask+metadata already covered" From fbfd64988d99e95e1242108a1ddc28bbc145cd9a Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Sun, 16 Aug 2026 12:20:37 -0400 Subject: [PATCH 5/9] Grow the surface_distance Dijkstra heap on demand (#3723) _dijkstra and _dijkstra_geodesic sized their lazy-deletion min-heap at height * width. A pixel is re-pushed every time its tentative distance improves, so live occupancy is bounded by the number of improving relaxations, not by the pixel count. On a 48x48 raster with 40% target pixels over 200 m of relief, peak occupancy reached 1.09x the cap and _heap_push wrote past the end of h_keys / h_rows / h_cols, aborting the interpreter with "double free or corruption". The heap now starts at one entry per pixel (which the seeding loop can never exceed) and doubles when full. cost_distance fixed the same defect with a static height * width * (n_neighbors + 1) bound, but that would multiply this module's heap allocation by 9 and force _BYTES_PER_PIXEL from 80 to about 272, tightening the memory guard on every caller. Adds regression tests for the numpy, dask iterative and geodesic paths, plus a SurfaceDistanceDenseTargets asv benchmark. The existing benchmark misses this regime: its integer source raster makes nearly every pixel a target, so almost no relaxation improves a distance. --- .claude/sweep-performance-state.csv | 2 +- benchmarks/benchmarks/surface_distance.py | 35 +++++++ xrspatial/surface_distance.py | 35 ++++++- xrspatial/tests/test_surface_distance.py | 121 ++++++++++++++++++++++ 4 files changed, 190 insertions(+), 3 deletions(-) diff --git a/.claude/sweep-performance-state.csv b/.claude/sweep-performance-state.csv index 77efd2a00..051a49ecb 100644 --- a/.claude/sweep-performance-state.csv +++ b/.claude/sweep-performance-state.csv @@ -43,7 +43,7 @@ resample,2026-04-15T12:00:00Z,SAFE,compute-bound,0,false-positive,Downgraded. GP sieve,2026-04-14T12:00:00Z,WILL OOM,memory-bound,0,false-positive,False positive. Memory guards already in place on both dask paths. CCL is inherently global — documented limitation. CuPy CPU fallback is deliberate and documented. sky_view_factor,2026-03-31T18:00:00Z,SAFE,compute-bound,0,, slope,2026-03-31T18:00:00Z,SAFE,compute-bound,0,, -surface_distance,2026-03-31T18:00:00Z,SAFE,memory-bound,0,1128,Memory guard added to dd_grid allocation. +surface_distance,2026-08-16,RISKY,compute-bound,0,3723,"CRITICAL #3723/PR: _dijkstra and _dijkstra_geodesic sized the lazy-deletion heap height*width; peak occupancy hit 1.09x cap at 40% target density, _heap_push wrote OOB (SIGABRT 'double free or corruption'); fixed by doubling on demand rather than cost_distance's static (n_neighbors+1) bound, which would have raised the 80 B/px guard to ~272. MEDIUM (documented, not fixed): dask iterative path recomputes each source and elev block 6x (2 in _preprocess_tiles_sd, 4 across sweeps, 1 in _assemble_sd) - measured on 40x40/80x80/120x120; MEDIUM: cupy Bellman-Ford relaxation is 6-27x slower than the numpy Dijkstra (512^2: 0.315s vs 0.051s, 373 kernel launches each with a device sync), and DIRECTION mode round-trips to host for _vectorized_calc_direction; MEDIUM: map_overlap depth guard tests pad < max(chunks), so depth may approach chunk size (~9x redundant work) and exceed the smallest chunk when chunking is uneven. Memory guard is correctly chunk-scoped on dask (test_dask_path_bounded_per_chunk) - the morphology #3401 false-MemoryError hazard is NOT present. _sd_relax_kernel has ~8 float64 locals, no register pressure. OOM verdict RISKY: memory stays chunk-bounded but the whole iterative Dijkstra runs eagerly and single-threaded on the client at graph-build time." terrain,2026-03-31T18:00:00Z,RISKY,compute-bound,0,, terrain_metrics,2026-03-31T18:00:00Z,SAFE,memory-bound,0,, viewshed,2026-04-05T12:00:00Z,SAFE,memory-bound,0,fixed-in-tree,Tier B memory estimate tightened from 280 to 368 bytes/pixel (accounts for lexsort double-alloc + computed raster). astype copy=False avoids needless float64 copy. diff --git a/benchmarks/benchmarks/surface_distance.py b/benchmarks/benchmarks/surface_distance.py index 1cf7aad60..bdf1602b3 100644 --- a/benchmarks/benchmarks/surface_distance.py +++ b/benchmarks/benchmarks/surface_distance.py @@ -1,4 +1,5 @@ import numpy as np +import xarray as xr from xrspatial.surface_distance import ( surface_distance, surface_allocation, surface_direction, @@ -24,3 +25,37 @@ def time_surface_allocation(self, nx, type): def time_surface_direction(self, nx, type): surface_direction(self.agg, self.elev) + + +class SurfaceDistanceDenseTargets: + """Dijkstra with a mid-density target raster over rugged relief. + + This is the regime where the lazy-deletion heap holds more than one + live entry per pixel (#3723). The existing SurfaceDistance benchmark + misses it: its integer source raster makes nearly every pixel a + target, so almost no relaxation ever improves a distance. + """ + + params = ([200, 400], [0.05, 0.2, 0.4]) + param_names = ("nx", "target_fraction") + + def setup(self, nx, target_fraction): + ny = nx // 2 + rng = np.random.default_rng(71942) + source = np.zeros((ny, nx), dtype=np.float64) + n_targets = max(1, int(ny * nx * target_fraction)) + source.flat[rng.choice(ny * nx, size=n_targets, replace=False)] = 1.0 + elev = rng.random((ny, nx)) * 200.0 + + coords = dict(y=np.arange(ny, dtype=np.float64), + x=np.arange(nx, dtype=np.float64)) + self.agg = xr.DataArray(source, coords=coords, dims=["y", "x"], + attrs={"res": (1.0, 1.0)}) + self.elev = xr.DataArray(elev, coords=coords, dims=["y", "x"], + attrs={"res": (1.0, 1.0)}) + + def time_surface_distance(self, nx, target_fraction): + surface_distance(self.agg, self.elev) + + def peakmem_surface_distance(self, nx, target_fraction): + surface_distance(self.agg, self.elev) diff --git a/xrspatial/surface_distance.py b/xrspatial/surface_distance.py index c91aeb13e..81b9aa98e 100644 --- a/xrspatial/surface_distance.py +++ b/xrspatial/surface_distance.py @@ -170,6 +170,28 @@ def _check_gpu_memory(rows, cols): # --------------------------------------------------------------------------- +@ngjit +def _heap_grow(keys, rows, cols, size): + """Return copies of the heap arrays with twice the capacity. + + The Dijkstra kernels below use a lazy-deletion min-heap: a pixel is + pushed again every time its tentative distance improves, and stale + entries linger until a pop skips them. Live occupancy is therefore + bounded by the number of improving relaxations, not by the pixel + count, and a fixed ``height * width`` heap can overflow (#3723). + Growing on demand keeps the usual footprint at one entry per pixel + without capping the number of pushes. + """ + cap = 2 * len(keys) + new_keys = np.empty(cap, dtype=np.float64) + new_rows = np.empty(cap, dtype=np.int64) + new_cols = np.empty(cap, dtype=np.int64) + new_keys[:size] = keys[:size] + new_rows[:size] = rows[:size] + new_cols[:size] = cols[:size] + return new_keys, new_rows, new_cols + + @ngjit def _seed_sources(source_data, elev_data, target_values, dist, alloc, src_row, src_col): @@ -211,7 +233,9 @@ def _dijkstra(elev_data, height, width, max_distance, """ n_neighbors = len(dy) - max_heap = height * width + # Starting capacity: one entry per pixel, which the seeding loop below + # can never exceed. Relaxations can, so _heap_grow doubles it there. + max_heap = max(height * width, 1) h_keys = np.empty(max_heap, dtype=np.float64) h_rows = np.empty(max_heap, dtype=np.int64) h_cols = np.empty(max_heap, dtype=np.int64) @@ -260,6 +284,9 @@ def _dijkstra(elev_data, height, width, max_distance, alloc[vr, vc] = alloc[ur, uc] src_row[vr, vc] = src_row[ur, uc] src_col[vr, vc] = src_col[ur, uc] + if h_size == len(h_keys): + h_keys, h_rows, h_cols = _heap_grow( + h_keys, h_rows, h_cols, h_size) h_size = _heap_push(h_keys, h_rows, h_cols, h_size, new_cost, vr, vc) @@ -274,7 +301,8 @@ def _dijkstra_geodesic(elev_data, height, width, max_distance, """ n_neighbors = len(dy) - max_heap = height * width + # See _dijkstra for the heap-capacity rationale. + max_heap = max(height * width, 1) h_keys = np.empty(max_heap, dtype=np.float64) h_rows = np.empty(max_heap, dtype=np.int64) h_cols = np.empty(max_heap, dtype=np.int64) @@ -322,6 +350,9 @@ def _dijkstra_geodesic(elev_data, height, width, max_distance, alloc[vr, vc] = alloc[ur, uc] src_row[vr, vc] = src_row[ur, uc] src_col[vr, vc] = src_col[ur, uc] + if h_size == len(h_keys): + h_keys, h_rows, h_cols = _heap_grow( + h_keys, h_rows, h_cols, h_size) h_size = _heap_push(h_keys, h_rows, h_cols, h_size, new_cost, vr, vc) diff --git a/xrspatial/tests/test_surface_distance.py b/xrspatial/tests/test_surface_distance.py index 74ae62fbd..9bc6d7d35 100644 --- a/xrspatial/tests/test_surface_distance.py +++ b/xrspatial/tests/test_surface_distance.py @@ -792,3 +792,124 @@ def test_error_message_mentions_grid_size(self): surface_distance(raster, elevation) with pytest.raises(MemoryError, match="dask"): surface_distance(raster, elevation) + + +# --------------------------------------------------------------------------- +# Heap capacity regression (#3723) +# --------------------------------------------------------------------------- + + +def _reference_surface_distance(source, elev, connectivity=8, cellsize=1.0): + """Pure-Python multi-source Dijkstra reference for surface distance.""" + import heapq + + h, w = source.shape + diag = np.sqrt(2.0) * cellsize + if connectivity == 8: + nbrs = [(-1, -1, diag), (-1, 0, cellsize), (-1, 1, diag), + (0, -1, cellsize), (0, 1, cellsize), + (1, -1, diag), (1, 0, cellsize), (1, 1, diag)] + else: + nbrs = [(0, -1, cellsize), (-1, 0, cellsize), + (1, 0, cellsize), (0, 1, cellsize)] + + dist = np.full((h, w), np.inf) + heap = [] + for r in range(h): + for c in range(w): + if (source[r, c] != 0 and np.isfinite(source[r, c]) + and np.isfinite(elev[r, c])): + dist[r, c] = 0.0 + heapq.heappush(heap, (0.0, r, c)) + + done = np.zeros((h, w), dtype=bool) + while heap: + d, r, c = heapq.heappop(heap) + if done[r, c]: + continue + done[r, c] = True + for dr, dc, hd in nbrs: + vr, vc = r + dr, c + dc + if not (0 <= vr < h and 0 <= vc < w) or done[vr, vc]: + continue + if not np.isfinite(elev[vr, vc]): + continue + dz = elev[vr, vc] - elev[r, c] + nd = d + np.sqrt(hd * hd + dz * dz) + if nd < dist[vr, vc]: + dist[vr, vc] = nd + heapq.heappush(heap, (nd, vr, vc)) + return dist + + +def _dense_target_scene(n=48, n_targets=921, relief=200.0, seed=1): + """Dense-target scene that overflowed the old height*width heap.""" + rng = np.random.default_rng(seed) + source = np.zeros((n, n), dtype=np.float64) + source.flat[rng.choice(n * n, size=n_targets, replace=False)] = 1.0 + elev = rng.random((n, n)) * relief + return source, elev + + +def test_dense_targets_do_not_overflow_the_heap(): + """A lazy-deletion heap can exceed height*width live entries. + + Before #3723 the heap arrays were sized height*width, so this scene + made _heap_push write past the end of them (SIGABRT without bounds + checking, IndexError with NUMBA_BOUNDSCHECK=1). + """ + source, elev = _dense_target_scene() + raster = _make_raster(source) + elevation = _make_raster(elev) + + result = _compute(surface_distance(raster, elevation, connectivity=8)) + expected = _reference_surface_distance(source, elev, connectivity=8) + + assert np.all(np.isfinite(result)) + np.testing.assert_allclose(result, expected.astype(np.float32), + rtol=1e-5, atol=1e-4) + + +def test_dense_targets_allocation_and_direction_do_not_overflow(): + """The allocation and direction modes share the same Dijkstra kernel.""" + source, elev = _dense_target_scene() + raster = _make_raster(source) + elevation = _make_raster(elev) + + alloc = _compute(surface_allocation(raster, elevation, connectivity=8)) + direction = _compute(surface_direction(raster, elevation, connectivity=8)) + + assert np.all(alloc == 1.0) + assert np.all(np.isfinite(direction)) + + +@pytest.mark.skipif(da is None, reason="dask not installed") +def test_dense_targets_dask_iterative_does_not_overflow(): + """The dask iterative path runs the same kernel per tile.""" + source, elev = _dense_target_scene() + raster_np = _make_raster(source) + elev_np = _make_raster(elev) + raster = _make_raster(source, backend='dask+numpy', chunks=(48, 48)) + elevation = _make_raster(elev, backend='dask+numpy', chunks=(48, 48)) + + np_result = _compute(surface_distance(raster_np, elev_np)) + with pytest.warns(UserWarning, match="iterative"): + dask_result = _compute(surface_distance(raster, elevation)) + + np.testing.assert_allclose(dask_result, np_result, rtol=1e-5, + equal_nan=True) + + +def test_geodesic_dense_targets_do_not_overflow(): + """_dijkstra_geodesic carries the same heap sizing.""" + source, elev = _dense_target_scene(n=32, n_targets=410, relief=200.0) + h, w = source.shape + coords = {'y': np.linspace(10.0, 10.0 + 0.01 * (h - 1), h), + 'x': np.linspace(20.0, 20.0 + 0.01 * (w - 1), w)} + raster = xr.DataArray(source, dims=['y', 'x'], coords=coords, + attrs={'res': (0.01, 0.01)}) + elevation = xr.DataArray(elev, dims=['y', 'x'], coords=coords, + attrs={'res': (0.01, 0.01)}) + + result = _compute(surface_distance(raster, elevation, method='geodesic')) + assert np.all(np.isfinite(result)) From 17dcefcc3c3cf2763b291b5e7fa80fb92fc350c9 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Sun, 16 Aug 2026 12:21:03 -0400 Subject: [PATCH 6/9] Pin surface_distance against scipy csgraph and document the allocation tie-break (#3726) Validated surface_distance/allocation/direction against an independent reference: scipy.sparse.csgraph.dijkstra over an explicitly built grid graph with sqrt(horizontal^2 + dz^2) edge weights. All four backends reproduce it bit for bit. Adds golden-value tests pinning that parity, plus a tilted-plane case checked against the closed-form surface distance, since the module had no reference test of any kind. Also documents that surface_allocation and surface_direction resolve exact distance ties differently on each backend (#3726). The distances agree; only the reported source varies. --- xrspatial/surface_distance.py | 15 +++ xrspatial/tests/test_surface_distance.py | 146 +++++++++++++++++++++++ 2 files changed, 161 insertions(+) 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) From 0babfb464039660e47af0c30c5445862d9d4ac25 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Sun, 16 Aug 2026 12:21:05 -0400 Subject: [PATCH 7/9] Record the surface_distance reference-validation sweep result --- .claude/sweep-reference-validation-state.csv | 1 + 1 file changed, 1 insertion(+) 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)." From 4afb223dd1ffdc1328aa31bdd15050bdfb10963d Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Sun, 16 Aug 2026 12:24:49 -0400 Subject: [PATCH 8/9] Address review findings on the heap fix (#3723) - Guard the seeding loops in both kernels as well. The old comment argued the loop can push at most one entry per pixel into a heap sized for one entry per pixel, which is true today, but it leaves a memory-safety invariant resting on a comment. Every push site now checks capacity. - Say in the _BYTES_PER_PIXEL comment that the heap terms are a starting capacity, not a bound, so 80 reads as a floor rather than an estimate that quietly went stale. - Move the rationale for the capacity choice from _heap_grow into _dijkstra, where the capacity is picked. - Cover the bounded map_overlap dask branch with a dense-target test. It calls _surface_distance_numpy per padded chunk, so it was already fixed, but only the iterative branch had a test. --- xrspatial/surface_distance.py | 31 +++++++++++++++--------- xrspatial/tests/test_surface_distance.py | 18 ++++++++++++++ 2 files changed, 37 insertions(+), 12 deletions(-) diff --git a/xrspatial/surface_distance.py b/xrspatial/surface_distance.py index 81b9aa98e..e682eaf11 100644 --- a/xrspatial/surface_distance.py +++ b/xrspatial/surface_distance.py @@ -89,6 +89,12 @@ class cupy: # type: ignore[no-redef] # output (float32) 4 # direction-mode temps ~16 # Total ~80 bytes/pixel. A 50000x50000 raster needs ~200 GB. +# +# The heap terms above are the starting capacity of one entry per pixel. +# Dense-target rasters push more than that and _heap_grow doubles the +# three arrays, so the real figure can exceed 80 (48 bytes/pixel of heap +# after one doubling, and ~72 transiently while the copy is in flight). +# This estimate is a floor, not a bound. _BYTES_PER_PIXEL = 80 # CuPy backend skips the explicit binary heap (parallel relaxation instead) @@ -172,16 +178,7 @@ def _check_gpu_memory(rows, cols): @ngjit def _heap_grow(keys, rows, cols, size): - """Return copies of the heap arrays with twice the capacity. - - The Dijkstra kernels below use a lazy-deletion min-heap: a pixel is - pushed again every time its tentative distance improves, and stale - entries linger until a pop skips them. Live occupancy is therefore - bounded by the number of improving relaxations, not by the pixel - count, and a fixed ``height * width`` heap can overflow (#3723). - Growing on demand keeps the usual footprint at one entry per pixel - without capping the number of pushes. - """ + """Return copies of the heap arrays with twice the capacity.""" cap = 2 * len(keys) new_keys = np.empty(cap, dtype=np.float64) new_rows = np.empty(cap, dtype=np.int64) @@ -233,8 +230,12 @@ def _dijkstra(elev_data, height, width, max_distance, """ n_neighbors = len(dy) - # Starting capacity: one entry per pixel, which the seeding loop below - # can never exceed. Relaxations can, so _heap_grow doubles it there. + # This is a lazy-deletion min-heap: a pixel is pushed again every time + # its tentative distance improves, and stale entries linger until a pop + # skips them. Live occupancy therefore tracks the number of improving + # relaxations, not the pixel count, so the old fixed height*width heap + # let _heap_push write past the end of these arrays (#3723). Start at + # one entry per pixel and double whenever a push would overflow. max_heap = max(height * width, 1) h_keys = np.empty(max_heap, dtype=np.float64) h_rows = np.empty(max_heap, dtype=np.int64) @@ -247,6 +248,9 @@ def _dijkstra(elev_data, height, width, max_distance, for r in range(height): for c in range(width): if dist[r, c] < np.inf: + if h_size == len(h_keys): + h_keys, h_rows, h_cols = _heap_grow( + h_keys, h_rows, h_cols, h_size) h_size = _heap_push(h_keys, h_rows, h_cols, h_size, dist[r, c], r, c) @@ -313,6 +317,9 @@ def _dijkstra_geodesic(elev_data, height, width, max_distance, for r in range(height): for c in range(width): if dist[r, c] < np.inf: + if h_size == len(h_keys): + h_keys, h_rows, h_cols = _heap_grow( + h_keys, h_rows, h_cols, h_size) h_size = _heap_push(h_keys, h_rows, h_cols, h_size, dist[r, c], r, c) diff --git a/xrspatial/tests/test_surface_distance.py b/xrspatial/tests/test_surface_distance.py index 9bc6d7d35..ffdf3e439 100644 --- a/xrspatial/tests/test_surface_distance.py +++ b/xrspatial/tests/test_surface_distance.py @@ -900,6 +900,24 @@ def test_dense_targets_dask_iterative_does_not_overflow(): equal_nan=True) +@pytest.mark.skipif(da is None, reason="dask not installed") +def test_dense_targets_dask_bounded_does_not_overflow(): + """The bounded map_overlap path runs the kernel per padded chunk.""" + source, elev = _dense_target_scene() + raster_np = _make_raster(source) + elev_np = _make_raster(elev) + raster = _make_raster(source, backend='dask+numpy', chunks=(16, 16)) + elevation = _make_raster(elev, backend='dask+numpy', chunks=(16, 16)) + + np_result = _compute(surface_distance(raster_np, elev_np, + max_distance=10.0)) + dask_result = _compute(surface_distance(raster, elevation, + max_distance=10.0)) + + np.testing.assert_allclose(dask_result, np_result, rtol=1e-5, + equal_nan=True) + + def test_geodesic_dense_targets_do_not_overflow(): """_dijkstra_geodesic carries the same heap sizing.""" source, elev = _dense_target_scene(n=32, n_targets=410, relief=200.0) From dd54912ca7a6f4f61ea477dd9f81364869bfe2e4 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Fri, 6 Nov 2026 10:17:58 -0500 Subject: [PATCH 9/9] Pin the surface_direction docstring example to the north-up bearings from #3731 --- xrspatial/surface_distance.py | 9 +++++---- xrspatial/tests/test_surface_distance.py | 4 ++-- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/xrspatial/surface_distance.py b/xrspatial/surface_distance.py index 809805be6..1d7e3e0e6 100644 --- a/xrspatial/surface_distance.py +++ b/xrspatial/surface_distance.py @@ -1771,12 +1771,13 @@ def surface_direction( >>> elev['x'] = np.arange(m) >>> surface_direction(raster, elev).values - array([[135., 180., 225.], + array([[ 45., 360., 315.], [ 90., 0., 270.], - [ 45., 360., 315.]], dtype=float32) + [135., 180., 225.]], dtype=float32) - The pixel east of the source reads 270 (west) because the bearing - points back at the source, and the source itself reads 0. + The raster is north-up (``y`` descends with the row index), so the + row above the source reads 360 and the row below reads 180, the same + convention :func:`xrspatial.direction` uses. The source itself reads 0. """ 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 75a8e8efa..bab6bb91c 100644 --- a/xrspatial/tests/test_surface_distance.py +++ b/xrspatial/tests/test_surface_distance.py @@ -1225,9 +1225,9 @@ def test_direction_docstring_example_matches_output(): ]) raster, elev = _docstring_example_rasters(source, np.zeros((3, 3))) - expected = np.array([[135., 180., 225.], + expected = np.array([[45., 360., 315.], [90., 0., 270.], - [45., 360., 315.]], dtype=np.float32) + [135., 180., 225.]], dtype=np.float32) result = surface_direction(raster, elev) assert result.dtype == np.float32 np.testing.assert_allclose(result.values, expected, rtol=1e-6)