Add surface_distance module: 3D terrain distance via multi-source Dijkstra (#918)

* Add surface_distance module: 3D terrain distance via multi-source Dijkstra

Adds surface_distance(), surface_allocation(), and surface_direction()
which compute distance, nearest-target allocation, and compass direction
along a 3D terrain surface (DEM). Unlike proximity (flat 2D Euclidean)
or cost_distance (scalar friction), edge costs account for elevation
relief: sqrt(horizontal_dist² + dz²).

Backends: NumPy (Dijkstra), CuPy (parallel Bellman-Ford relaxation),
Dask+NumPy bounded (map_overlap) and unbounded (iterative boundary
Dijkstra), Dask+CuPy. Geodesic mode (lat/lon → great-circle horizontal
distances) supported for NumPy.

Includes 32 tests, ASV benchmarks, __init__ exports, and .xrs accessor
methods on both DataArray and Dataset.

* Add surface distance functions to README proximity table

Author: brendancol

README.md

@@ -206,6 +206,9 @@ In the GIS world, rasters are used for representing continuous phenomena (e.g. e
 | [Cost Distance](xrspatial/cost_distance.py) | Computes minimum accumulated cost to the nearest source through a friction surface | ✅️ | ✅ | ✅️ | ✅️ |
 | [Direction](xrspatial/proximity.py) | Computes the direction from each cell to the nearest source feature | ✅️ | ✅ | ✅️ | ✅️ |
 | [Proximity](xrspatial/proximity.py) | Computes the distance from each cell to the nearest source feature | ✅️ | ✅ | ✅️ | ✅️ |
+| [Surface Distance](xrspatial/surface_distance.py) | Computes distance along the 3D terrain surface to the nearest source | ✅️ | ✅ | ✅️ | ✅️ |
+| [Surface Allocation](xrspatial/surface_distance.py) | Assigns each cell to the nearest source by terrain surface distance | ✅️ | ✅ | ✅️ | ✅️ |
+| [Surface Direction](xrspatial/surface_distance.py) | Computes compass direction to the nearest source by terrain surface distance | ✅️ | ✅ | ✅️ | ✅️ |
 
 --------
 

benchmarks/benchmarks/surface_distance.py

@@ -0,0 +1,26 @@
+import numpy as np
+
+from xrspatial.surface_distance import (
+    surface_distance, surface_allocation, surface_direction,
+)
+
+from .common import get_xr_dataarray
+
+
+class SurfaceDistance:
+    params = ([100, 300, 1000], ["numpy", "dask"])
+    param_names = ("nx", "type")
+
+    def setup(self, nx, type):
+        ny = nx // 2
+        self.agg = get_xr_dataarray((ny, nx), type, is_int=True)
+        self.elev = get_xr_dataarray((ny, nx), type)
+
+    def time_surface_distance(self, nx, type):
+        surface_distance(self.agg, self.elev)
+
+    def time_surface_allocation(self, nx, type):
+        surface_allocation(self.agg, self.elev)
+
+    def time_surface_direction(self, nx, type):
+        surface_direction(self.agg, self.elev)

xrspatial/__init__.py

@@ -32,6 +32,9 @@
 from xrspatial.proximity import manhattan_distance  # noqa
 from xrspatial.proximity import proximity  # noqa
 from xrspatial.slope import slope  # noqa
+from xrspatial.surface_distance import surface_allocation  # noqa
+from xrspatial.surface_distance import surface_direction  # noqa
+from xrspatial.surface_distance import surface_distance  # noqa
 from xrspatial.terrain import generate_terrain  # noqa
 from xrspatial.polygonize import polygonize  # noqa
 from xrspatial.viewshed import viewshed  # noqa

xrspatial/accessor.py

@@ -97,6 +97,18 @@ def cost_distance(self, friction, **kwargs):
         from .cost_distance import cost_distance
         return cost_distance(self._obj, friction, **kwargs)
 
+    def surface_distance(self, elevation, **kwargs):
+        from .surface_distance import surface_distance
+        return surface_distance(self._obj, elevation, **kwargs)
+
+    def surface_allocation(self, elevation, **kwargs):
+        from .surface_distance import surface_allocation
+        return surface_allocation(self._obj, elevation, **kwargs)
+
+    def surface_direction(self, elevation, **kwargs):
+        from .surface_distance import surface_direction
+        return surface_direction(self._obj, elevation, **kwargs)
+
     # ---- Pathfinding ----
 
     def a_star_search(self, start, goal, **kwargs):
@@ -259,6 +271,18 @@ def cost_distance(self, friction, **kwargs):
         from .cost_distance import cost_distance
         return cost_distance(self._obj, friction, **kwargs)
 
+    def surface_distance(self, elevation, **kwargs):
+        from .surface_distance import surface_distance
+        return surface_distance(self._obj, elevation, **kwargs)
+
+    def surface_allocation(self, elevation, **kwargs):
+        from .surface_distance import surface_allocation
+        return surface_allocation(self._obj, elevation, **kwargs)
+
+    def surface_direction(self, elevation, **kwargs):
+        from .surface_distance import surface_direction
+        return surface_direction(self._obj, elevation, **kwargs)
+
     # ---- Multispectral (band mapping via kwargs) ----
 
     def ndvi(self, nir, red, **kwargs):