Added Foundational Hydrology Tools (#921)

* Add D8 flow direction module: foundational hydrology operation

Implements D8 flow direction using ESRI power-of-2 encoding (compatible
with GDAL/ArcGIS) across all four backends (numpy, cupy, dask+numpy,
dask+cupy). Follows the slope.py pattern with explicit cellsize-
parameterized backend wrappers.

New files:
- xrspatial/flow_direction.py — CPU kernel (@ngjit), GPU device/global
  kernels, four backend wrappers, public API with @supports_dataset
- xrspatial/tests/test_flow_direction.py — 66 tests (flat, cardinal,
  diagonal, known bowl, NaN handling, cross-backend, boundary modes,
  dtype acceptance, dataset support, cellsize effect, valid codes)
- benchmarks/benchmarks/flow_direction.py — ASV benchmark

Modified files:
- xrspatial/__init__.py — export flow_direction
- xrspatial/accessor.py — add Hydrology section to both accessors
- README.md — add Hydrology section with Flow Direction entry

* Add stream order module: Strahler and Shreve ordering of drainage networks

* Add BoundaryStore: memmap-backed boundary strip storage for dask tile sweeps

* Add D8 flow accumulation module with iterative dask tile sweep

* Add depression fill module using Planchon-Darboux iterative flooding

* Add sink identification module with connected-component labeling

* Add watershed delineation and drainage basin labeling module

* Add hydrology table entries to README, fix minor docstring wording

* Add standalone basin module, fix BoundaryStore temp file leaks

Extract basin delineation from watershed.py into its own basin.py
module with a basin() public function. The old basins() stays as a
backward-compatible wrapper. Add basin to __init__.py and both
xarray accessors.

Add BoundarySnapshot to _boundary_store.py -- a read-only in-memory
copy of converged boundary strips. BoundaryStore.snapshot() copies
strip data to plain numpy arrays and closes the underlying memmap
temp files. Apply this pattern across all hydrology dask backends
(flow_accumulation, watershed, fill, stream_order) so temp
directories are cleaned up before the lazy dask result is returned.
Previously, BoundaryStore objects captured in map_blocks closures
kept temp files alive indefinitely.

* Add stream link segmentation module with all four backends

Assigns unique position-based IDs to each stream segment between
junctions, headwaters, and outlets. Supports numpy, cupy, dask+numpy,
and dask+cupy backends.

* Add snap pour point module with all four backends

Moves each pour point to the highest flow-accumulation cell within a
circular search radius so that watershed delineation starts from the
actual drainage channel.  Dask backend extracts sparse pour points
chunk-by-chunk (map_blocks flag pass + selective load) to keep memory
bounded regardless of grid size.

* Add flow path tracing module with all four backends

* Add D-infinity flow direction module with all four backends

* Fix flow_path Dask OOM: vectorized chunk grouping, numpy buffers, adaptive cache

Replace three memory/performance bottlenecks in _flow_path_dask:

- Path tracing (phase 3): growable numpy buffers (~24 bytes/cell)
  instead of list-of-tuples (~164 bytes/cell)
- Output assembly (phase 4): pre-group cells by chunk via vectorized
  searchsorted + stable argsort, then O(1) dict lookup per block
  instead of O(n_chunks * n_cells) linear scan
- LRU cache: adaptive sizing capped at ~512 MB instead of fixed 32
  entries regardless of chunk size

* Add flow length, TWI, and HAND modules with tests and notebook

* Add benchmarks for flow_length, TWI, and HAND

Author: brendancol

README.md

@@ -239,6 +239,22 @@ In the GIS world, rasters are used for representing continuous phenomena (e.g. e
 
 -----------
 
+### **Hydrology**
+
+| Name | Description | NumPy xr.DataArray | Dask xr.DataArray | CuPy GPU xr.DataArray | Dask GPU xr.DataArray |
+|:----------:|:------------|:----------------------:|:--------------------:|:-------------------:|:------:|
+| [Flow Direction (D8)](xrspatial/flow_direction.py) | Computes D8 flow direction from each cell toward the steepest downhill neighbor | ✅️ | ✅️ | ✅️ | ✅️ |
+| [Flow Direction (Dinf)](xrspatial/flow_direction_dinf.py) | Computes D-infinity flow direction as a continuous angle toward the steepest downslope facet | ✅️ | ✅️ | ✅️ | ✅️ |
+| [Flow Accumulation (D8)](xrspatial/flow_accumulation.py) | Counts upstream cells draining through each cell in a D8 flow direction grid | ✅️ | ✅️ | ✅️ | ✅️ |
+| [Watershed](xrspatial/watershed.py) | Labels each cell with the pour point it drains to via D8 flow direction | ✅️ | ✅️ | ✅️ | ✅️ |
+| [Basins](xrspatial/watershed.py) | Delineates drainage basins by labeling each cell with its outlet ID | ✅️ | ✅️ | ✅️ | ✅️ |
+| [Stream Order](xrspatial/stream_order.py) | Assigns Strahler or Shreve stream order to cells in a drainage network | ✅️ | ✅️ | ✅️ | ✅️ |
+| [Stream Link](xrspatial/stream_link.py) | Assigns unique IDs to each stream segment between junctions | ✅️ | ✅️ | ✅️ | ✅️ |
+| [Snap Pour Point](xrspatial/snap_pour_point.py) | Snaps pour points to the highest-accumulation cell within a search radius | ✅️ | ✅️ | ✅️ | ✅️ |
+| [Flow Path](xrspatial/flow_path.py) | Traces downstream flow paths from start points through a D8 direction grid | ✅️ | ✅️ | ✅️ | ✅️ |
+
+-----------
+
 ### **Zonal**
 
 | Name | Description | NumPy xr.DataArray | Dask xr.DataArray | CuPy GPU xr.DataArray | Dask GPU xr.DataArray |

benchmarks/benchmarks/flow_accumulation.py

@@ -0,0 +1,28 @@
+from xrspatial import flow_direction, flow_accumulation
+
+from .common import Benchmarking, get_xr_dataarray
+
+
+class FlowAccumulation:
+    params = ([100, 300, 1000], ["numpy", "dask"])
+    param_names = ("nx", "type")
+
+    def setup(self, nx, type):
+        ny = nx // 2
+        elev = get_xr_dataarray((ny, nx), type)
+        # Compute flow_dir from elevation (materialise for benchmarking)
+        self.flow_dir = flow_direction(elev)
+        if hasattr(self.flow_dir.data, 'compute'):
+            self.flow_dir.data = self.flow_dir.data.compute()
+            # Re-chunk for dask benchmark
+            if type == "dask":
+                import dask.array as da
+                self.flow_dir.data = da.from_array(
+                    self.flow_dir.data,
+                    chunks=(max(1, ny // 2), max(1, nx // 2)),
+                )
+
+    def time_flow_accumulation(self, nx, type):
+        result = flow_accumulation(self.flow_dir)
+        if hasattr(result.data, 'compute'):
+            result.data.compute()

benchmarks/benchmarks/flow_direction.py

@@ -0,0 +1,11 @@
+from xrspatial import flow_direction
+
+from .common import Benchmarking
+
+
+class FlowDirection(Benchmarking):
+    def __init__(self):
+        super().__init__(func=flow_direction)
+
+    def time_flow_direction(self, nx, type):
+        return self.time(nx, type)

benchmarks/benchmarks/flow_length.py

@@ -0,0 +1,31 @@
+from xrspatial import flow_direction, flow_length
+
+from .common import get_xr_dataarray
+
+
+class FlowLength:
+    params = ([100, 300, 1000], ["numpy", "dask"])
+    param_names = ("nx", "type")
+
+    def setup(self, nx, type):
+        ny = nx // 2
+        elev = get_xr_dataarray((ny, nx), "numpy")
+        self.flow_dir = flow_direction(elev)
+        fd_data = self.flow_dir.data
+
+        if type == "dask":
+            import dask.array as da
+            self.flow_dir.data = da.from_array(
+                fd_data,
+                chunks=(max(1, ny // 2), max(1, nx // 2)),
+            )
+
+    def time_flow_length_downstream(self, nx, type):
+        result = flow_length(self.flow_dir, direction='downstream')
+        if hasattr(result.data, 'compute'):
+            result.data.compute()
+
+    def time_flow_length_upstream(self, nx, type):
+        result = flow_length(self.flow_dir, direction='upstream')
+        if hasattr(result.data, 'compute'):
+            result.data.compute()

benchmarks/benchmarks/hand.py

@@ -0,0 +1,31 @@
+from xrspatial import flow_direction, flow_accumulation, hand
+
+from .common import get_xr_dataarray
+
+
+class HAND:
+    params = ([100, 300, 1000], ["numpy", "dask"])
+    param_names = ("nx", "type")
+
+    def setup(self, nx, type):
+        ny = nx // 2
+        self.elev = get_xr_dataarray((ny, nx), "numpy")
+        flow_dir = flow_direction(self.elev)
+        self.flow_dir = flow_dir
+        self.flow_accum = flow_accumulation(flow_dir)
+
+        if type == "dask":
+            import dask.array as da
+            chunks = (max(1, ny // 2), max(1, nx // 2))
+            fd_data = self.flow_dir.data
+            fa_data = self.flow_accum.data
+            elev_data = self.elev.data
+
+            self.flow_dir.data = da.from_array(fd_data, chunks=chunks)
+            self.flow_accum.data = da.from_array(fa_data, chunks=chunks)
+            self.elev.data = da.from_array(elev_data, chunks=chunks)
+
+    def time_hand(self, nx, type):
+        result = hand(self.flow_dir, self.flow_accum, self.elev)
+        if hasattr(result.data, 'compute'):
+            result.data.compute()

benchmarks/benchmarks/twi.py

@@ -0,0 +1,30 @@
+from xrspatial import flow_direction, flow_accumulation, slope, twi
+
+from .common import get_xr_dataarray
+
+
+class TWI:
+    params = ([100, 300, 1000], ["numpy", "dask"])
+    param_names = ("nx", "type")
+
+    def setup(self, nx, type):
+        ny = nx // 2
+        elev = get_xr_dataarray((ny, nx), "numpy")
+        flow_dir = flow_direction(elev)
+        self.flow_accum = flow_accumulation(flow_dir)
+        self.slope_agg = slope(elev)
+
+        if type == "dask":
+            import dask.array as da
+            chunks = (max(1, ny // 2), max(1, nx // 2))
+            self.flow_accum.data = da.from_array(
+                self.flow_accum.data, chunks=chunks,
+            )
+            self.slope_agg.data = da.from_array(
+                self.slope_agg.data, chunks=chunks,
+            )
+
+    def time_twi(self, nx, type):
+        result = twi(self.flow_accum, self.slope_agg)
+        if hasattr(result.data, 'compute'):
+            result.data.compute()

benchmarks/benchmarks/watershed.py

@@ -0,0 +1,64 @@
+from xrspatial import flow_direction, watershed, basins
+
+from .common import get_xr_dataarray
+
+import numpy as np
+
+
+class Watershed:
+    params = ([100, 300, 1000], ["numpy", "dask"])
+    param_names = ("nx", "type")
+
+    def setup(self, nx, type):
+        ny = nx // 2
+        elev = get_xr_dataarray((ny, nx), "numpy")
+        self.flow_dir = flow_direction(elev)
+        fd_data = self.flow_dir.data
+
+        # Create pour_points at pits (code 0)
+        pp = np.full_like(fd_data, np.nan)
+        pit_mask = fd_data == 0
+        pp[pit_mask] = np.arange(1, pit_mask.sum() + 1, dtype=np.float64)
+
+        if type == "dask":
+            import dask.array as da
+            self.flow_dir.data = da.from_array(
+                fd_data,
+                chunks=(max(1, ny // 2), max(1, nx // 2)),
+            )
+            self.pour_points = self.flow_dir.copy()
+            self.pour_points.data = da.from_array(
+                pp,
+                chunks=(max(1, ny // 2), max(1, nx // 2)),
+            )
+        else:
+            self.pour_points = self.flow_dir.copy()
+            self.pour_points.data = pp
+
+    def time_watershed(self, nx, type):
+        result = watershed(self.flow_dir, self.pour_points)
+        if hasattr(result.data, 'compute'):
+            result.data.compute()
+
+
+class Basins:
+    params = ([100, 300, 1000], ["numpy", "dask"])
+    param_names = ("nx", "type")
+
+    def setup(self, nx, type):
+        ny = nx // 2
+        elev = get_xr_dataarray((ny, nx), "numpy")
+        self.flow_dir = flow_direction(elev)
+        fd_data = self.flow_dir.data
+
+        if type == "dask":
+            import dask.array as da
+            self.flow_dir.data = da.from_array(
+                fd_data,
+                chunks=(max(1, ny // 2), max(1, nx // 2)),
+            )
+
+    def time_basins(self, nx, type):
+        result = basins(self.flow_dir)
+        if hasattr(result.data, 'compute'):
+            result.data.compute()

examples/user_guide/11_Hydrology.ipynb

@@ -14,7 +14,14 @@
 from xrspatial.classify import reclassify  # noqa
 from xrspatial.curvature import curvature  # noqa
 from xrspatial.emerging_hotspots import emerging_hotspots  # noqa
+from xrspatial.fill import fill  # noqa
+from xrspatial.flow_accumulation import flow_accumulation  # noqa
+from xrspatial.flow_direction import flow_direction  # noqa
+from xrspatial.flow_direction_dinf import flow_direction_dinf  # noqa
+from xrspatial.flow_length import flow_length  # noqa
+from xrspatial.flow_path import flow_path  # noqa
 from xrspatial.focal import mean  # noqa
+from xrspatial.hand import hand  # noqa
 from xrspatial.hillshade import hillshade  # noqa
 from xrspatial.mahalanobis import mahalanobis  # noqa
 from xrspatial.multispectral import arvi  # noqa
@@ -31,6 +38,10 @@
 from xrspatial.proximity import great_circle_distance  # noqa
 from xrspatial.proximity import manhattan_distance  # noqa
 from xrspatial.proximity import proximity  # noqa
+from xrspatial.sink import sink  # noqa
+from xrspatial.snap_pour_point import snap_pour_point  # noqa
+from xrspatial.stream_link import stream_link  # noqa
+from xrspatial.stream_order import stream_order  # noqa
 from xrspatial.slope import slope  # noqa
 from xrspatial.surface_distance import surface_allocation  # noqa
 from xrspatial.surface_distance import surface_direction  # noqa
@@ -39,8 +50,12 @@
 from xrspatial.terrain_metrics import roughness  # noqa
 from xrspatial.terrain_metrics import tpi  # noqa
 from xrspatial.terrain_metrics import tri  # noqa
+from xrspatial.twi import twi  # noqa
 from xrspatial.polygonize import polygonize  # noqa
 from xrspatial.viewshed import viewshed  # noqa
+from xrspatial.basin import basin  # noqa
+from xrspatial.watershed import basins  # noqa
+from xrspatial.watershed import watershed  # noqa
 from xrspatial.zonal import apply as zonal_apply  # noqa
 from xrspatial.zonal import crop  # noqa
 from xrspatial.zonal import trim  # noqa