surface_distance: max_distance=nan and negative max_distance give backend-dependent results with no error

[brendancol]

surface_distance, surface_allocation and surface_direction never validate
max_distance or the shape of target_values. Both gaps show up as
backend-dependent behavior on inputs that a caller can easily produce by
accident.

1. max_distance=nan gives three different answers on three backends

The numpy kernel bounds the search with a rejection test
(if cost_u > max_distance: break, surface_distance.py:237) which is False for
NaN, so it runs to completion and returns a full unbounded distance surface.
The CUDA kernel uses the acceptance form (if best < current and best <= max_distance, surface_distance.py:533) which is also False for NaN, so nothing
relaxes and only the seed pixels survive. dask follows the numpy kernel but
also emits the "max_distance is infinite" warning, which is not what the caller
passed.

Nothing raises and nothing warns. The output looks like a normal distance
raster in every case.

2. max_distance negative is silently honoured on numpy/cupy and a dask internal error otherwise

numpy and cupy return an all-NaN raster (the source pixels at distance 0 are
masked too, since 0 > -5 is False but _finalize_dist compares
dist > max_distance). dask raises ValueError: length should not be negative
out of da.map_overlap, whose message never mentions max_distance.

Repro

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

src = np.zeros((6, 6)); src[0, 0] = 1.0
elev = np.zeros((6, 6))
coords = {'y': np.arange(6.0), 'x': np.arange(6.0)}

r_np = xr.DataArray(src, dims=('y', 'x'), coords=coords)
e_np = xr.DataArray(elev, dims=('y', 'x'), coords=coords)
r_cp = xr.DataArray(cp.asarray(src), dims=('y', 'x'), coords=coords)
e_cp = xr.DataArray(cp.asarray(elev), dims=('y', 'x'), coords=coords)
r_dk = xr.DataArray(dsa.from_array(src, chunks=(3, 3)), dims=('y', 'x'), coords=coords)
e_dk = xr.DataArray(dsa.from_array(elev, chunks=(3, 3)), dims=('y', 'x'), coords=coords)

surface_distance(r_np, e_np, max_distance=np.nan)          # 36/36 finite, max 7.07
surface_distance(r_cp, e_cp, max_distance=np.nan)          #  1/36 finite, max 0.0
surface_distance(r_dk, e_dk, max_distance=np.nan).compute()# 36/36 finite, max 7.07

surface_distance(r_np, e_np, max_distance=-5.0)            # all NaN
surface_distance(r_cp, e_cp, max_distance=-5.0)            # all NaN
surface_distance(r_dk, e_dk, max_distance=-5.0).compute()  # ValueError: length should not be negative

Observed on this host (CUDA available, numba 0.x, dask 2026.x):

--- max_distance=np.nan ---
numpy                        finite=36/36 max=7.071067810058594
cupy                         finite= 1/36 max=0.0
dask+numpy                   finite=36/36 max=7.071067810058594
--- max_distance=-5.0 ---
numpy                        finite= 0/36 max=n/a
cupy                         finite= 0/36 max=n/a
dask+numpy                   ValueError: length should not be negative

3. Non-1D target_values crashes inside numba

target_values is passed through np.asarray(..., dtype=np.float64) with no
shape check (surface_distance.py:1318). A scalar (target_values=1 instead of
[1]) reaches _seed_sources as a 0-d array and fails there:

TypingError: Failed in nopython mode pipeline (step: nopython frontend)
No implementation of function Function(<built-in function getitem>) found for signature:
 >>> getitem(array(float64, 0d, C), int64)
There are 22 candidate implementations:
 ...

followed by roughly a hundred lines of overload candidates. The cupy path fails
differently, with TypeError: len() of unsized object. Neither message names
the parameter.

Precedent

proximity() takes the same two parameters and already validates both, in
xrspatial/proximity.py:1396:

if np.isnan(max_distance) or max_distance < 0:
    raise ValueError(
        "max_distance must be non-negative, got {0!r}.".format(max_distance)
    )

with the finite-target_values check just above it (added for #2850, the same
"backend-dependent output from an unvalidated parameter" problem). The surface
functions were written against the same API and did not pick up the guards.

Proposed fix

Validate in _compute(), next to the existing connectivity and method
checks, and reuse proximity's wording:

• ValueError when max_distance is NaN or negative
• ValueError when target_values is not 1-D, naming the parameter

Both are new exceptions on inputs that previously produced silent
backend-dependent output or an internal error, so no documented usage changes.
The existing 37 tests stay green.

Note for a separate change: proximity(raster, target_values=1) also fails with
TypeError: len() of unsized object, so the 0-d target_values rough edge is
shared. Only the surface functions are touched here.

Found by /sweep-error-handling.