surface_distance tests never touch dask+cupy and only call surface_direction on numpy

Start with xrspatial/tests/test_surface_distance.py and the referenced branches in surface_distance.py, especially _surface_distance_dask_cupy(), surface_direction(), and the dask iterative helpers. Run the stated NUMBA_DISABLE_JIT=1 pytest coverage command, then add backend, edge-case, parameter, error-path, Dataset, and metadata tests until the listed gaps are exercised.

Summary

xrspatial/tests/test_surface_distance.py has 37 tests covering three
public functions across four backends. The matrix has large holes: the
dask+cupy backend is never exercised at all, and surface_direction() is
only ever called on eager numpy. Measured branch coverage of
xrspatial.surface_distance is 80% (570/695 statements).

The surface_direction() hole is not theoretical. It hid #3713, in which
the function returns wrong compass bearings on every dask raster that
takes the iterative path.

What is missing

Backend coverage (Cat 1)

_surface_distance_dask_cupy() (surface_distance.py:1232-1289) has zero
coverage. Neither of its two branches is reached: the bounded
map_overlap branch nor the unbounded branch that converts to dask+numpy
and back.

surface_direction() is called by four tests, all on eager numpy. No test
runs it on cupy, dask+numpy, or dask+cupy. As a result the cupy DIRECTION
branch (:602-616), which round-trips through
_vectorized_calc_direction() on the host, is unexercised, and so is the
dask assembly path that #3713 breaks.

Two cupy seeding paths are also unreached: the target_values branch
(:560-561) and the early return when a chunk holds no sources
(:571-573).

NaN / Inf edge cases (Cat 2)

test_nan_barrier covers a NaN wall on numpy and dask+numpy only. Nothing
tests Inf elevation, which the np.isfinite() guards treat as a barrier
the same way, and nothing tests an entirely NaN elevation surface.

Geometric edge cases (Cat 3)

No 1x1 raster test. No Nx1 column-strip test; the existing 1xN row strips
(test_45_degree_slope, test_target_values,
test_multiple_sources_nearest_wins) only exercise the horizontal
degenerate case.

Parameter coverage (Cat 4)

connectivity=4 is tested once, on eager numpy. The dask iterative path
branches on connectivity in _compute_seeds_sd() and _can_skip_sd()
(:919-932), and those branches never run with connectivity=4.

The dask iterative path is never run with target_values set, leaving
_preprocess_tiles_sd() :698-701 and _run_tile() :956-959 uncovered.

Three error paths have no test: the raster.dims should be rejection
(:1306) and the three per-backend NotImplementedErrors for geodesic
mode on cupy, dask, and dask+cupy (:1371, :1379, :1401).

All three public functions carry @supports_dataset, and none is ever
called with a Dataset. This is the same gap the proximity trio has.

Metadata preservation (Cat 5)

No test asserts that input attrs, coords, or dims survive the call.
_compute() reads res off the input attrs via
get_dataarray_resolution() to build its edge costs, so a chained call
depends on them being carried through.

Reproduction

Coverage before, on this host:

NUMBA_DISABLE_JIT=1 python -m pytest xrspatial/tests/test_surface_distance.py \
    --cov=xrspatial.surface_distance --cov-branch -q

80% (570/695 statements, 208/276 branches)

NUMBA_DISABLE_JIT=1 is needed because pytest-cov and numba's extension
registration collide, and it also makes the GPU tests fail, so the figure
understates real coverage. It is still the right before/after yardstick.

Proof that dask+cupy is reachable and currently untested:

import numpy as np, xarray as xr, dask.array as da, cupy
from xrspatial import surface_distance

source = np.zeros((8, 10)); source[2, 3] = 1.0; source[6, 7] = 2.0
elev = np.random.default_rng(42).uniform(0, 100, (8, 10))


def mk(a, chunks=None, gpu=False):
    d = xr.DataArray(a, dims=['y', 'x'],
                     coords={'y': np.arange(a.shape[0], dtype=float),
                             'x': np.arange(a.shape[1], dtype=float)},
                     attrs={'res': (1.0, 1.0)})
    if gpu:
        d.data = cupy.asarray(d.data)
    if chunks:
        d.data = da.from_array(d.data, chunks=chunks)
    return d


ref = surface_distance(mk(source), mk(elev), max_distance=15.0).values
got = surface_distance(mk(source, (4, 5), True), mk(elev, (4, 5), True),
                       max_distance=15.0).compute().data.get()
print(np.nanmax(np.abs(got - ref)))

0.0

The backend works. Nothing was checking that it does.

Proposed fix

Add tests for each hole above, parametrised over all four backends where
the path exists. CUDA is available on this host, so the cupy and dask+cupy
tests run rather than skip.

Provenance

Found by /sweep-test-coverage on surface_distance.