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

[brendancol]

Closes #3709.

benchmarks/benchmarks/surface_distance.py was 26 lines against a 1561-line module. It timed one of six compute paths, and it timed that one on input where the algorithm barely runs.

What changed

Benchmark file only. No source changes.

cupy and dask+cupy are now parameterized. _surface_distance_cupy (with the @cuda.jit _sd_relax_kernel) and _surface_distance_dask_cupy had no coverage at all. params now reads ["numpy", "cupy", "dask", "dask+cupy"], matching pathfinding.py::AStarSearch. The dask+cupy type came into common.get_xr_dataarray with #3706.

The bounded dask branch is now reached. _surface_distance_dask picks between map_overlap and the iterative tile Dijkstra based on max_distance; every existing benchmark call used the default np.inf and landed on the iterative fallback, emitting using iterative tile Dijkstra each time. time_surface_distance_bounded picks a max_distance whose pixel radius stays inside one chunk, following the arithmetic and the comment already in cost_distance.py. Verified by asserting the warning fires for the unbounded method and not for the bounded one:

time_surface_distance_bounded    iterative-branch warning=False (expected False)
time_surface_distance            iterative-branch warning=True (expected True)

Dask results are computed. The old file never called .compute(), so the dask numbers were real only because _sd_dask_iterative happens to evaluate eagerly while building its graph. A bounded benchmark added the same way would have measured 5 ms of graph construction. flood.py, twi.py, interpolate.py and flow_length.py all use the if hasattr(result.data, 'compute') guard; this file now does too.

Source pixels are sparse. get_xr_dataarray(is_int=True) draws rng.integers(-nx, nx), and _seed_sources treats any non-zero finite pixel as a target, so at nx=1000 499769 of 500000 pixels seeded the search at distance zero. new_cost < dist[vr, vc] was never true and the relaxation body never executed. Maximum distance found was 3.62 map units on a 360x180 degree grid. With scattered point sources it is 160.44, and the call is 6.4x slower.

Geodesic has its own class. _dijkstra_geodesic and _precompute_dd_grid were untouched; they cost about 2.2x the planar path. numpy only, since the module raises NotImplementedError for geodesic on cupy, dask, and dask+cupy.

Verification

All 50 parameter combinations executed on this host (which has CUDA and cupy), 13.01s for a single pass:

numpy      nx= 1000 time_surface_distance              0.1064s
numpy      nx= 1000 time_surface_distance_bounded      0.0101s
cupy       nx= 1000 time_surface_distance              0.2155s
cupy       nx= 1000 time_surface_distance_bounded      0.0130s
dask       nx= 1000 time_surface_distance              1.2930s
dask       nx= 1000 time_surface_distance_bounded      0.0261s
dask+cupy  nx= 1000 time_surface_distance              1.3228s
dask+cupy  nx= 1000 time_surface_distance_bounded      0.1611s
numpy      nx= 1000 time_surface_distance_geodesic     0.2828s

xrspatial/tests/test_surface_distance.py is 37 passed. flake8 clean.

Verification was by direct invocation of the benchmark classes rather than asv run, because asv discover currently fails suite-wide in this conda environment on an asv_runner metadata bug under Python 3.14 (dist.metadata["Name"] is None for some installed distribution). It fails identically on unmodified files such as CostDistance, so it is an environment problem, not something this branch introduces.

One thing to expect

Changing the source pattern will show as a one-time step change in the asv history for time_surface_distance, time_surface_allocation and time_surface_direction. The old numbers were timing a heap drain.

Not fixed here

nx=100 gives a 50x100 grid where setup dominates. The ladder up to nx=1000 still shows scaling, so it is recorded in the sweep state rather than changed.

Found by /sweep-benchmarks.