surface_allocation/surface_direction tie-break is backend-dependent and undocumented

Locate the docstrings for surface_allocation() and surface_direction(), then compare their wording with proximity.allocation's documented tie-break rule. Confirm the intended choice with maintainers: document the backend-dependent, unspecified tie-break or adopt the lowest-flat-index rule; done means the selected behavior is clearly stated in the relevant docstrings.

surface_allocation() and surface_direction() pick a different source
on each backend when two targets are exactly equidistant, and neither
docstring says the tie-break is unspecified.

Found while validating surface_distance against
scipy.sparse.csgraph.dijkstra. The distances agree bit for bit on every
backend; only the reported nearest source differs. proximity.allocation
documents its rule ("among equidistant targets, the lowest flat index
wins"); the surface trio documents nothing, so a caller has no way to know
the answer is backend-dependent.

Reproduction

A Gaussian hill with two sources placed symmetrically about the main
diagonal, so a band of pixels is exactly equidistant from both:

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

n = 64
y, x = np.mgrid[0:n, 0:n].astype(float)
c = (n - 1) / 2
elev = 400.0 * np.exp(-(((x - c) ** 2 + (y - c) ** 2) / (2 * 12.0 ** 2)))
src = np.zeros((n, n)); src[2, 2] = 1.0; src[n - 3, n - 3] = 2.0

def da(a):
    return xr.DataArray(a, dims=('y', 'x'),
                        coords={'y': np.arange(n)[::-1] * 30.0,
                                'x': np.arange(n) * 30.0})

a_np = surface_allocation(da(src), da(elev)).data
a_cu = surface_allocation(da(cp.asarray(src)), da(cp.asarray(elev))).data.get()
a_dk = surface_allocation(da(dsa.from_array(src, chunks=(32, 32))),
                          da(dsa.from_array(elev, chunks=(32, 32)))).data.compute()

for name, a in [("cupy", a_cu), ("dask", a_dk)]:
    d = (a != a_np) & ~(np.isnan(a) & np.isnan(a_np))
    print(f"numpy vs {name}: {int(d.sum())} differing allocation px of {n * n}")

numpy vs cupy: 29 differing allocation px of 4096
numpy vs dask: 25 differing allocation px of 4096

Every differing pixel is an exact tie. Computing the distance to each
source separately gives max |d1 - d2| == 0.0 over the differing pixels,
so this is only the tie-break and not a distance error.

Why it happens

The three backends resolve ties by three unrelated mechanisms:

• numpy / dask: _dijkstra copies alloc, src_row, src_col from
  whichever neighbour relaxed the pixel first, which follows the binary
  heap's pop order among equal keys.
• cupy: _sd_relax_kernel uses a strict new_cost < best, so the winner
  is whichever neighbour happened to be relaxed first in the pass that
  converged.
• the dask iterative path additionally depends on tile sweep order.

None of these is a stated contract.

Suggested fix

Documentation is enough. Either state that the tie-break is unspecified
and backend-dependent, or adopt proximity.allocation's lowest-flat-index
rule across all four backends and document that. The first is a one-line
docstring change; the second changes results and needs a decision on
whether matching proximity is worth the GPU cost.

Severity: MEDIUM. Distances are correct everywhere; only the reported
source label is unstable, and only on exact ties.