Add D-inf and MFD variants of stream_order and stream_link

[brendancol]

Author of Proposal: (self)

Reason or problem

stream_order and stream_link only accept D8 flow direction grids. We already have D-infinity (Tarboton 1997) and Multiple Flow Direction (Freeman/Quinn) routing for flow direction and accumulation, but stream network analysis is D8-only. If you compute D-inf or MFD routing, you have to convert back to D8 before you can get stream orders or link IDs, which throws away the information those models provide.

Proposal

Add four new public functions:

• stream_order_dinf -- Strahler/Shreve ordering on D-inf angle grids
• stream_order_mfd -- Strahler/Shreve ordering on MFD fraction grids
• stream_link_dinf -- link segmentation on D-inf angle grids
• stream_link_mfd -- link segmentation on MFD fraction grids

Design:

Same Kahn's BFS topological sort as the D8 implementations. Only the topology graph construction differs:

• D-inf: The angle theta in [0, 2pi) picks two bracketing neighbors with proportional contributions. A cell flows to any neighbor receiving fraction > 0. In-degree counts how many upstream D-inf angles point at the current cell.

• MFD: The (8, H, W) fraction grid encodes flow proportions directly. A cell flows to every neighbor with fraction > 0. In-degree counts incoming non-zero fractions from upstream neighbors, same as flow_accumulation_mfd.

After the topology graph is built, Strahler/Shreve ordering and link segmentation work the same as D8.

Four new modules using the existing _dinf / _mfd naming convention:

• xrspatial/stream_order_dinf.py
• xrspatial/stream_order_mfd.py
• xrspatial/stream_link_dinf.py
• xrspatial/stream_link_mfd.py

All four backends: NumPy, CuPy, Dask+NumPy, Dask+CuPy.

Usage:

from xrspatial import flow_direction_dinf, flow_accumulation
from xrspatial import stream_order_dinf, stream_link_dinf

angles = flow_direction_dinf(elevation)
accum = flow_accumulation(elevation)
order = stream_order_dinf(angles, accum, threshold=500, method='strahler')
links = stream_link_dinf(angles, accum, threshold=500)

from xrspatial import flow_direction_mfd, flow_accumulation_mfd
from xrspatial import stream_order_mfd, stream_link_mfd

fractions = flow_direction_mfd(elevation)
accum = flow_accumulation_mfd(fractions)
order = stream_order_mfd(fractions, accum, threshold=500, method='strahler')
links = stream_link_mfd(fractions, accum, threshold=500)

Lets users run stream analysis end-to-end with D-inf or MFD routing without falling back to D8.

Stakeholders and impacts

Hydrology users who already use D-inf or MFD flow routing. No existing APIs change.

Drawbacks

Four new modules with shared logic means some code duplication with the D8 versions.

Alternatives

1. Add a routing parameter to the existing stream_order / stream_link functions. Those files are already 1500+ lines each and the input types differ (scalar codes vs. angles vs. 3D fraction arrays), so separate modules make more sense.

2. Do nothing; users convert D-inf/MFD to D8 first. Loses the distributional information.

Unresolved questions

Whether D-inf accumulation needs its own flow_accumulation_dinf function or if the existing D8 accumulation is fine. Any accumulation grid works as input since the threshold comparison is routing-agnostic.

Additional notes

The MFD topology engine can reuse the neighbor offset arrays (_DY, _DX, _OPPOSITE) from flow_accumulation_mfd.py. The D-inf angle-to-neighbor conversion follows the same facet logic as flow_direction_dinf.py.