Skip to content

Commit 001db21

Browse files
brendancolbrendancol
authored
Add D-inf and MFD stream ordering and link segmentation (#983) (#984)
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 All four support NumPy, CuPy, Dask+NumPy, and Dask+CuPy backends.
1 parent 11b3379 commit 001db21

File tree

12 files changed

+6628
-0
lines changed

12 files changed

+6628
-0
lines changed

‎README.md‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -296,7 +296,11 @@ In the GIS world, rasters are used for representing continuous phenomena (e.g. e
296296
| [Watershed](xrspatial/watershed.py) | Labels each cell with the pour point it drains to via D8 flow direction | Standard (D8 tracing) | ✅️ | ✅️ | ✅️ | ✅️ |
297297
| [Basins](xrspatial/watershed.py) | Delineates drainage basins by labeling each cell with its outlet ID | Standard (D8 tracing) | ✅️ | ✅️ | ✅️ | ✅️ |
298298
| [Stream Order](xrspatial/stream_order.py) | Assigns Strahler or Shreve stream order to cells in a drainage network | Strahler 1957, Shreve 1966 | ✅️ | ✅️ | ✅️ | ✅️ |
299+
| [Stream Order (Dinf)](xrspatial/stream_order_dinf.py) | Strahler/Shreve stream ordering on D-infinity flow direction grids | Tarboton 1997 | ✅️ | ✅️ | ✅️ | ✅️ |
300+
| [Stream Order (MFD)](xrspatial/stream_order_mfd.py) | Strahler/Shreve stream ordering on MFD fraction grids | Freeman 1991 | ✅️ | ✅️ | ✅️ | ✅️ |
299301
| [Stream Link](xrspatial/stream_link.py) | Assigns unique IDs to each stream segment between junctions | Standard | ✅️ | ✅️ | ✅️ | ✅️ |
302+
| [Stream Link (Dinf)](xrspatial/stream_link_dinf.py) | Stream link segmentation on D-infinity flow direction grids | Tarboton 1997 | ✅️ | ✅️ | ✅️ | ✅️ |
303+
| [Stream Link (MFD)](xrspatial/stream_link_mfd.py) | Stream link segmentation on MFD fraction grids | Freeman 1991 | ✅️ | ✅️ | ✅️ | ✅️ |
300304
| [Snap Pour Point](xrspatial/snap_pour_point.py) | Snaps pour points to the highest-accumulation cell within a search radius | Custom | ✅️ | ✅️ | ✅️ | ✅️ |
301305
| [Flow Path](xrspatial/flow_path.py) | Traces downstream flow paths from start points through a D8 direction grid | Standard (D8 tracing) | ✅️ | ✅️ | 🔄 | 🔄 |
302306
| [HAND](xrspatial/hand.py) | Computes Height Above Nearest Drainage by tracing D8 flow to the nearest stream cell | Nobre et al. 2011 | ✅️ | ✅️ | 🔄 | 🔄 |

‎docs/source/reference/hydrology.rst‎

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -96,13 +96,41 @@ Stream Link
9696

9797
xrspatial.stream_link.stream_link
9898

99+
Stream Link (D-infinity)
100+
========================
101+
.. autosummary::
102+
:toctree: _autosummary
103+
104+
xrspatial.stream_link_dinf.stream_link_dinf
105+
106+
Stream Link (MFD)
107+
=================
108+
.. autosummary::
109+
:toctree: _autosummary
110+
111+
xrspatial.stream_link_mfd.stream_link_mfd
112+
99113
Stream Order
100114
============
101115
.. autosummary::
102116
:toctree: _autosummary
103117

104118
xrspatial.stream_order.stream_order
105119

120+
Stream Order (D-infinity)
121+
=========================
122+
.. autosummary::
123+
:toctree: _autosummary
124+
125+
xrspatial.stream_order_dinf.stream_order_dinf
126+
127+
Stream Order (MFD)
128+
==================
129+
.. autosummary::
130+
:toctree: _autosummary
131+
132+
xrspatial.stream_order_mfd.stream_order_mfd
133+
106134
Height Above Nearest Drainage (HAND)
107135
=====================================
108136
.. autosummary::

‎examples/user_guide/25_Stream_Analysis_Dinf_MFD.ipynb‎

Lines changed: 249 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,249 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"# Stream ordering and link segmentation with D-inf and MFD routing\n",
8+
"\n",
9+
"xarray-spatial provides stream network analysis (ordering and link segmentation) for three flow routing models:\n",
10+
"\n",
11+
"- **D8**: single steepest-descent neighbor (existing `stream_order`, `stream_link`)\n",
12+
"- **D-infinity**: continuous angle distributing flow between two neighbors (Tarboton 1997)\n",
13+
"- **MFD**: flow fractions distributed to all downslope neighbors (Freeman/Quinn)\n",
14+
"\n",
15+
"This notebook shows how to use the D-inf and MFD variants."
16+
]
17+
},
18+
{
19+
"cell_type": "code",
20+
"execution_count": null,
21+
"metadata": {},
22+
"outputs": [],
23+
"source": [
24+
"import numpy as np\n",
25+
"import xarray as xr\n",
26+
"import matplotlib.pyplot as plt\n",
27+
"\n",
28+
"from xrspatial import generate_terrain\n",
29+
"from xrspatial import flow_direction, flow_direction_dinf, flow_direction_mfd\n",
30+
"from xrspatial import flow_accumulation, flow_accumulation_mfd\n",
31+
"from xrspatial import stream_order, stream_link\n",
32+
"from xrspatial import stream_order_dinf, stream_link_dinf\n",
33+
"from xrspatial import stream_order_mfd, stream_link_mfd"
34+
]
35+
},
36+
{
37+
"cell_type": "markdown",
38+
"metadata": {},
39+
"source": [
40+
"## Generate synthetic terrain"
41+
]
42+
},
43+
{
44+
"cell_type": "code",
45+
"execution_count": null,
46+
"metadata": {},
47+
"outputs": [],
48+
"source": [
49+
"W, H = 400, 400\n",
50+
"cvs_x = np.linspace(0, 1000, W)\n",
51+
"cvs_y = np.linspace(0, 1000, H)\n",
52+
"\n",
53+
"terrain = generate_terrain(x_range=(0, 1000), y_range=(0, 1000),\n",
54+
" width=W, height=H)\n",
55+
"\n",
56+
"fig, ax = plt.subplots(figsize=(6, 5))\n",
57+
"terrain.plot(ax=ax, cmap='terrain')\n",
58+
"ax.set_title('Synthetic elevation')\n",
59+
"ax.set_aspect('equal')\n",
60+
"plt.tight_layout()\n",
61+
"plt.show()"
62+
]
63+
},
64+
{
65+
"cell_type": "markdown",
66+
"metadata": {},
67+
"source": [
68+
"## Compute flow directions and accumulation\n",
69+
"\n",
70+
"We compute all three routing models from the same elevation surface."
71+
]
72+
},
73+
{
74+
"cell_type": "code",
75+
"execution_count": null,
76+
"metadata": {},
77+
"outputs": [],
78+
"source": [
79+
"# D8\n",
80+
"fd_d8 = flow_direction(terrain)\n",
81+
"fa_d8 = flow_accumulation(fd_d8)\n",
82+
"\n",
83+
"# D-infinity\n",
84+
"fd_dinf = flow_direction_dinf(terrain)\n",
85+
"\n",
86+
"# MFD\n",
87+
"fd_mfd = flow_direction_mfd(terrain)\n",
88+
"fa_mfd = flow_accumulation_mfd(fd_mfd)\n",
89+
"\n",
90+
"print(f'D8 flow dir shape: {fd_d8.shape}')\n",
91+
"print(f'D-inf angles shape: {fd_dinf.shape}')\n",
92+
"print(f'MFD fractions shape: {fd_mfd.shape}')"
93+
]
94+
},
95+
{
96+
"cell_type": "markdown",
97+
"metadata": {},
98+
"source": [
99+
"## Stream ordering: D8 vs D-inf vs MFD\n",
100+
"\n",
101+
"We extract stream networks using a threshold and compare Strahler ordering across routing models."
102+
]
103+
},
104+
{
105+
"cell_type": "code",
106+
"execution_count": null,
107+
"metadata": {},
108+
"outputs": [],
109+
"source": [
110+
"threshold = 200\n",
111+
"\n",
112+
"# D8 stream order (uses D8 accumulation)\n",
113+
"so_d8 = stream_order(fd_d8, fa_d8, threshold=threshold, method='strahler')\n",
114+
"\n",
115+
"# D-inf stream order (uses D8 accumulation for thresholding)\n",
116+
"so_dinf = stream_order_dinf(fd_dinf, fa_d8, threshold=threshold, method='strahler')\n",
117+
"\n",
118+
"# MFD stream order (uses MFD accumulation for thresholding)\n",
119+
"so_mfd = stream_order_mfd(fd_mfd, fa_mfd, threshold=threshold, method='strahler')"
120+
]
121+
},
122+
{
123+
"cell_type": "code",
124+
"execution_count": null,
125+
"metadata": {},
126+
"outputs": [],
127+
"source": [
128+
"fig, axes = plt.subplots(1, 3, figsize=(15, 5))\n",
129+
"\n",
130+
"for ax, data, title in zip(axes,\n",
131+
" [so_d8, so_dinf, so_mfd],\n",
132+
" ['D8', 'D-infinity', 'MFD']):\n",
133+
" im = ax.imshow(np.where(np.isnan(data.values), 0, data.values),\n",
134+
" cmap='Blues', interpolation='nearest',\n",
135+
" vmin=0, vmax=max(np.nanmax(so_d8.values),\n",
136+
" np.nanmax(so_dinf.values),\n",
137+
" np.nanmax(so_mfd.values)))\n",
138+
" ax.set_title(f'Strahler order ({title})')\n",
139+
" ax.set_aspect('equal')\n",
140+
"\n",
141+
"fig.colorbar(im, ax=axes, shrink=0.6, label='Stream order')\n",
142+
"plt.tight_layout()\n",
143+
"plt.show()"
144+
]
145+
},
146+
{
147+
"cell_type": "markdown",
148+
"metadata": {},
149+
"source": [
150+
"## Shreve magnitude comparison"
151+
]
152+
},
153+
{
154+
"cell_type": "code",
155+
"execution_count": null,
156+
"metadata": {},
157+
"outputs": [],
158+
"source": [
159+
"sv_d8 = stream_order(fd_d8, fa_d8, threshold=threshold, method='shreve')\n",
160+
"sv_dinf = stream_order_dinf(fd_dinf, fa_d8, threshold=threshold, method='shreve')\n",
161+
"sv_mfd = stream_order_mfd(fd_mfd, fa_mfd, threshold=threshold, method='shreve')\n",
162+
"\n",
163+
"fig, axes = plt.subplots(1, 3, figsize=(15, 5))\n",
164+
"\n",
165+
"for ax, data, title in zip(axes,\n",
166+
" [sv_d8, sv_dinf, sv_mfd],\n",
167+
" ['D8', 'D-infinity', 'MFD']):\n",
168+
" vals = np.where(np.isnan(data.values), 0, data.values)\n",
169+
" im = ax.imshow(np.log1p(vals), cmap='viridis', interpolation='nearest')\n",
170+
" ax.set_title(f'Shreve magnitude ({title})')\n",
171+
" ax.set_aspect('equal')\n",
172+
"\n",
173+
"fig.colorbar(im, ax=axes, shrink=0.6, label='log(1 + magnitude)')\n",
174+
"plt.tight_layout()\n",
175+
"plt.show()"
176+
]
177+
},
178+
{
179+
"cell_type": "markdown",
180+
"metadata": {},
181+
"source": [
182+
"## Stream link segmentation"
183+
]
184+
},
185+
{
186+
"cell_type": "code",
187+
"execution_count": null,
188+
"metadata": {},
189+
"outputs": [],
190+
"source": [
191+
"sl_d8 = stream_link(fd_d8, fa_d8, threshold=threshold)\n",
192+
"sl_dinf = stream_link_dinf(fd_dinf, fa_d8, threshold=threshold)\n",
193+
"sl_mfd = stream_link_mfd(fd_mfd, fa_mfd, threshold=threshold)\n",
194+
"\n",
195+
"fig, axes = plt.subplots(1, 3, figsize=(15, 5))\n",
196+
"\n",
197+
"for ax, data, title in zip(axes,\n",
198+
" [sl_d8, sl_dinf, sl_mfd],\n",
199+
" ['D8', 'D-infinity', 'MFD']):\n",
200+
" vals = data.values.copy()\n",
201+
" # Color by link ID mod a palette size for visibility\n",
202+
" vals_display = np.where(np.isnan(vals), 0,\n",
203+
" np.mod(vals, 20) + 1)\n",
204+
" ax.imshow(vals_display, cmap='tab20', interpolation='nearest')\n",
205+
" ax.set_title(f'Stream links ({title})')\n",
206+
" ax.set_aspect('equal')\n",
207+
"\n",
208+
"plt.tight_layout()\n",
209+
"plt.show()"
210+
]
211+
},
212+
{
213+
"cell_type": "markdown",
214+
"metadata": {},
215+
"source": [
216+
"## Counting stream segments\n",
217+
"\n",
218+
"Each routing model produces a slightly different stream network topology.\n",
219+
"MFD and D-inf can produce more junction points than D8 because flow\n",
220+
"splits across multiple neighbors."
221+
]
222+
},
223+
{
224+
"cell_type": "code",
225+
"execution_count": null,
226+
"metadata": {},
227+
"outputs": [],
228+
"source": [
229+
"for label, data in [('D8', sl_d8), ('D-inf', sl_dinf), ('MFD', sl_mfd)]:\n",
230+
" n_links = len(np.unique(data.values[~np.isnan(data.values)]))\n",
231+
" n_stream = int(np.sum(~np.isnan(data.values)))\n",
232+
" print(f'{label:6s}: {n_links:4d} links, {n_stream:6d} stream cells')"
233+
]
234+
}
235+
],
236+
"metadata": {
237+
"kernelspec": {
238+
"display_name": "Python 3",
239+
"language": "python",
240+
"name": "python3"
241+
},
242+
"language_info": {
243+
"name": "python",
244+
"version": "3.10.0"
245+
}
246+
},
247+
"nbformat": 4,
248+
"nbformat_minor": 4
249+
}

‎xrspatial/__init__.py‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -74,7 +74,11 @@
7474
from xrspatial.sink import sink # noqa
7575
from xrspatial.snap_pour_point import snap_pour_point # noqa
7676
from xrspatial.stream_link import stream_link # noqa
77+
from xrspatial.stream_link_dinf import stream_link_dinf # noqa
78+
from xrspatial.stream_link_mfd import stream_link_mfd # noqa
7779
from xrspatial.stream_order import stream_order # noqa
80+
from xrspatial.stream_order_dinf import stream_order_dinf # noqa
81+
from xrspatial.stream_order_mfd import stream_order_mfd # noqa
7882
from xrspatial.sky_view_factor import sky_view_factor # noqa
7983
from xrspatial.slope import slope # noqa
8084
from xrspatial.surface_distance import surface_allocation # noqa

0 commit comments

Comments
 (0)