surface_distance/allocation/direction docstrings have no examples, no backend statement, and no geodesic caveat

Start in xrspatial/surface_distance.py with the docstrings for surface_distance(), surface_allocation(), and surface_direction(), comparing the nearby cost_distance() examples and backend wording. Add runnable NumPy examples, backend support text, and the geodesic limitation, then add the requested guards in xrspatial/tests/test_surface_distance.py and run those tests.

Summary

The three public functions in xrspatial/surface_distance.py —
surface_distance(), surface_allocation(), surface_direction() — are
the only proximity-family functions with no runnable example, no statement
of which array backends they support, and no note that method='geodesic'
works on one backend out of four. grep -c '>>>' xrspatial/surface_distance.py
returns 0.

Found by /sweep-documentation. Three separate gaps, all in the same three
docstrings.

1. No Examples section (Cat 1)

None of the three public functions has an Examples section. Their nearest
sibling, cost_distance(), has one (xrspatial/cost_distance.py:1264-1294),
and so do proximity(), allocation(), and direction(). A user reading
the surface-distance API reference page gets a parameter table and nothing
to copy.

2. No backend-support statement (Cat 5)

cost_distance() opens with:

Cost-distance supports NumPy, CuPy, Dask with NumPy, and Dask with CuPy
backed xarray DataArray. The return values of cost_distance are of the
same type as the input type [...]

surface_distance() and friends say nothing, even though _compute()
(xrspatial/surface_distance.py:1362-1408) dispatches to all four. Verified
on a CUDA host — all twelve combinations of the three functions and the four
backends run:

numpy        surface_distance     declared=float32 computed=float32
numpy        surface_allocation   declared=float32 computed=float32
numpy        surface_direction    declared=float32 computed=float32
cupy         surface_distance     declared=float32 computed=float32
...
dask+cupy    surface_direction    declared=float32 computed=float32

xrspatial/tests/test_proximity.py::test_docstring_states_all_backends and
test_cost_distance.py::test_docstring_states_all_backends pin this wording
for the neighbouring modules; test_surface_distance.py has no equivalent.

3. method='geodesic' is documented without its backend limitation (Cat 5)

The method parameter is described as:

method : str, default='planar'
    ``'planar'`` uses cell sizes in map units.
    ``'geodesic'`` computes great-circle horizontal distances
    from lat/lon coordinates (elevation in meters).

No caveat. In practice geodesic runs on numpy only; the other three
backends raise:

numpy        surface_distance     OK dtype=float32
cupy         surface_distance     RAISED NotImplementedError: geodesic mode is not yet supported for CuPy arrays
dask+numpy   surface_distance     RAISED NotImplementedError: geodesic mode is not yet supported for Dask arrays
dask+cupy    surface_distance     RAISED NotImplementedError: geodesic mode is not yet supported for Dask+CuPy arrays

(xrspatial/surface_distance.py:1370-1372, :1378-1380, :1400-1402.)

Large terrain rasters are the case people reach for dask on, so this is the
combination most likely to be tried and the one the docstring is silent
about.

Proposed fix

Docstrings and reference page only, no behaviour change:

• Add a backend paragraph to all three functions, matching the
  cost_distance() wording.
• Note on method that geodesic requires a NumPy-backed DataArray and name
  what the others raise.
• Add a .. sourcecode:: python Examples block to each function with output
  pinned to a real run.
• Add docstring guard tests to test_surface_distance.py in the shape of
  the existing test_docstring_states_all_backends /
  test_docstring_has_examples_section tests.

Not in scope: surface_direction()'s wrong dask result, filed separately as
#3713. The example added here uses the numpy backend, which is correct.