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