From 835b2d983c001dc177f4542a293f0708644f1b46 Mon Sep 17 00:00:00 2001
From: Brendan Collins <brendancol@gmail.com>
Date: Thu, 25 Jun 2026 12:58:11 -0700
Subject: [PATCH 1/4] Make routing wrappers the sole public hydrology API
 (#3528)

Drop the suffixed *_d8/_dinf/_mfd hydrology functions from the top-level
xrspatial namespace; they stay importable from xrspatial.hydro where the
routing wrappers dispatch to them. Give each _RoutingDispatch wrapper a
__name__ and a docstring describing the routing options. Repoint the one
top-level suffixed import in test_dask_task_names.py to xrspatial.hydro.

Adds test_routing_public_api.py covering the public surface, the routing
dispatch (previously untested), default routing, unknown-routing errors,
and stream_order ordering pass-through.
---
 xrspatial/__init__.py                         |  14 ---
 xrspatial/hydro/__init__.py                   |  62 +++++++++--
 .../hydro/tests/test_routing_public_api.py    | 100 ++++++++++++++++++
 xrspatial/tests/test_dask_task_names.py       |   3 +-
 4 files changed, 154 insertions(+), 25 deletions(-)
 create mode 100644 xrspatial/hydro/tests/test_routing_public_api.py

diff --git a/xrspatial/__init__.py b/xrspatial/__init__.py
index 198337490..cbb79ce6a 100644
--- a/xrspatial/__init__.py
+++ b/xrspatial/__init__.py
@@ -31,7 +31,6 @@
 from xrspatial.emerging_hotspots import emerging_hotspots  # noqa
 from xrspatial.erosion import erode  # noqa
 from xrspatial.hydro import fill  # noqa: unified wrapper
-from xrspatial.hydro import fill_d8  # noqa
 from xrspatial.interpolate import idw  # noqa
 from xrspatial.interpolate import kriging  # noqa
 from xrspatial.interpolate import spline  # noqa
@@ -52,13 +51,9 @@
 from xrspatial.flood import vegetation_curve_number  # noqa
 from xrspatial.flood import vegetation_roughness  # noqa
 from xrspatial.hydro import flow_accumulation  # noqa: unified wrapper
-from xrspatial.hydro import flow_accumulation_d8, flow_accumulation_dinf, flow_accumulation_mfd  # noqa
 from xrspatial.hydro import flow_direction  # noqa: unified wrapper
-from xrspatial.hydro import flow_direction_d8, flow_direction_dinf, flow_direction_mfd  # noqa
 from xrspatial.hydro import flow_length  # noqa: unified wrapper
-from xrspatial.hydro import flow_length_d8, flow_length_dinf, flow_length_mfd  # noqa
 from xrspatial.hydro import flow_path  # noqa: unified wrapper
-from xrspatial.hydro import flow_path_d8, flow_path_dinf, flow_path_mfd  # noqa
 from xrspatial.focal import mean  # noqa
 from xrspatial.glcm import glcm_texture  # noqa
 from xrspatial.morphology import morph_black_tophat  # noqa
@@ -69,7 +64,6 @@
 from xrspatial.morphology import morph_opening  # noqa
 from xrspatial.morphology import morph_white_tophat  # noqa
 from xrspatial.hydro import hand  # noqa: unified wrapper
-from xrspatial.hydro import hand_d8, hand_dinf, hand_mfd  # noqa
 from xrspatial.hillshade import hillshade  # noqa
 from xrspatial.mahalanobis import mahalanobis  # noqa
 from xrspatial.multispectral import arvi  # noqa
@@ -100,13 +94,9 @@
 from xrspatial.proximity import manhattan_distance  # noqa
 from xrspatial.proximity import proximity  # noqa
 from xrspatial.hydro import sink  # noqa: unified wrapper
-from xrspatial.hydro import sink_d8  # noqa
 from xrspatial.hydro import snap_pour_point  # noqa: unified wrapper
-from xrspatial.hydro import snap_pour_point_d8  # noqa
 from xrspatial.hydro import stream_link  # noqa: unified wrapper
-from xrspatial.hydro import stream_link_d8, stream_link_dinf, stream_link_mfd  # noqa
 from xrspatial.hydro import stream_order  # noqa: unified wrapper
-from xrspatial.hydro import stream_order_d8, stream_order_dinf, stream_order_mfd  # noqa
 from xrspatial.sieve import sieve  # noqa
 from xrspatial.sky_view_factor import sky_view_factor  # noqa
 from xrspatial.slope import slope  # noqa
@@ -122,7 +112,6 @@
 from xrspatial.terrain_metrics import tri  # noqa
 from xrspatial.validate import validate  # noqa
 from xrspatial.hydro import twi  # noqa: unified wrapper
-from xrspatial.hydro import twi_d8  # noqa
 from xrspatial.polygon_clip import clip_polygon  # noqa
 from xrspatial.polygonize import polygonize  # noqa
 from xrspatial.viewshed import viewshed  # noqa
@@ -130,11 +119,8 @@
 from xrspatial.visibility import line_of_sight  # noqa
 from xrspatial.visibility import visibility_frequency  # noqa
 from xrspatial.hydro import basin  # noqa: unified wrapper
-from xrspatial.hydro import basin_d8  # noqa
 from xrspatial.hydro import basins  # noqa: backward-compat alias
-from xrspatial.hydro import basins_d8  # noqa
 from xrspatial.hydro import watershed  # noqa: unified wrapper
-from xrspatial.hydro import watershed_d8, watershed_dinf, watershed_mfd  # noqa
 from xrspatial.zonal import apply as zonal_apply  # noqa
 from xrspatial.zonal import crop  # noqa
 from xrspatial.zonal import trim  # noqa
diff --git a/xrspatial/hydro/__init__.py b/xrspatial/hydro/__init__.py
index 1cf71cb15..2912db836 100644
--- a/xrspatial/hydro/__init__.py
+++ b/xrspatial/hydro/__init__.py
@@ -58,11 +58,21 @@ class _RoutingDispatch:
     (``'d8'``, ``'dinf'``, ``'mfd'``) rather than array type.
     """
 
-    __slots__ = ('_name', '_impls')
-
-    def __init__(self, name, **impls):
+    def __init__(self, name, summary=None, **impls):
         self._name = name
         self._impls = impls
+        self.__name__ = name
+        if summary is not None:
+            keys = ', '.join(repr(k) for k in impls)
+            variants = ', '.join(f"``{fn.__name__}``" for fn in impls.values())
+            self.__doc__ = (
+                f"{summary}\n\n"
+                f"Select the routing algorithm with the ``routing`` keyword "
+                f"(default ``'d8'``). Valid values: {keys}.\n\n"
+                f"Dispatches to the matching implementation in "
+                f"``xrspatial.hydro`` ({variants}), where the full parameter "
+                f"list and algorithm references live."
+            )
 
     def __call__(self, *args, routing='d8', **kwargs):
         try:
@@ -80,37 +90,45 @@ def __call__(self, *args, routing='d8', **kwargs):
 
 flow_direction = _RoutingDispatch(
     'flow_direction',
+    summary="Direction of steepest descent out of each cell.",
     d8=flow_direction_d8, dinf=flow_direction_dinf, mfd=flow_direction_mfd,
 )
 
 flow_accumulation = _RoutingDispatch(
     'flow_accumulation',
+    summary="Upstream cells or area draining through each cell.",
     d8=flow_accumulation_d8, dinf=flow_accumulation_dinf,
     mfd=flow_accumulation_mfd,
 )
 
 flow_length = _RoutingDispatch(
     'flow_length',
+    summary="Distance along the flow path to the outlet or from the divide.",
     d8=flow_length_d8, dinf=flow_length_dinf, mfd=flow_length_mfd,
 )
 
 flow_path = _RoutingDispatch(
     'flow_path',
+    summary="Trace downstream flow paths from a set of start points.",
     d8=flow_path_d8, dinf=flow_path_dinf, mfd=flow_path_mfd,
 )
 
 watershed = _RoutingDispatch(
     'watershed',
+    summary="Label each cell with the pour point it drains to.",
     d8=watershed_d8, dinf=watershed_dinf, mfd=watershed_mfd,
 )
 
 hand = _RoutingDispatch(
     'hand',
+    summary="Height above the nearest drainage (HAND).",
     d8=hand_d8, dinf=hand_dinf, mfd=hand_mfd,
 )
 
 stream_link = _RoutingDispatch(
     'stream_link',
+    summary="Assign unique IDs to stream segments above a flow-accumulation "
+            "threshold.",
     d8=stream_link_d8, dinf=stream_link_dinf, mfd=stream_link_mfd,
 )
 
@@ -135,15 +153,41 @@ def __call__(self, *args, routing='d8', ordering='strahler', **kwargs):
 
 stream_order = _StreamOrderDispatch(
     'stream_order',
+    summary="Strahler or Shreve stream ordering of the stream network.",
     d8=stream_order_d8, dinf=stream_order_dinf, mfd=stream_order_mfd,
 )
 
 
 # -- 5 D8-only functions (future-proofed with routing param) --------------
 
-basin = _RoutingDispatch('basin', d8=basin_d8)
-basins = _RoutingDispatch('basins', d8=basins_d8)
-sink = _RoutingDispatch('sink', d8=sink_d8)
-snap_pour_point = _RoutingDispatch('snap_pour_point', d8=snap_pour_point_d8)
-fill = _RoutingDispatch('fill', d8=fill_d8)
-twi = _RoutingDispatch('twi', d8=twi_d8)
+basin = _RoutingDispatch(
+    'basin',
+    summary="Label each cell with the outlet of the basin it drains to.",
+    d8=basin_d8,
+)
+basins = _RoutingDispatch(
+    'basins',
+    summary="Deprecated alias for :func:`basin`.",
+    d8=basins_d8,
+)
+sink = _RoutingDispatch(
+    'sink',
+    summary="Label depression cells that have no downstream neighbor.",
+    d8=sink_d8,
+)
+snap_pour_point = _RoutingDispatch(
+    'snap_pour_point',
+    summary="Snap pour points to the nearby cell of highest flow "
+            "accumulation.",
+    d8=snap_pour_point_d8,
+)
+fill = _RoutingDispatch(
+    'fill',
+    summary="Fill depressions in a DEM so every cell can drain.",
+    d8=fill_d8,
+)
+twi = _RoutingDispatch(
+    'twi',
+    summary="Topographic wetness index, ln(a / tan(beta)).",
+    d8=twi_d8,
+)
diff --git a/xrspatial/hydro/tests/test_routing_public_api.py b/xrspatial/hydro/tests/test_routing_public_api.py
new file mode 100644
index 000000000..de3962ba7
--- /dev/null
+++ b/xrspatial/hydro/tests/test_routing_public_api.py
@@ -0,0 +1,100 @@
+"""Public hydrology API surface: routing wrappers only (#3528).
+
+The three routing flavors (d8, dinf, mfd) are selected with the ``routing``
+keyword on a single public wrapper per family. The suffixed implementations
+(``flow_direction_d8`` etc.) are internal to ``xrspatial.hydro`` and are not
+exported at the top level.
+"""
+
+import numpy as np
+import pytest
+
+import xrspatial
+from xrspatial import hydro
+from xrspatial.tests.general_checks import create_test_raster
+
+# One public wrapper per family.
+WRAPPERS = [
+    'flow_direction', 'flow_accumulation', 'flow_length', 'flow_path',
+    'watershed', 'hand', 'stream_link', 'stream_order',
+    'basin', 'basins', 'sink', 'snap_pour_point', 'fill', 'twi',
+]
+
+# Representative suffixed implementations across every family.
+SUFFIXED = [
+    'flow_direction_d8', 'flow_direction_dinf', 'flow_direction_mfd',
+    'flow_accumulation_d8', 'flow_accumulation_dinf', 'flow_accumulation_mfd',
+    'flow_length_d8', 'flow_path_mfd', 'watershed_dinf', 'hand_mfd',
+    'stream_link_d8', 'stream_order_dinf',
+    'basin_d8', 'basins_d8', 'sink_d8', 'snap_pour_point_d8', 'fill_d8',
+    'twi_d8',
+]
+
+
+def _ramp():
+    """A plane z = i + j so dinf/mfd produce real (non-pit) flow."""
+    data = np.add.outer(np.arange(6.0), np.arange(6.0))
+    return create_test_raster(data)
+
+
+@pytest.mark.parametrize('name', WRAPPERS)
+def test_wrapper_is_public(name):
+    assert hasattr(xrspatial, name)
+    # the top-level name is the same object exported from xrspatial.hydro
+    assert getattr(xrspatial, name) is getattr(hydro, name)
+
+
+@pytest.mark.parametrize('name', SUFFIXED)
+def test_suffixed_not_public_top_level(name):
+    assert not hasattr(xrspatial, name)
+
+
+@pytest.mark.parametrize('name', SUFFIXED)
+def test_suffixed_importable_from_hydro(name):
+    # still reachable internally, where the wrappers dispatch to them
+    assert hasattr(hydro, name)
+
+
+@pytest.mark.parametrize('routing', ['d8', 'dinf', 'mfd'])
+def test_dispatch_matches_direct_impl(routing):
+    agg = _ramp()
+    direct = getattr(hydro, f'flow_direction_{routing}')
+    expected = direct(agg)
+    result = xrspatial.flow_direction(agg, routing=routing)
+    np.testing.assert_array_equal(
+        np.asarray(result.data), np.asarray(expected.data)
+    )
+
+
+def test_default_routing_is_d8():
+    agg = _ramp()
+    default = xrspatial.flow_direction(agg)
+    d8 = xrspatial.flow_direction(agg, routing='d8')
+    np.testing.assert_array_equal(
+        np.asarray(default.data), np.asarray(d8.data)
+    )
+
+
+def test_unknown_routing_raises():
+    agg = _ramp()
+    with pytest.raises(ValueError, match="Unknown routing"):
+        xrspatial.flow_direction(agg, routing='bogus')
+
+
+def test_stream_order_threads_ordering():
+    # the d8 dispatch must forward `ordering` to stream_order_d8
+    agg = _ramp()
+    fdir = xrspatial.flow_direction(agg, routing='d8')
+    accum = xrspatial.flow_accumulation(fdir, routing='d8')
+    via_wrapper = xrspatial.stream_order(
+        fdir, accum, threshold=1, ordering='shreve'
+    )
+    direct = hydro.stream_order_d8(fdir, accum, threshold=1, ordering='shreve')
+    np.testing.assert_array_equal(
+        np.asarray(via_wrapper.data), np.asarray(direct.data)
+    )
+
+
+def test_wrapper_is_self_describing():
+    assert xrspatial.flow_direction.__name__ == 'flow_direction'
+    assert 'routing' in xrspatial.flow_direction.__doc__
diff --git a/xrspatial/tests/test_dask_task_names.py b/xrspatial/tests/test_dask_task_names.py
index 2b82203b2..61eda3e2c 100644
--- a/xrspatial/tests/test_dask_task_names.py
+++ b/xrspatial/tests/test_dask_task_names.py
@@ -9,8 +9,6 @@
     aspect,
     bilateral,
     curvature,
-    flow_direction_mfd,
-    flow_length_mfd,
     generate_terrain,
     hillshade,
     preview,
@@ -19,6 +17,7 @@
     surface_direction,
     surface_distance,
 )
+from xrspatial.hydro import flow_direction_mfd, flow_length_mfd
 from xrspatial.convolution import convolve_2d
 from xrspatial.tests.general_checks import create_test_raster
 

From 50c116952a464df64aaf5f78026539197ae4a0d3 Mon Sep 17 00:00:00 2001
From: Brendan Collins <brendancol@gmail.com>
Date: Thu, 25 Jun 2026 13:00:14 -0700
Subject: [PATCH 2/4] Document hydrology routing wrappers; collapse README
 matrix (#3528)

Restructure hydrology.rst so each family leads with the unified wrapper and
its routing parameter, listing the d8/dinf/mfd implementations as routing
variants. Collapse the README feature matrix from 29 per-variant rows to 13
family rows, noting the routing options.
---
 README.md                           |  42 ++----
 docs/source/reference/hydrology.rst | 214 +++++++++++-----------------
 2 files changed, 97 insertions(+), 159 deletions(-)

diff --git a/README.md b/README.md
index 597377c22..8dc957cc6 100644
--- a/README.md
+++ b/README.md
@@ -310,35 +310,19 @@ Built-in Numba JIT and CUDA projection kernels bypass pyproj for per-pixel coord
 
 | Name | Description | Source | NumPy xr.DataArray | Dask xr.DataArray | CuPy GPU xr.DataArray | Dask GPU xr.DataArray |
 |:----------:|:------------|:------:|:----------------------:|:--------------------:|:-------------------:|:------:|
-| [Flow Direction (D8)](xrspatial/hydro/flow_direction_d8.py) | Computes D8 flow direction from each cell toward the steepest downhill neighbor | O'Callaghan & Mark 1984 | ✅ | 🔼 | 🔼 | 🔼 |
-| [Flow Direction (Dinf)](xrspatial/hydro/flow_direction_dinf.py) | Computes D-infinity flow direction as a continuous angle toward the steepest downslope facet | Tarboton 1997 | ✅ | 🔼 | 🔼 | 🔼 |
-| [Flow Direction (MFD)](xrspatial/hydro/flow_direction_mfd.py) | Partitions flow to all downslope neighbors with an adaptive exponent (Qin et al. 2007) | Qin et al. 2007 | ✅ | 🔼 | 🔼 | 🔼 |
-| [Flow Accumulation (D8)](xrspatial/hydro/flow_accumulation_d8.py) | Counts upstream cells draining through each cell in a D8 flow direction grid | Jenson & Domingue 1988 | ✅ | 🔼 | 🔼 | 🔼 |
-| [Flow Accumulation (Dinf)](xrspatial/hydro/flow_accumulation_dinf.py) | Accumulates upstream area by splitting flow proportionally between two neighbors (Tarboton 1997) | Tarboton 1997 | ✅ | 🔼 | 🔼 | 🔼 |
-| [Flow Accumulation (MFD)](xrspatial/hydro/flow_accumulation_mfd.py) | Accumulates upstream area through all MFD flow paths weighted by directional fractions | Qin et al. 2007 | ✅ | 🔼 | 🔼 | 🔼 |
-| [Flow Length (D8)](xrspatial/hydro/flow_length_d8.py) | Computes D8 flow path length from each cell to outlet (downstream) or from divide (upstream) | Standard (D8 tracing) | 🔼 | 🔼 | 🔼 | 🔼 |
-| [Flow Length (Dinf)](xrspatial/hydro/flow_length_dinf.py) | Proportion-weighted flow path length using D-inf angle decomposition (downstream or upstream) | Tarboton 1997 | 🔼 | 🔼 | 🔼 | 🔼 |
-| [Flow Length (MFD)](xrspatial/hydro/flow_length_mfd.py) | Proportion-weighted flow path length using MFD fractions (downstream or upstream) | Qin et al. 2007 | 🔼 | 🔼 | 🔼 | 🔼 |
-| [Fill (D8)](xrspatial/hydro/fill_d8.py) | Fills depressions in a DEM using Planchon-Darboux iterative flooding | Planchon & Darboux 2002 | ✅ | 🔼 | 🔼 | 🔼 |
-| [Sink (D8)](xrspatial/hydro/sink_d8.py) | Identifies and labels depression cells in a D8 flow direction grid | Standard (D8 tracing) | ✅ | 🔼 | 🔼 | 🔼 |
-| [Watershed (D8)](xrspatial/hydro/watershed_d8.py) | Labels each cell with the pour point it drains to via D8 flow direction | Standard (D8 tracing) | ✅ | 🔼 | 🔼 | 🔼 |
-| [Watershed (Dinf)](xrspatial/hydro/watershed_dinf.py) | Labels each cell with the pour point it drains to via D-infinity flow direction | Tarboton 1997 | ✅ | 🔼 | 🔼 | 🔼 |
-| [Watershed (MFD)](xrspatial/hydro/watershed_mfd.py) | Labels each cell with the pour point it drains to via MFD fractions | Qin et al. 2007 | ✅ | 🔼 | 🔼 | 🔼 |
-| [Basins](xrspatial/hydro/watershed_d8.py) | Delineates drainage basins by labeling each cell with its outlet ID | Standard (D8 tracing) | ✅ | 🔼 | 🔼 | 🔼 |
-| [Stream Order (D8)](xrspatial/hydro/stream_order_d8.py) | Assigns Strahler or Shreve stream order to cells in a drainage network | Strahler 1957, Shreve 1966 | ✅ | 🔼 | 🔼 | 🔼 |
-| [Stream Order (Dinf)](xrspatial/hydro/stream_order_dinf.py) | Strahler/Shreve stream ordering on D-infinity flow direction grids | Tarboton 1997 | ✅ | 🔼 | 🔼 | 🔼 |
-| [Stream Order (MFD)](xrspatial/hydro/stream_order_mfd.py) | Strahler/Shreve stream ordering on MFD fraction grids | Freeman 1991 | ✅ | 🔼 | 🔼 | 🔼 |
-| [Stream Link (D8)](xrspatial/hydro/stream_link_d8.py) | Assigns unique IDs to each stream segment between junctions | Standard | ✅ | 🔼 | 🔼 | 🔼 |
-| [Stream Link (Dinf)](xrspatial/hydro/stream_link_dinf.py) | Stream link segmentation on D-infinity flow direction grids | Tarboton 1997 | ✅ | 🔼 | 🔼 | 🔼 |
-| [Stream Link (MFD)](xrspatial/hydro/stream_link_mfd.py) | Stream link segmentation on MFD fraction grids | Freeman 1991 | ✅ | 🔼 | 🔼 | 🔼 |
-| [Snap Pour Point](xrspatial/hydro/snap_pour_point_d8.py) | Snaps pour points to the highest-accumulation cell within a search radius | Custom | ✅ | 🔼 | 🔼 | 🔼 |
-| [Flow Path (D8)](xrspatial/hydro/flow_path_d8.py) | Traces downstream flow paths from start points through a D8 direction grid | Standard (D8 tracing) | ✅ | 🔼 | 🔼 | 🔼 |
-| [Flow Path (Dinf)](xrspatial/hydro/flow_path_dinf.py) | Traces downstream flow paths using D-infinity dominant neighbor | Tarboton 1997 | 🔼 | 🔼 | 🔼 | 🔼 |
-| [Flow Path (MFD)](xrspatial/hydro/flow_path_mfd.py) | Traces downstream flow paths through MFD fraction-weighted neighbors | Qin et al. 2007 | 🔼 | 🔼 | 🔼 | 🔼 |
-| [HAND (D8)](xrspatial/hydro/hand_d8.py) | Computes Height Above Nearest Drainage by tracing D8 flow to the nearest stream cell | Nobre et al. 2011 | ✅ | 🔼 | 🔼 | 🔼 |
-| [HAND (Dinf)](xrspatial/hydro/hand_dinf.py) | Computes Height Above Nearest Drainage using D-infinity flow direction | Nobre et al. 2011 | ✅ | 🔼 | 🔼 | 🔼 |
-| [HAND (MFD)](xrspatial/hydro/hand_mfd.py) | Computes Height Above Nearest Drainage using MFD fractions | Nobre et al. 2011 | ✅ | 🔼 | 🔼 | 🔼 |
-| [TWI](xrspatial/hydro/twi_d8.py) | Topographic Wetness Index: ln(specific catchment area / tan(slope)) | Beven & Kirkby 1979 | ✅ | 🔼 | 🔼 | 🔼 |
+| [Flow Direction](xrspatial/hydro/flow_direction_d8.py) | Direction of steepest descent out of each cell (D8 · Dinf · MFD via `routing=`) | O'Callaghan & Mark 1984; Tarboton 1997; Qin et al. 2007 | ✅ | 🔼 | 🔼 | 🔼 |
+| [Flow Accumulation](xrspatial/hydro/flow_accumulation_d8.py) | Upstream cells or area draining through each cell (D8 · Dinf · MFD via `routing=`) | Jenson & Domingue 1988; Tarboton 1997; Qin et al. 2007 | ✅ | 🔼 | 🔼 | 🔼 |
+| [Flow Length](xrspatial/hydro/flow_length_d8.py) | Flow path length to the outlet or from the divide (D8 · Dinf · MFD via `routing=`) | Tarboton 1997; Qin et al. 2007 | 🔼 | 🔼 | 🔼 | 🔼 |
+| [Flow Path](xrspatial/hydro/flow_path_d8.py) | Traces downstream flow paths from start points (D8 · Dinf · MFD via `routing=`) | Tarboton 1997; Qin et al. 2007 | 🔼 | 🔼 | 🔼 | 🔼 |
+| [Watershed](xrspatial/hydro/watershed_d8.py) | Labels each cell with the pour point it drains to (D8 · Dinf · MFD via `routing=`) | Tarboton 1997; Qin et al. 2007 | ✅ | 🔼 | 🔼 | 🔼 |
+| [Stream Link](xrspatial/hydro/stream_link_d8.py) | Assigns unique IDs to stream segments above a threshold (D8 · Dinf · MFD via `routing=`) | Tarboton 1997; Freeman 1991 | ✅ | 🔼 | 🔼 | 🔼 |
+| [Stream Order](xrspatial/hydro/stream_order_d8.py) | Strahler or Shreve stream ordering of the network (D8 · Dinf · MFD via `routing=`) | Strahler 1957, Shreve 1966 | ✅ | 🔼 | 🔼 | 🔼 |
+| [HAND](xrspatial/hydro/hand_d8.py) | Height Above Nearest Drainage (D8 · Dinf · MFD via `routing=`) | Nobre et al. 2011 | ✅ | 🔼 | 🔼 | 🔼 |
+| [Fill](xrspatial/hydro/fill_d8.py) | Fills depressions in a DEM using Planchon-Darboux iterative flooding (D8) | Planchon & Darboux 2002 | ✅ | 🔼 | 🔼 | 🔼 |
+| [Sink](xrspatial/hydro/sink_d8.py) | Identifies and labels depression cells (D8) | Standard (D8 tracing) | ✅ | 🔼 | 🔼 | 🔼 |
+| [Basin](xrspatial/hydro/basin_d8.py) | Labels each cell with the outlet of the basin it drains to (D8) | Standard (D8 tracing) | ✅ | 🔼 | 🔼 | 🔼 |
+| [Snap Pour Point](xrspatial/hydro/snap_pour_point_d8.py) | Snaps pour points to the highest-accumulation cell within a search radius (D8) | Custom | ✅ | 🔼 | 🔼 | 🔼 |
+| [TWI](xrspatial/hydro/twi_d8.py) | Topographic Wetness Index: ln(specific catchment area / tan(slope)) (D8) | Beven & Kirkby 1979 | ✅ | 🔼 | 🔼 | 🔼 |
 
 -----------
 
diff --git a/docs/source/reference/hydrology.rst b/docs/source/reference/hydrology.rst
index 3d6a512dc..a053def8b 100644
--- a/docs/source/reference/hydrology.rst
+++ b/docs/source/reference/hydrology.rst
@@ -11,206 +11,160 @@ Hydrology
    masked out), fill or interpolate them first, or expect disconnected
    drainage networks.
 
-Flow Direction (D8)
-===================
-.. autosummary::
-    :toctree: _autosummary
+Each family exposes a single public function. The routing algorithm
+(``'d8'``, ``'dinf'``, or ``'mfd'``) is chosen with the ``routing`` keyword,
+which defaults to ``'d8'``::
 
-    xrspatial.hydro.flow_direction_d8.flow_direction_d8
+    import xrspatial
+    fdir = xrspatial.flow_direction(dem, routing='dinf')
+    acc = xrspatial.flow_accumulation(fdir, routing='dinf')
 
-Flow Direction (D-infinity)
-===========================
-.. autosummary::
-    :toctree: _autosummary
+The wrapper dispatches to the per-routing implementations listed under each
+family. Those implementations live in ``xrspatial.hydro`` and carry the full
+parameter list and algorithm references.
 
-    xrspatial.hydro.flow_direction_dinf.flow_direction_dinf
+Flow Direction
+==============
+.. py:function:: xrspatial.flow_direction(agg, *, routing='d8', **kwargs)
 
-Flow Direction (MFD)
-====================
-.. autosummary::
-    :toctree: _autosummary
+   Direction of steepest descent out of each cell.  ``routing`` selects
+   ``'d8'``, ``'dinf'``, or ``'mfd'`` (default ``'d8'``).
 
-    xrspatial.hydro.flow_direction_mfd.flow_direction_mfd
+Routing variants:
 
-Flow Accumulation (D8)
-======================
 .. autosummary::
     :toctree: _autosummary
 
-    xrspatial.hydro.flow_accumulation_d8.flow_accumulation_d8
+    xrspatial.hydro.flow_direction_d8.flow_direction_d8
+    xrspatial.hydro.flow_direction_dinf.flow_direction_dinf
+    xrspatial.hydro.flow_direction_mfd.flow_direction_mfd
 
-Flow Accumulation (D-infinity)
-===============================
-.. autosummary::
-    :toctree: _autosummary
+Flow Accumulation
+=================
+.. py:function:: xrspatial.flow_accumulation(flow_dir, *, routing='d8', **kwargs)
 
-    xrspatial.hydro.flow_accumulation_dinf.flow_accumulation_dinf
+   Upstream cells or area draining through each cell.  ``routing`` selects
+   ``'d8'``, ``'dinf'``, or ``'mfd'`` (default ``'d8'``).
+
+Routing variants:
 
-Flow Accumulation (MFD)
-=======================
 .. autosummary::
     :toctree: _autosummary
 
+    xrspatial.hydro.flow_accumulation_d8.flow_accumulation_d8
+    xrspatial.hydro.flow_accumulation_dinf.flow_accumulation_dinf
     xrspatial.hydro.flow_accumulation_mfd.flow_accumulation_mfd
 
-Flow Length (D8)
-================
-.. autosummary::
-    :toctree: _autosummary
-
-    xrspatial.hydro.flow_length_d8.flow_length_d8
+Flow Length
+===========
+.. py:function:: xrspatial.flow_length(flow_dir, *, routing='d8', **kwargs)
 
-Flow Length (D-infinity)
-========================
-.. autosummary::
-    :toctree: _autosummary
+   Distance along the flow path to the outlet or from the divide.
+   ``routing`` selects ``'d8'``, ``'dinf'``, or ``'mfd'`` (default ``'d8'``).
 
-    xrspatial.hydro.flow_length_dinf.flow_length_dinf
+Routing variants:
 
-Flow Length (MFD)
-=================
 .. autosummary::
     :toctree: _autosummary
 
+    xrspatial.hydro.flow_length_d8.flow_length_d8
+    xrspatial.hydro.flow_length_dinf.flow_length_dinf
     xrspatial.hydro.flow_length_mfd.flow_length_mfd
 
-Flow Path (D8)
-==============
-.. autosummary::
-    :toctree: _autosummary
-
-    xrspatial.hydro.flow_path_d8.flow_path_d8
+Flow Path
+=========
+.. py:function:: xrspatial.flow_path(flow_dir, start_points, *, routing='d8', **kwargs)
 
-Flow Path (D-infinity)
-======================
-.. autosummary::
-    :toctree: _autosummary
+   Trace downstream flow paths from a set of start points.  ``routing``
+   selects ``'d8'``, ``'dinf'``, or ``'mfd'`` (default ``'d8'``).
 
-    xrspatial.hydro.flow_path_dinf.flow_path_dinf
+Routing variants:
 
-Flow Path (MFD)
-===============
 .. autosummary::
     :toctree: _autosummary
 
+    xrspatial.hydro.flow_path_d8.flow_path_d8
+    xrspatial.hydro.flow_path_dinf.flow_path_dinf
     xrspatial.hydro.flow_path_mfd.flow_path_mfd
 
-Fill
-====
-.. autosummary::
-    :toctree: _autosummary
+Watershed
+=========
+.. py:function:: xrspatial.watershed(flow_dir, pour_points, *, routing='d8', **kwargs)
 
-    xrspatial.hydro.fill_d8.fill_d8
+   Label each cell with the pour point it drains to.  ``routing`` selects
+   ``'d8'``, ``'dinf'``, or ``'mfd'`` (default ``'d8'``).
 
-Sink
-====
-.. autosummary::
-    :toctree: _autosummary
+Routing variants:
 
-    xrspatial.hydro.sink_d8.sink_d8
-
-Basin
-=====
-.. autosummary::
-    :toctree: _autosummary
-
-    xrspatial.hydro.basin_d8.basin_d8
-
-Watershed (D8)
-==============
 .. autosummary::
     :toctree: _autosummary
 
     xrspatial.hydro.watershed_d8.watershed_d8
-    xrspatial.hydro.watershed_d8.basins_d8
-
-Watershed (D-infinity)
-======================
-.. autosummary::
-    :toctree: _autosummary
-
     xrspatial.hydro.watershed_dinf.watershed_dinf
-
-Watershed (MFD)
-===============
-.. autosummary::
-    :toctree: _autosummary
-
     xrspatial.hydro.watershed_mfd.watershed_mfd
+    xrspatial.hydro.watershed_d8.basins_d8
 
-Snap Pour Point
-===============
-.. autosummary::
-    :toctree: _autosummary
-
-    xrspatial.hydro.snap_pour_point_d8.snap_pour_point_d8
+Stream Link
+===========
+.. py:function:: xrspatial.stream_link(flow_dir, flow_accum, *, routing='d8', threshold=100, **kwargs)
 
-Stream Link (D8)
-================
-.. autosummary::
-    :toctree: _autosummary
+   Assign unique IDs to stream segments above a flow-accumulation threshold.
+   ``routing`` selects ``'d8'``, ``'dinf'``, or ``'mfd'`` (default ``'d8'``).
 
-    xrspatial.hydro.stream_link_d8.stream_link_d8
+Routing variants:
 
-Stream Link (D-infinity)
-========================
 .. autosummary::
     :toctree: _autosummary
 
+    xrspatial.hydro.stream_link_d8.stream_link_d8
     xrspatial.hydro.stream_link_dinf.stream_link_dinf
-
-Stream Link (MFD)
-=================
-.. autosummary::
-    :toctree: _autosummary
-
     xrspatial.hydro.stream_link_mfd.stream_link_mfd
 
-Stream Order (D8)
-=================
-.. autosummary::
-    :toctree: _autosummary
+Stream Order
+============
+.. py:function:: xrspatial.stream_order(flow_dir, flow_accum, *, routing='d8', ordering='strahler', threshold=100, **kwargs)
 
-    xrspatial.hydro.stream_order_d8.stream_order_d8
+   Strahler or Shreve stream ordering of the stream network.  ``routing``
+   selects ``'d8'``, ``'dinf'``, or ``'mfd'`` (default ``'d8'``); ``ordering``
+   selects ``'strahler'`` or ``'shreve'``.
 
-Stream Order (D-infinity)
-=========================
-.. autosummary::
-    :toctree: _autosummary
-
-    xrspatial.hydro.stream_order_dinf.stream_order_dinf
+Routing variants:
 
-Stream Order (MFD)
-==================
 .. autosummary::
     :toctree: _autosummary
 
+    xrspatial.hydro.stream_order_d8.stream_order_d8
+    xrspatial.hydro.stream_order_dinf.stream_order_dinf
     xrspatial.hydro.stream_order_mfd.stream_order_mfd
 
-Height Above Nearest Drainage (D8)
-===================================
-.. autosummary::
-    :toctree: _autosummary
+Height Above Nearest Drainage (HAND)
+====================================
+.. py:function:: xrspatial.hand(flow_dir, flow_accum, elevation, *, routing='d8', threshold=100, **kwargs)
 
-    xrspatial.hydro.hand_d8.hand_d8
+   Height above the nearest drainage.  ``routing`` selects ``'d8'``,
+   ``'dinf'``, or ``'mfd'`` (default ``'d8'``).
+
+Routing variants:
 
-Height Above Nearest Drainage (D-infinity)
-==========================================
 .. autosummary::
     :toctree: _autosummary
 
+    xrspatial.hydro.hand_d8.hand_d8
     xrspatial.hydro.hand_dinf.hand_dinf
+    xrspatial.hydro.hand_mfd.hand_mfd
 
-Height Above Nearest Drainage (MFD)
-====================================
-.. autosummary::
-    :toctree: _autosummary
+D8-only functions
+=================
 
-    xrspatial.hydro.hand_mfd.hand_mfd
+These families currently implement D8 routing only. They take the same
+``routing`` keyword for forward compatibility, where ``'d8'`` is the only
+accepted value.
 
-Topographic Wetness Index (TWI)
-===============================
 .. autosummary::
     :toctree: _autosummary
 
+    xrspatial.hydro.fill_d8.fill_d8
+    xrspatial.hydro.sink_d8.sink_d8
+    xrspatial.hydro.basin_d8.basin_d8
+    xrspatial.hydro.snap_pour_point_d8.snap_pour_point_d8
     xrspatial.hydro.twi_d8.twi_d8

From bd495ecef6e5290df65afdbbd365a31812038756 Mon Sep 17 00:00:00 2001
From: Brendan Collins <brendancol@gmail.com>
Date: Thu, 25 Jun 2026 13:08:57 -0700
Subject: [PATCH 3/4] Update example notebooks to routing= API (#3528)

Address review blocker: three user-guide notebooks imported/called the
suffixed hydrology names from top-level xrspatial, which this PR removes.
Switch them to the wrapper + routing= form (flow_direction(routing='mfd'),
flow_accumulation(routing='mfd'), flow_length(routing='mfd', ...),
stream_order(routing='dinf', ordering=...), stream_link(routing='mfd', ...)).
Also convert the d8 stream_order calls from the fragile method= kwarg to
ordering= for consistency, and update two prose cells that named the removed
functions.
---
 examples/user_guide/11_Hydrology.ipynb        |  85 +-----------
 .../27_Stream_Analysis_Dinf_MFD.ipynb         | 126 +-----------------
 examples/user_guide/32_Flow_Length_MFD.ipynb  |  43 +-----
 3 files changed, 17 insertions(+), 237 deletions(-)

diff --git a/examples/user_guide/11_Hydrology.ipynb b/examples/user_guide/11_Hydrology.ipynb
index 13202463b..1812e21a4 100644
--- a/examples/user_guide/11_Hydrology.ipynb
+++ b/examples/user_guide/11_Hydrology.ipynb
@@ -319,13 +319,7 @@
    "cell_type": "markdown",
    "id": "nbozw06cw1",
    "metadata": {},
-   "source": [
-    "## Multiple Flow Direction\n",
-    "\n",
-    "D8 sends all flow to a single neighbor. The `flow_direction_mfd` function partitions flow to all downslope neighbors, using an adaptive exponent from [Qin et al. (2007)](https://doi.org/10.1080/13658810601168578) that concentrates flow on steep slopes and disperses it on gentle ones. The output is a 3-D DataArray `(8, H, W)` with fractional weights that sum to 1.0 at each cell.\n",
-    "\n",
-    "The plot below shows the maximum fraction per cell: values near 1.0 mean flow is concentrated (D8-like), while lower values mean it spreads across multiple neighbors. See [notebook 27](27_Stream_Analysis_Dinf_MFD.ipynb) for a deeper comparison of D8, D-infinity, and MFD routing."
-   ]
+   "source": "## Multiple Flow Direction\n\nD8 sends all flow to a single neighbor. Passing `routing='mfd'` to `flow_direction` partitions flow to all downslope neighbors, using an adaptive exponent from [Qin et al. (2007)](https://doi.org/10.1080/13658810601168578) that concentrates flow on steep slopes and disperses it on gentle ones. The output is a 3-D DataArray `(8, H, W)` with fractional weights that sum to 1.0 at each cell.\n\nThe plot below shows the maximum fraction per cell: values near 1.0 mean flow is concentrated (D8-like), while lower values mean it spreads across multiple neighbors. See [notebook 27](27_Stream_Analysis_Dinf_MFD.ipynb) for a deeper comparison of D8, D-infinity, and MFD routing."
   },
   {
    "cell_type": "code",
@@ -340,20 +334,7 @@
     }
    },
    "outputs": [],
-   "source": [
-    "mfd_fracs = xrspatial.flow_direction_mfd(dem_filled)\n",
-    "print(f\"MFD output shape: {mfd_fracs.shape}\")\n",
-    "print(f\"Dims: {mfd_fracs.dims}\")\n",
-    "\n",
-    "max_frac = mfd_fracs.max(dim='neighbor')\n",
-    "\n",
-    "fig, ax = plt.subplots(figsize=(10, 7.5))\n",
-    "plot_basemap(ax)\n",
-    "max_frac.plot.imshow(ax=ax, cmap='YlOrRd', alpha=0.8, add_colorbar=True,\n",
-    "                     cbar_kwargs={'label': 'Max flow fraction (1.0 = single receiver)',\n",
-    "                                  'shrink': 0.7})\n",
-    "ax.set_axis_off()"
-   ]
+   "source": "mfd_fracs = xrspatial.flow_direction(dem_filled, routing='mfd')\nprint(f\"MFD output shape: {mfd_fracs.shape}\")\nprint(f\"Dims: {mfd_fracs.dims}\")\n\nmax_frac = mfd_fracs.max(dim='neighbor')\n\nfig, ax = plt.subplots(figsize=(10, 7.5))\nplot_basemap(ax)\nmax_frac.plot.imshow(ax=ax, cmap='YlOrRd', alpha=0.8, add_colorbar=True,\n                     cbar_kwargs={'label': 'Max flow fraction (1.0 = single receiver)',\n                                  'shrink': 0.7})\nax.set_axis_off()"
   },
   {
    "cell_type": "markdown",
@@ -439,13 +420,7 @@
    "cell_type": "markdown",
    "id": "httgn03vcco",
    "metadata": {},
-   "source": [
-    "## Flow Accumulation (MFD)\n",
-    "\n",
-    "`flow_accumulation_mfd` routes upstream contributing area through all downslope paths at once, using the fractional weights from `flow_direction_mfd`. Where D8 accumulation produces sharp single-pixel drainage lines, MFD accumulation spreads flow and produces smoother contributing-area fields.\n",
-    "\n",
-    "The side-by-side comparison shows D8 (left) with crisp channels and MFD (right) with diffuse drainage patterns."
-   ]
+   "source": "## Flow Accumulation (MFD)\n\n`flow_accumulation` with `routing='mfd'` routes upstream contributing area through all downslope paths at once, using the fractional weights from `flow_direction(routing='mfd')`. Where D8 accumulation produces sharp single-pixel drainage lines, MFD accumulation spreads flow and produces smoother contributing-area fields.\n\nThe side-by-side comparison shows D8 (left) with crisp channels and MFD (right) with diffuse drainage patterns."
   },
   {
    "cell_type": "code",
@@ -460,28 +435,7 @@
     }
    },
    "outputs": [],
-   "source": [
-    "flow_accum_mfd = xrspatial.flow_accumulation_mfd(mfd_fracs)\n",
-    "\n",
-    "# Log-transform for visualization (log1p handles zeros)\n",
-    "log_d8 = xr.DataArray(np.log1p(flow_accum.values),\n",
-    "                      dims=flow_accum.dims, coords=flow_accum.coords)\n",
-    "log_mfd = xr.DataArray(np.log1p(flow_accum_mfd.values),\n",
-    "                       dims=flow_accum_mfd.dims, coords=flow_accum_mfd.coords)\n",
-    "\n",
-    "fig, axes = plt.subplots(1, 2, figsize=(14, 5))\n",
-    "fig.patch.set_facecolor('white')\n",
-    "\n",
-    "for ax, data, title in zip(axes, [log_d8, log_mfd],\n",
-    "                            ['D8 flow accumulation', 'MFD flow accumulation']):\n",
-    "    ax.set_facecolor('white')\n",
-    "    data.plot.imshow(ax=ax, cmap='Blues', add_colorbar=True,\n",
-    "                     cbar_kwargs={'label': 'log(1 + accum)', 'shrink': 0.7})\n",
-    "    ax.set_title(title)\n",
-    "    ax.set_axis_off()\n",
-    "\n",
-    "plt.tight_layout()"
-   ]
+   "source": "flow_accum_mfd = xrspatial.flow_accumulation(mfd_fracs, routing='mfd')\n\n# Log-transform for visualization (log1p handles zeros)\nlog_d8 = xr.DataArray(np.log1p(flow_accum.values),\n                      dims=flow_accum.dims, coords=flow_accum.coords)\nlog_mfd = xr.DataArray(np.log1p(flow_accum_mfd.values),\n                       dims=flow_accum_mfd.dims, coords=flow_accum_mfd.coords)\n\nfig, axes = plt.subplots(1, 2, figsize=(14, 5))\nfig.patch.set_facecolor('white')\n\nfor ax, data, title in zip(axes, [log_d8, log_mfd],\n                            ['D8 flow accumulation', 'MFD flow accumulation']):\n    ax.set_facecolor('white')\n    data.plot.imshow(ax=ax, cmap='Blues', add_colorbar=True,\n                     cbar_kwargs={'label': 'log(1 + accum)', 'shrink': 0.7})\n    ax.set_title(title)\n    ax.set_axis_off()\n\nplt.tight_layout()"
   },
   {
    "cell_type": "markdown",
@@ -554,20 +508,7 @@
     }
    },
    "outputs": [],
-   "source": [
-    "strahler = xrspatial.stream_order(\n",
-    "    flow_dir, flow_accum, threshold=200, method='strahler'\n",
-    ")\n",
-    "max_order = int(np.nanmax(strahler.values))\n",
-    "print(f\"Max Strahler order: {max_order}\")\n",
-    "\n",
-    "fig, ax = plt.subplots(figsize=(10, 7.5))\n",
-    "plot_basemap(ax)\n",
-    "strahler.plot.imshow(ax=ax, cmap=stream_cmap, alpha=0.95,\n",
-    "                     add_colorbar=True,\n",
-    "                     cbar_kwargs={'label': 'Strahler order', 'shrink': 0.7})\n",
-    "ax.set_axis_off()"
-   ]
+   "source": "strahler = xrspatial.stream_order(\n    flow_dir, flow_accum, threshold=200, ordering='strahler'\n)\nmax_order = int(np.nanmax(strahler.values))\nprint(f\"Max Strahler order: {max_order}\")\n\nfig, ax = plt.subplots(figsize=(10, 7.5))\nplot_basemap(ax)\nstrahler.plot.imshow(ax=ax, cmap=stream_cmap, alpha=0.95,\n                     add_colorbar=True,\n                     cbar_kwargs={'label': 'Strahler order', 'shrink': 0.7})\nax.set_axis_off()"
   },
   {
    "cell_type": "code",
@@ -582,19 +523,7 @@
     }
    },
    "outputs": [],
-   "source": [
-    "shreve = xrspatial.stream_order(\n",
-    "    flow_dir, flow_accum, threshold=200, method='shreve'\n",
-    ")\n",
-    "print(f\"Max Shreve magnitude: {int(np.nanmax(shreve.values))}\")\n",
-    "\n",
-    "fig, ax = plt.subplots(figsize=(10, 7.5))\n",
-    "plot_basemap(ax)\n",
-    "shreve.plot.imshow(ax=ax, cmap=stream_cmap, alpha=0.95,\n",
-    "                   add_colorbar=True,\n",
-    "                   cbar_kwargs={'label': 'Shreve magnitude', 'shrink': 0.7})\n",
-    "ax.set_axis_off()"
-   ]
+   "source": "shreve = xrspatial.stream_order(\n    flow_dir, flow_accum, threshold=200, ordering='shreve'\n)\nprint(f\"Max Shreve magnitude: {int(np.nanmax(shreve.values))}\")\n\nfig, ax = plt.subplots(figsize=(10, 7.5))\nplot_basemap(ax)\nshreve.plot.imshow(ax=ax, cmap=stream_cmap, alpha=0.95,\n                   add_colorbar=True,\n                   cbar_kwargs={'label': 'Shreve magnitude', 'shrink': 0.7})\nax.set_axis_off()"
   },
   {
    "cell_type": "markdown",
@@ -958,4 +887,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 5
-}
+}
\ No newline at end of file
diff --git a/examples/user_guide/27_Stream_Analysis_Dinf_MFD.ipynb b/examples/user_guide/27_Stream_Analysis_Dinf_MFD.ipynb
index c176e9c27..46275355c 100644
--- a/examples/user_guide/27_Stream_Analysis_Dinf_MFD.ipynb
+++ b/examples/user_guide/27_Stream_Analysis_Dinf_MFD.ipynb
@@ -45,23 +45,7 @@
     }
    },
    "outputs": [],
-   "source": [
-    "%matplotlib inline\n",
-    "import numpy as np\n",
-    "import pandas as pd\n",
-    "import xarray as xr\n",
-    "\n",
-    "import matplotlib.pyplot as plt\n",
-    "from matplotlib.colors import LinearSegmentedColormap, ListedColormap, LogNorm\n",
-    "from matplotlib.patches import Patch\n",
-    "\n",
-    "from xrspatial import generate_terrain, fill, erode, hillshade\n",
-    "from xrspatial import flow_direction, flow_direction_dinf, flow_direction_mfd\n",
-    "from xrspatial import flow_accumulation, flow_accumulation_mfd\n",
-    "from xrspatial import stream_order, stream_link\n",
-    "from xrspatial import stream_order_dinf, stream_link_dinf\n",
-    "from xrspatial import stream_order_mfd, stream_link_mfd"
-   ]
+   "source": "%matplotlib inline\nimport numpy as np\nimport pandas as pd\nimport xarray as xr\n\nimport matplotlib.pyplot as plt\nfrom matplotlib.colors import LinearSegmentedColormap, ListedColormap, LogNorm\nfrom matplotlib.patches import Patch\n\nfrom xrspatial import generate_terrain, fill, erode, hillshade\nfrom xrspatial import flow_direction, flow_accumulation\nfrom xrspatial import stream_order, stream_link"
   },
   {
    "cell_type": "markdown",
@@ -143,37 +127,7 @@
     }
    },
    "outputs": [],
-   "source": [
-    "# D8\n",
-    "fd_d8 = flow_direction(dem)\n",
-    "fa_d8 = flow_accumulation(fd_d8)\n",
-    "\n",
-    "# D-infinity\n",
-    "fd_dinf = flow_direction_dinf(dem)\n",
-    "\n",
-    "# MFD\n",
-    "fd_mfd = flow_direction_mfd(dem)\n",
-    "fa_mfd = flow_accumulation_mfd(fd_mfd)\n",
-    "\n",
-    "print(f'D8 flow dir shape:   {fd_d8.shape}')\n",
-    "print(f'D-inf angles shape:  {fd_dinf.shape}')\n",
-    "print(f'MFD fractions shape: {fd_mfd.shape}')\n",
-    "print(f'D8 accumulation:     1 to {np.nanmax(fa_d8.values):.0f}')\n",
-    "print(f'MFD accumulation:    1 to {np.nanmax(fa_mfd.values):.0f}')\n",
-    "\n",
-    "water_cmap = LinearSegmentedColormap.from_list('water', ['#6baed6', '#08306b'])\n",
-    "\n",
-    "fig, axes = plt.subplots(1, 2, figsize=(14, 5))\n",
-    "for ax, data, title in zip(axes, [fa_d8, fa_mfd], ['D8', 'MFD']):\n",
-    "    plot_basemap(ax)\n",
-    "    data.plot.imshow(ax=ax, cmap=water_cmap, alpha=0.85,\n",
-    "                     norm=LogNorm(vmin=1, vmax=float(np.nanmax(data.values))),\n",
-    "                     add_colorbar=True,\n",
-    "                     cbar_kwargs={'label': 'Upstream cells (log)', 'shrink': 0.7})\n",
-    "    ax.set_title(f'{title} flow accumulation')\n",
-    "    ax.set_axis_off()\n",
-    "plt.tight_layout()"
-   ]
+   "source": "# D8\nfd_d8 = flow_direction(dem)\nfa_d8 = flow_accumulation(fd_d8)\n\n# D-infinity\nfd_dinf = flow_direction(dem, routing='dinf')\n\n# MFD\nfd_mfd = flow_direction(dem, routing='mfd')\nfa_mfd = flow_accumulation(fd_mfd, routing='mfd')\n\nprint(f'D8 flow dir shape:   {fd_d8.shape}')\nprint(f'D-inf angles shape:  {fd_dinf.shape}')\nprint(f'MFD fractions shape: {fd_mfd.shape}')\nprint(f'D8 accumulation:     1 to {np.nanmax(fa_d8.values):.0f}')\nprint(f'MFD accumulation:    1 to {np.nanmax(fa_mfd.values):.0f}')\n\nwater_cmap = LinearSegmentedColormap.from_list('water', ['#6baed6', '#08306b'])\n\nfig, axes = plt.subplots(1, 2, figsize=(14, 5))\nfor ax, data, title in zip(axes, [fa_d8, fa_mfd], ['D8', 'MFD']):\n    plot_basemap(ax)\n    data.plot.imshow(ax=ax, cmap=water_cmap, alpha=0.85,\n                     norm=LogNorm(vmin=1, vmax=float(np.nanmax(data.values))),\n                     add_colorbar=True,\n                     cbar_kwargs={'label': 'Upstream cells (log)', 'shrink': 0.7})\n    ax.set_title(f'{title} flow accumulation')\n    ax.set_axis_off()\nplt.tight_layout()"
   },
   {
    "cell_type": "markdown",
@@ -196,23 +150,7 @@
     }
    },
    "outputs": [],
-   "source": [
-    "threshold = 200\n",
-    "\n",
-    "# D8 stream order\n",
-    "so_d8 = stream_order(fd_d8, fa_d8, threshold=threshold, method='strahler')\n",
-    "\n",
-    "# D-inf stream order\n",
-    "so_dinf = stream_order_dinf(fd_dinf, fa_d8, threshold=threshold, method='strahler')\n",
-    "\n",
-    "# MFD stream order\n",
-    "so_mfd = stream_order_mfd(fd_mfd, fa_mfd, threshold=threshold, method='strahler')\n",
-    "\n",
-    "for label, data in [('D8', so_d8), ('D-inf', so_dinf), ('MFD', so_mfd)]:\n",
-    "    n = int(np.sum(~np.isnan(data.values)))\n",
-    "    mx = int(np.nanmax(data.values)) if n > 0 else 0\n",
-    "    print(f'{label:6s}: {n:5d} stream cells, max Strahler order {mx}')"
-   ]
+   "source": "threshold = 200\n\n# D8 stream order\nso_d8 = stream_order(fd_d8, fa_d8, threshold=threshold, ordering='strahler')\n\n# D-inf stream order\nso_dinf = stream_order(fd_dinf, fa_d8, routing='dinf', threshold=threshold, ordering='strahler')\n\n# MFD stream order\nso_mfd = stream_order(fd_mfd, fa_mfd, routing='mfd', threshold=threshold, ordering='strahler')\n\nfor label, data in [('D8', so_d8), ('D-inf', so_dinf), ('MFD', so_mfd)]:\n    n = int(np.sum(~np.isnan(data.values)))\n    mx = int(np.nanmax(data.values)) if n > 0 else 0\n    print(f'{label:6s}: {n:5d} stream cells, max Strahler order {mx}')"
   },
   {
    "cell_type": "code",
@@ -280,34 +218,7 @@
     }
    },
    "outputs": [],
-   "source": [
-    "sv_d8 = stream_order(fd_d8, fa_d8, threshold=threshold, method='shreve')\n",
-    "sv_dinf = stream_order_dinf(fd_dinf, fa_d8, threshold=threshold, method='shreve')\n",
-    "sv_mfd = stream_order_mfd(fd_mfd, fa_mfd, threshold=threshold, method='shreve')\n",
-    "\n",
-    "for label, data in [('D8', sv_d8), ('D-inf', sv_dinf), ('MFD', sv_mfd)]:\n",
-    "    n = int(np.sum(~np.isnan(data.values)))\n",
-    "    mx = int(np.nanmax(data.values)) if n > 0 else 0\n",
-    "    print(f'{label:6s}: {n:5d} stream cells, max Shreve magnitude {mx}')\n",
-    "\n",
-    "fig, axes = plt.subplots(1, 3, figsize=(18, 5.5))\n",
-    "\n",
-    "for ax, data, title in zip(axes, [sv_d8, sv_dinf, sv_mfd],\n",
-    "                            ['D8', 'D-infinity', 'MFD']):\n",
-    "    plot_basemap(ax)\n",
-    "    vmax = float(np.nanmax(data.values))\n",
-    "    display = xr.where(data.isnull(), np.nan, data)\n",
-    "    log_display = xr.DataArray(np.log1p(display.values),\n",
-    "                               dims=data.dims, coords=data.coords)\n",
-    "    log_display.plot.imshow(ax=ax, cmap='viridis', alpha=0.9,\n",
-    "                            add_colorbar=True,\n",
-    "                            cbar_kwargs={'label': 'log(1 + magnitude)',\n",
-    "                                         'shrink': 0.6})\n",
-    "    ax.set_title(f'Shreve magnitude ({title}), max={int(vmax)}')\n",
-    "    ax.set_axis_off()\n",
-    "\n",
-    "plt.tight_layout()"
-   ]
+   "source": "sv_d8 = stream_order(fd_d8, fa_d8, threshold=threshold, ordering='shreve')\nsv_dinf = stream_order(fd_dinf, fa_d8, routing='dinf', threshold=threshold, ordering='shreve')\nsv_mfd = stream_order(fd_mfd, fa_mfd, routing='mfd', threshold=threshold, ordering='shreve')\n\nfor label, data in [('D8', sv_d8), ('D-inf', sv_dinf), ('MFD', sv_mfd)]:\n    n = int(np.sum(~np.isnan(data.values)))\n    mx = int(np.nanmax(data.values)) if n > 0 else 0\n    print(f'{label:6s}: {n:5d} stream cells, max Shreve magnitude {mx}')\n\nfig, axes = plt.subplots(1, 3, figsize=(18, 5.5))\n\nfor ax, data, title in zip(axes, [sv_d8, sv_dinf, sv_mfd],\n                            ['D8', 'D-infinity', 'MFD']):\n    plot_basemap(ax)\n    vmax = float(np.nanmax(data.values))\n    display = xr.where(data.isnull(), np.nan, data)\n    log_display = xr.DataArray(np.log1p(display.values),\n                               dims=data.dims, coords=data.coords)\n    log_display.plot.imshow(ax=ax, cmap='viridis', alpha=0.9,\n                            add_colorbar=True,\n                            cbar_kwargs={'label': 'log(1 + magnitude)',\n                                         'shrink': 0.6})\n    ax.set_title(f'Shreve magnitude ({title}), max={int(vmax)}')\n    ax.set_axis_off()\n\nplt.tight_layout()"
   },
   {
    "cell_type": "markdown",
@@ -330,32 +241,7 @@
     }
    },
    "outputs": [],
-   "source": [
-    "sl_d8 = stream_link(fd_d8, fa_d8, threshold=threshold)\n",
-    "sl_dinf = stream_link_dinf(fd_dinf, fa_d8, threshold=threshold)\n",
-    "sl_mfd = stream_link_mfd(fd_mfd, fa_mfd, threshold=threshold)\n",
-    "\n",
-    "fig, axes = plt.subplots(1, 3, figsize=(18, 5.5))\n",
-    "\n",
-    "for ax, data, title in zip(axes, [sl_d8, sl_dinf, sl_mfd],\n",
-    "                            ['D8', 'D-infinity', 'MFD']):\n",
-    "    plot_basemap(ax)\n",
-    "    # Multiply by a prime before mod to break spatial color banding\n",
-    "    # (position-based IDs create vertical stripes with plain mod 20)\n",
-    "    display = xr.where(data.isnull(), np.nan, np.mod(data * 7919, 20) + 1)\n",
-    "    display.plot.imshow(ax=ax, cmap='tab20', alpha=0.9, add_colorbar=False)\n",
-    "    ax.set_title(f'Stream links ({title})')\n",
-    "    ax.set_axis_off()\n",
-    "\n",
-    "plt.tight_layout()\n",
-    "\n",
-    "for label, data in [('D8', sl_d8), ('D-inf', sl_dinf), ('MFD', sl_mfd)]:\n",
-    "    n_links = len(np.unique(data.values[~np.isnan(data.values)]))\n",
-    "    n_stream = int(np.sum(~np.isnan(data.values)))\n",
-    "    avg_len = n_stream / max(n_links, 1)\n",
-    "    print(f'{label:6s}: {n_links:4d} links, {n_stream:6d} stream cells, '\n",
-    "          f'{avg_len:.1f} cells/link avg')"
-   ]
+   "source": "sl_d8 = stream_link(fd_d8, fa_d8, threshold=threshold)\nsl_dinf = stream_link(fd_dinf, fa_d8, routing='dinf', threshold=threshold)\nsl_mfd = stream_link(fd_mfd, fa_mfd, routing='mfd', threshold=threshold)\n\nfig, axes = plt.subplots(1, 3, figsize=(18, 5.5))\n\nfor ax, data, title in zip(axes, [sl_d8, sl_dinf, sl_mfd],\n                            ['D8', 'D-infinity', 'MFD']):\n    plot_basemap(ax)\n    # Multiply by a prime before mod to break spatial color banding\n    # (position-based IDs create vertical stripes with plain mod 20)\n    display = xr.where(data.isnull(), np.nan, np.mod(data * 7919, 20) + 1)\n    display.plot.imshow(ax=ax, cmap='tab20', alpha=0.9, add_colorbar=False)\n    ax.set_title(f'Stream links ({title})')\n    ax.set_axis_off()\n\nplt.tight_layout()\n\nfor label, data in [('D8', sl_d8), ('D-inf', sl_dinf), ('MFD', sl_mfd)]:\n    n_links = len(np.unique(data.values[~np.isnan(data.values)]))\n    n_stream = int(np.sum(~np.isnan(data.values)))\n    avg_len = n_stream / max(n_links, 1)\n    print(f'{label:6s}: {n_links:4d} links, {n_stream:6d} stream cells, '\n          f'{avg_len:.1f} cells/link avg')"
   },
   {
    "cell_type": "markdown",
@@ -406,4 +292,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 4
-}
+}
\ No newline at end of file
diff --git a/examples/user_guide/32_Flow_Length_MFD.ipynb b/examples/user_guide/32_Flow_Length_MFD.ipynb
index b07f8e417..40f0566e6 100644
--- a/examples/user_guide/32_Flow_Length_MFD.ipynb
+++ b/examples/user_guide/32_Flow_Length_MFD.ipynb
@@ -20,14 +20,7 @@
    "id": "e5f6a7b8",
    "metadata": {},
    "outputs": [],
-   "source": [
-    "import numpy as np\n",
-    "import matplotlib.pyplot as plt\n",
-    "import xarray as xr\n",
-    "\n",
-    "from xrspatial import flow_direction_mfd, flow_length_mfd, flow_accumulation_mfd\n",
-    "from xrspatial.terrain import generate_terrain"
-   ]
+   "source": "import numpy as np\nimport matplotlib.pyplot as plt\nimport xarray as xr\n\nfrom xrspatial import flow_direction, flow_length, flow_accumulation\nfrom xrspatial.terrain import generate_terrain"
   },
   {
    "cell_type": "markdown",
@@ -72,17 +65,7 @@
    "id": "c1d2e3f4",
    "metadata": {},
    "outputs": [],
-   "source": [
-    "mfd_dir = flow_direction_mfd(terrain)\n",
-    "mfd_acc = flow_accumulation_mfd(mfd_dir)\n",
-    "\n",
-    "fig, ax = plt.subplots(figsize=(8, 6))\n",
-    "im = ax.imshow(np.log1p(mfd_acc.values), cmap='Blues', origin='lower')\n",
-    "ax.set_title('MFD flow accumulation (log scale)')\n",
-    "plt.colorbar(im, ax=ax, label='log(1 + accumulation)')\n",
-    "plt.tight_layout()\n",
-    "plt.show()"
-   ]
+   "source": "mfd_dir = flow_direction(terrain, routing='mfd')\nmfd_acc = flow_accumulation(mfd_dir, routing='mfd')\n\nfig, ax = plt.subplots(figsize=(8, 6))\nim = ax.imshow(np.log1p(mfd_acc.values), cmap='Blues', origin='lower')\nax.set_title('MFD flow accumulation (log scale)')\nplt.colorbar(im, ax=ax, label='log(1 + accumulation)')\nplt.tight_layout()\nplt.show()"
   },
   {
    "cell_type": "markdown",
@@ -100,16 +83,7 @@
    "id": "e9f0a1b2",
    "metadata": {},
    "outputs": [],
-   "source": [
-    "downstream = flow_length_mfd(mfd_dir, direction='downstream')\n",
-    "\n",
-    "fig, ax = plt.subplots(figsize=(8, 6))\n",
-    "im = ax.imshow(downstream.values, cmap='viridis', origin='lower')\n",
-    "ax.set_title('MFD downstream flow length')\n",
-    "plt.colorbar(im, ax=ax, label='distance (coordinate units)')\n",
-    "plt.tight_layout()\n",
-    "plt.show()"
-   ]
+   "source": "downstream = flow_length(mfd_dir, routing='mfd', direction='downstream')\n\nfig, ax = plt.subplots(figsize=(8, 6))\nim = ax.imshow(downstream.values, cmap='viridis', origin='lower')\nax.set_title('MFD downstream flow length')\nplt.colorbar(im, ax=ax, label='distance (coordinate units)')\nplt.tight_layout()\nplt.show()"
   },
   {
    "cell_type": "markdown",
@@ -127,16 +101,7 @@
    "id": "a7b8c9d0",
    "metadata": {},
    "outputs": [],
-   "source": [
-    "upstream = flow_length_mfd(mfd_dir, direction='upstream')\n",
-    "\n",
-    "fig, ax = plt.subplots(figsize=(8, 6))\n",
-    "im = ax.imshow(upstream.values, cmap='magma', origin='lower')\n",
-    "ax.set_title('MFD upstream flow length')\n",
-    "plt.colorbar(im, ax=ax, label='distance (coordinate units)')\n",
-    "plt.tight_layout()\n",
-    "plt.show()"
-   ]
+   "source": "upstream = flow_length(mfd_dir, routing='mfd', direction='upstream')\n\nfig, ax = plt.subplots(figsize=(8, 6))\nim = ax.imshow(upstream.values, cmap='magma', origin='lower')\nax.set_title('MFD upstream flow length')\nplt.colorbar(im, ax=ax, label='distance (coordinate units)')\nplt.tight_layout()\nplt.show()"
   },
   {
    "cell_type": "markdown",

From 1643e6a58198c37bda24f84ac36a4821c201eb70 Mon Sep 17 00:00:00 2001
From: Brendan Collins <brendancol@gmail.com>
Date: Thu, 25 Jun 2026 13:33:45 -0700
Subject: [PATCH 4/4] Resolve hydro *_d8 accessor doc sources from
 xrspatial.hydro (#3528)

CI caught two accessor sites that looked up the suffixed hydrology names on
the top-level xrspatial namespace by string, which this PR no longer exports:
- accessor._delegated_doc surfaces the documented *_d8 variant docstring on
  the .xrs hydro methods; fall back to xrspatial.hydro when the name is not
  top-level so help text is restored.
- test_accessor_docstring_matches_source resolves its fill_d8/watershed_d8
  source the same way.

Fixes the 4 test_accessor.py failures (docstring-match + every-method-
documented drift guard) on the full test lane.
---
 xrspatial/accessor.py            | 10 ++++++++--
 xrspatial/tests/test_accessor.py |  3 ++-
 2 files changed, 10 insertions(+), 3 deletions(-)

diff --git a/xrspatial/accessor.py b/xrspatial/accessor.py
index 8c1454c2a..cff20a752 100644
--- a/xrspatial/accessor.py
+++ b/xrspatial/accessor.py
@@ -2218,8 +2218,9 @@ def validate(self, *, raise_on_error=False):
 # ``help(da.xrs.slope)`` shows the same documentation as ``help(slope)``.
 # ---------------------------------------------------------------------------
 
-# Accessor method name -> name of the standalone function (in the top-level
-# ``xrspatial`` namespace) whose docstring should be surfaced. Only needed
+# Accessor method name -> name of the standalone function whose docstring
+# should be surfaced (resolved from the top-level ``xrspatial`` namespace, or
+# from ``xrspatial.hydro`` for the internal ``*_d8`` variants). Only needed
 # when the method name differs from the function name, or when the direct
 # delegate's docstring is a generic dispatcher: the hydrology unified wrappers
 # route by ``routing=`` and carry only a stub docstring, so their help text is
@@ -2252,6 +2253,11 @@ def _delegated_doc(method_name):
     import xrspatial
     source_name = _DOC_SOURCE_OVERRIDES.get(method_name, method_name)
     func = getattr(xrspatial, source_name, None)
+    if func is None:
+        # The hydrology ``*_d8`` doc sources are internal to xrspatial.hydro;
+        # only the routing wrappers are exported at the top level.
+        from . import hydro
+        func = getattr(hydro, source_name, None)
     return func.__doc__ if func is not None else None
 
 
diff --git a/xrspatial/tests/test_accessor.py b/xrspatial/tests/test_accessor.py
index 329ded8b4..2901cef88 100644
--- a/xrspatial/tests/test_accessor.py
+++ b/xrspatial/tests/test_accessor.py
@@ -473,7 +473,8 @@ def test_ds_plot_without_matplotlib_raises(monkeypatch, elevation):
     ('watershed', 'watershed_d8'),
 ])
 def test_accessor_docstring_matches_source(method_name, source):
-    func = getattr(xrspatial, source)
+    # hydro *_d8 doc sources are internal to xrspatial.hydro, not top level
+    func = getattr(xrspatial, source, None) or getattr(xrspatial.hydro, source)
     method = getattr(XrsSpatialDataArrayAccessor, method_name)
     assert inspect.getdoc(method), f'{method_name} has no docstring'
     assert inspect.getdoc(method) == inspect.getdoc(func)