Document surface_distance backends, geodesic limits, and add runnable examples

.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

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,

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)