Benchmark the surface_distance paths that were never timed (#3709) (#3717)

The benchmark file covered one of six compute paths, and the path it did
cover ran on a source raster where 99.95% of pixels were sources, so the
Dijkstra relaxation body never executed.

- parameterize cupy and dask+cupy alongside numpy and dask, matching
  pathfinding.py
- replace the near-uniform integer source raster with scattered point
  sources so the frontier crosses the grid
- add time_surface_distance_bounded with a finite max_distance whose
  pixel radius stays inside one chunk, which routes the dask backends
  through map_overlap instead of the iterative tile fallback
- call .compute() on dask results, following the convention in flood.py,
  twi.py and interpolate.py
- add a numpy-only geodesic class for _dijkstra_geodesic and
  _precompute_dd_grid

Benchmark file only, no source changes. All 50 parameter combinations
were run on this host, 13s for a single pass.

Author: brendancol

.claude/sweep-benchmarks-state.csv

@@ -2,3 +2,4 @@ module,last_inspected,issue,severity_max,categories_found,notes
 edge_detection,2026-07-18,3672,MEDIUM,1,"No bench file existed; all 5 public funcs (sobel_x/y, prewitt_x/y, laplacian) uncovered. Compute delegates to convolve_2d (directly benchmarked, but only at 5x5/25x25 kernels, never 3x3, never via the wrappers) so MEDIUM not HIGH. Added benchmarks/benchmarks/edge_detection.py: EdgeDetection class, nx in [300,3000], numpy/cupy/dask, one timing method per public func. All 30 combos executed locally incl. cupy (GPU host). Cat 2/3/4 N/A (module implements no backends itself; no pre-existing bench to be broken)."
 geotiff,2026-07-02,3603,HIGH,1;2,"No benchmark existed for geotiff; open_geotiff/to_geotiff had zero asv coverage across numpy/dask/cupy. Added benchmarks/benchmarks/geotiff.py: WriteGeoTIFF (numpy/dask/cupy streaming), WriteCOG (numpy/cupy overview pyramid), ReadGeoTIFF (numpy/cupy decode), ReadGeoTIFFChunked (dask). All classes executed locally via direct call; cupy paths run on this GPU host. asv check discover fails suite-wide from an asv_runner + py3.14 metadata bug, unrelated to this file."
 pathfinding,2026-07-08,3645,HIGH,1;2;3,"Bench covered only numpy a_star_search at nx<=300; module also ships dask (separate sparse-Python A* + LRU chunk cache), cupy fallback, and public multi_stop_search with zero coverage. Extended AStarSearch to numpy/cupy/dask with nx up to 1000 (dask capped at 300, ~4s/call at 1000) and added MultiStopSearch (ordered + optimize_order). All combos executed locally incl. cupy (GPU host). LOW noted, not fixed: open-grid no-barrier/no-friction input is A* best case. dask+cupy not parameterized anywhere in suite (common.get_xr_dataarray has no such type). Existing bench imports/runs fine (Cat 4 clean)."
+surface_distance,2026-08-16,3709,HIGH,1;2;3,"Bench (26 lines) covered 1 of 6 compute paths. Cat 2 HIGH: cupy (_sd_relax_kernel) and dask+cupy (_surface_distance_dask_cupy) never parameterized; both run fine on this GPU host. Cat 2 HIGH: default max_distance=inf routed every dask call to the iterative tile fallback, so the bounded map_overlap branch the docstring recommends was never timed, and no .compute() was called (worked only because _sd_dask_iterative computes eagerly). Cat 3 MEDIUM: get_xr_dataarray(is_int=True) made 499769/500000 pixels sources at nx=1000, so the Dijkstra relaxation body never ran; max distance 3.6 vs 160 with sparse sources, 6.4x timing gap. Cat 1 MEDIUM: method='geodesic' (_dijkstra_geodesic + _precompute_dd_grid) unbenchmarked, ~2.2x planar cost. Cat 4 clean: file imported and ran. Fixed: 4 backends, sparse point sources, bounded method, .compute(), numpy-only geodesic class. All 50 combos executed locally, 13s single pass. LOW not fixed: nx=100 gives a 50x100 grid. Note: asv discover is broken suite-wide in this conda env (asv_runner dist-metadata bug on py3.14), so verification was by direct class invocation."

benchmarks/benchmarks/surface_distance.py

@@ -1,26 +1,105 @@
 import numpy as np
+import xarray as xr
 
 from xrspatial.surface_distance import (
     surface_distance, surface_allocation, surface_direction,
 )
+from xrspatial.utils import has_cuda_and_cupy
 
 from .common import get_xr_dataarray
 
 
+def _sparse_source_raster(ny, nx, type):
+    """Source raster with a handful of scattered target pixels.
+
+    ``get_xr_dataarray(is_int=True)`` draws integers over ``[-nx, nx)``,
+    and surface_distance treats every non-zero finite pixel as a source,
+    so all but roughly 1 in ``2 * nx`` pixels seed the search at distance
+    zero.  The Dijkstra relaxation body then never runs and the benchmark
+    times a heap drain instead of distance propagation.  Scattered point
+    sources make the frontier cross the whole grid, which is what these
+    functions are for.
+    """
+    rng = np.random.default_rng(71942)
+    z = np.zeros((ny, nx), dtype=np.float32)
+    n_sources = max(4, (ny * nx) // 20000)
+    rows = rng.integers(0, ny, n_sources)
+    cols = rng.integers(0, nx, n_sources)
+    z[rows, cols] = np.arange(1, n_sources + 1, dtype=np.float32)
+
+    chunks = (max(1, ny // 2), max(1, nx // 2))
+    if type == "cupy":
+        if not has_cuda_and_cupy():
+            raise NotImplementedError()
+        import cupy
+        z = cupy.asarray(z)
+    elif type == "dask":
+        import dask.array as da
+        z = da.from_array(z, chunks=chunks)
+    elif type == "dask+cupy":
+        if not has_cuda_and_cupy():
+            raise NotImplementedError()
+        import cupy
+        import dask.array as da
+        z = da.from_array(cupy.asarray(z), chunks=chunks)
+    elif type != "numpy":
+        raise RuntimeError(f"Unrecognised type {type}")
+
+    y = np.linspace(-90, 90, ny)
+    x = np.linspace(-180, 180, nx)
+    return xr.DataArray(z, coords=dict(y=y, x=x), dims=["y", "x"])
+
+
+def _compute(result):
+    if hasattr(result.data, "compute"):
+        result.data.compute()
+
+
 class SurfaceDistance:
-    params = ([100, 300, 1000], ["numpy", "dask"])
+    params = ([100, 300, 1000], ["numpy", "cupy", "dask", "dask+cupy"])
     param_names = ("nx", "type")
 
     def setup(self, nx, type):
         ny = nx // 2
-        self.agg = get_xr_dataarray((ny, nx), type, is_int=True)
+        self.agg = _sparse_source_raster(ny, nx, type)
         self.elev = get_xr_dataarray((ny, nx), type)
 
+        # A finite max_distance whose pixel radius stays inside one chunk
+        # (chunks are ny//2 x nx//2) routes the dask backends through the
+        # bounded map_overlap branch instead of the iterative tile one.
+        cellsize = min(360.0 / (nx - 1), 180.0 / (ny - 1))
+        self.max_distance = 20 * cellsize
+
     def time_surface_distance(self, nx, type):
-        surface_distance(self.agg, self.elev)
+        _compute(surface_distance(self.agg, self.elev))
+
+    def time_surface_distance_bounded(self, nx, type):
+        _compute(surface_distance(
+            self.agg, self.elev, max_distance=self.max_distance))
 
     def time_surface_allocation(self, nx, type):
-        surface_allocation(self.agg, self.elev)
+        _compute(surface_allocation(self.agg, self.elev))
 
     def time_surface_direction(self, nx, type):
-        surface_direction(self.agg, self.elev)
+        _compute(surface_direction(self.agg, self.elev))
+
+
+class SurfaceDistanceGeodesic:
+    """Great-circle horizontal distances from lat/lon coordinates.
+
+    Runs a separate numba kernel (``_dijkstra_geodesic``) behind a
+    precomputed per-pixel neighbour-distance grid, and costs about twice
+    the planar path.  numpy only: the module raises NotImplementedError
+    for geodesic on cupy, dask, and dask+cupy.
+    """
+
+    params = [300, 1000]
+    param_names = ("nx",)
+
+    def setup(self, nx):
+        ny = nx // 2
+        self.agg = _sparse_source_raster(ny, nx, "numpy")
+        self.elev = get_xr_dataarray((ny, nx), "numpy")
+
+    def time_surface_distance_geodesic(self, nx):
+        surface_distance(self.agg, self.elev, method="geodesic")