Hydrology pipeline¶

Real-world DEM-to-streams workflow covering:

DEM
  → full_hydro_pipeline   (fill → flow_direction → accumulate → streams)
  → order(method="hack")
  → order(method="topological")
  → to_vector(...) carrying sinuosity
  → main_stem(...)
  → prune_short(...)
  → hand(method="euclidean")

We build a synthetic mountain catchment with a confluence and exercise every stage on it.

In [1]:

import matplotlib.pyplot as plt
plt.rcParams["figure.dpi"] = 80
plt.rcParams["savefig.dpi"] = 80

import matplotlib.pyplot as plt plt.rcParams["figure.dpi"] = 80 plt.rcParams["savefig.dpi"] = 80

In [2]:

import numpy as np
from pyramids.dataset import Dataset
from digitalrivers import DEM

rows, cols = 20, 20
z = np.full((rows, cols), 100.0, dtype=np.float32)
for r in range(rows):
    z[r, 10] = float(40 - 2 * r)  # main stem descends south down col 10
for c in range(10):
    z[8, c] = float(35 - c)        # west tributary into (8, 10)
for c in range(11, cols):
    z[13, c] = float(25 - (c - 11))  # east tributary into (13, 10)

ds = Dataset.create_from_array(
    z, top_left_corner=(0.0, 0.0), cell_size=30.0, epsg=32618, no_data_value=-9999.0,
)
dem = DEM(ds.raster)
print(f"DEM shape={z.shape}  elevation range=[{z.min():.1f}, {z.max():.1f}] m")

import numpy as np from pyramids.dataset import Dataset from digitalrivers import DEM rows, cols = 20, 20 z = np.full((rows, cols), 100.0, dtype=np.float32) for r in range(rows): z[r, 10] = float(40 - 2 * r) # main stem descends south down col 10 for c in range(10): z[8, c] = float(35 - c) # west tributary into (8, 10) for c in range(11, cols): z[13, c] = float(25 - (c - 11)) # east tributary into (13, 10) ds = Dataset.create_from_array( z, top_left_corner=(0.0, 0.0), cell_size=30.0, epsg=32618, no_data_value=-9999.0, ) dem = DEM(ds.raster) print(f"DEM shape={z.shape} elevation range=[{z.min():.1f}, {z.max():.1f}] m")

2026-06-09 23:45:36 | INFO | pyramids.base.config | Logging is configured.

DEM shape=(20, 20)  elevation range=[2.0, 100.0] m

Visualise the synthetic catchment¶

In [3]:

import matplotlib.pyplot as plt
from matplotlib.colors import LightSource

elev = dem.values
ls = LightSource(azdeg=315, altdeg=45)
shade = ls.hillshade(elev, vert_exag=2.0, dx=30.0, dy=30.0)

fig, axes = plt.subplots(1, 2, figsize=(10, 4))
dem.plot(band=0, ax=axes[0], title="Synthetic DEM (m)", cmap="terrain")
axes[0].set_xticks([])
axes[0].set_yticks([])

axes[1].imshow(shade, cmap="gray", origin="upper")
axes[1].imshow(elev, cmap="terrain", origin="upper", alpha=0.55)
axes[1].set_title("Hillshade + elevation"); axes[1].axis("off")
fig.tight_layout()
plt.show()

import matplotlib.pyplot as plt from matplotlib.colors import LightSource elev = dem.values ls = LightSource(azdeg=315, altdeg=45) shade = ls.hillshade(elev, vert_exag=2.0, dx=30.0, dy=30.0) fig, axes = plt.subplots(1, 2, figsize=(10, 4)) dem.plot(band=0, ax=axes[0], title="Synthetic DEM (m)", cmap="terrain") axes[0].set_xticks([]) axes[0].set_yticks([]) axes[1].imshow(shade, cmap="gray", origin="upper") axes[1].imshow(elev, cmap="terrain", origin="upper", alpha=0.55) axes[1].set_title("Hillshade + elevation"); axes[1].axis("off") fig.tight_layout() plt.show()

C:\Users\main\AppData\Local\Temp\ipykernel_38692\2999521973.py:16: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  fig.tight_layout()

composite hydro pipeline¶

DEM.full_hydro_pipeline(...) chains fill → flow_direction → accumulate → streams in one call and returns a dict of typed results.

In [4]:

bundle = dem.full_hydro_pipeline(
    fill_method="priority_flood",
    flow_method="d8",
    stream_threshold_cells=4,
)
filled = bundle["filled_dem"]
fdir = bundle["flow_direction"]
acc = bundle["accumulation"]
streams = bundle["streams"]

n_stream = int(streams.read_array().astype(bool).sum())
print(f"filled_dem  type={type(filled).__name__}")
print(f"fdir        routing={fdir.routing}")
print(f"acc         max={int(acc.read_array().max())}  routing={acc.routing}")
print(f"streams     cells={n_stream}  threshold={streams.threshold}")

bundle = dem.full_hydro_pipeline( fill_method="priority_flood", flow_method="d8", stream_threshold_cells=4, ) filled = bundle["filled_dem"] fdir = bundle["flow_direction"] acc = bundle["accumulation"] streams = bundle["streams"] n_stream = int(streams.read_array().astype(bool).sum()) print(f"filled_dem type={type(filled).__name__}") print(f"fdir routing={fdir.routing}") print(f"acc max={int(acc.read_array().max())} routing={acc.routing}") print(f"streams cells={n_stream} threshold={streams.threshold}")

filled_dem  type=DEM
fdir        routing=d8
acc         max=86  routing=d8
streams     cells=35  threshold=4.0

Visualise the hydro-pipeline outputs¶

In [5]:

import matplotlib.pyplot as plt
from pyramids.dataset import Dataset

sm = streams.read_array().astype(bool)
fd_arr = fdir.read_array()
acc_arr = acc.read_array()

fig, axes = plt.subplots(2, 2, figsize=(10, 9))
filled.plot(band=0, ax=axes[0, 0], title="Filled DEM", cmap="terrain")
axes[0, 0].set_xticks([])
axes[0, 0].set_yticks([])

fd_show = np.where(fd_arr < 0, np.nan, fd_arr).astype(np.float32)
Dataset.dataset_like(filled, fd_show).plot(
    band=0, ax=axes[0, 1], title="D8 flow direction", cmap="tab10")
axes[0, 1].set_xticks([])
axes[0, 1].set_yticks([])

acc_show = np.where(acc_arr < 1, np.nan, acc_arr).astype(np.float32)
Dataset.dataset_like(filled, acc_show).plot(
    band=0, ax=axes[1, 0], title="Flow accumulation (log)", cmap="viridis",
    color_scale="power")
axes[1, 0].set_xticks([])
axes[1, 0].set_yticks([])

stream_show = np.where(sm, 1.0, np.nan).astype(np.float32)
Dataset.dataset_like(filled, stream_show).plot(
    band=0, ax=axes[1, 1], title="Streams", cmap="Blues")
axes[1, 1].set_xticks([])
axes[1, 1].set_yticks([])
fig.tight_layout()
plt.show()

import matplotlib.pyplot as plt from pyramids.dataset import Dataset sm = streams.read_array().astype(bool) fd_arr = fdir.read_array() acc_arr = acc.read_array() fig, axes = plt.subplots(2, 2, figsize=(10, 9)) filled.plot(band=0, ax=axes[0, 0], title="Filled DEM", cmap="terrain") axes[0, 0].set_xticks([]) axes[0, 0].set_yticks([]) fd_show = np.where(fd_arr < 0, np.nan, fd_arr).astype(np.float32) Dataset.dataset_like(filled, fd_show).plot( band=0, ax=axes[0, 1], title="D8 flow direction", cmap="tab10") axes[0, 1].set_xticks([]) axes[0, 1].set_yticks([]) acc_show = np.where(acc_arr < 1, np.nan, acc_arr).astype(np.float32) Dataset.dataset_like(filled, acc_show).plot( band=0, ax=axes[1, 0], title="Flow accumulation (log)", cmap="viridis", color_scale="power") axes[1, 0].set_xticks([]) axes[1, 0].set_yticks([]) stream_show = np.where(sm, 1.0, np.nan).astype(np.float32) Dataset.dataset_like(filled, stream_show).plot( band=0, ax=axes[1, 1], title="Streams", cmap="Blues") axes[1, 1].set_xticks([]) axes[1, 1].set_yticks([]) fig.tight_layout() plt.show()

C:\Users\main\AppData\Local\Temp\ipykernel_38692\2282631296.py:31: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  fig.tight_layout()

Hack stream ordering¶

Order 1 is assigned to the main stem (longest source-to-outlet path). Every tributary joining the main stem becomes order 2; tributaries of those order 3; and so on.

In [6]:

hack = streams.order(method="hack", flow_direction=fdir)
hack_arr = hack.read_array()
sm = streams.read_array().astype(bool)
print(f"Hack orders on stream cells: {sorted(set(int(v) for v in hack_arr[sm]))}")
print(f"Order-1 (main stem) cells: {int((hack_arr == 1).sum())}")
print(f"Order-2 (1st tributary) cells: {int((hack_arr == 2).sum())}")

hack = streams.order(method="hack", flow_direction=fdir) hack_arr = hack.read_array() sm = streams.read_array().astype(bool) print(f"Hack orders on stream cells: {sorted(set(int(v) for v in hack_arr[sm]))}") print(f"Order-1 (main stem) cells: {int((hack_arr == 1).sum())}") print(f"Order-2 (1st tributary) cells: {int((hack_arr == 2).sum())}")

Hack orders on stream cells: [1, 2]
Order-1 (main stem) cells: 27
Order-2 (1st tributary) cells: 8

Topological stream ordering¶

Kahn-sort sequential numbering. Heads are enumerated in row-major order first, then progressively downstream. The outlet's index equals the total stream-cell count.

In [7]:

topo = streams.order(method="topological", flow_direction=fdir)
topo_arr = topo.read_array()
print(f"Topological max (outlet index): {int(topo_arr.max())}")
print(f"Stream cell count:              {n_stream}")
assert int(topo_arr.max()) == n_stream

topo = streams.order(method="topological", flow_direction=fdir) topo_arr = topo.read_array() print(f"Topological max (outlet index): {int(topo_arr.max())}") print(f"Stream cell count: {n_stream}") assert int(topo_arr.max()) == n_stream

Topological max (outlet index): 35
Stream cell count:              35

Visualise: Hack vs Topological ordering¶

In [8]:

import matplotlib.pyplot as plt
from pyramids.dataset import Dataset

fig, axes = plt.subplots(1, 2, figsize=(10, 4))
for ax, arr, title in zip(
    axes, (hack_arr, topo_arr), ("Hack order", "Topological index"),
):
    show = np.where(sm, arr.astype(np.float32), np.nan).astype(np.float32)
    Dataset.dataset_like(filled, show).plot(
        band=0, ax=ax, title=title, cmap="plasma")
    ax.set_xticks([])
    ax.set_yticks([])
fig.tight_layout()
plt.show()

import matplotlib.pyplot as plt from pyramids.dataset import Dataset fig, axes = plt.subplots(1, 2, figsize=(10, 4)) for ax, arr, title in zip( axes, (hack_arr, topo_arr), ("Hack order", "Topological index"), ): show = np.where(sm, arr.astype(np.float32), np.nan).astype(np.float32) Dataset.dataset_like(filled, show).plot( band=0, ax=ax, title=title, cmap="plasma") ax.set_xticks([]) ax.set_yticks([]) fig.tight_layout() plt.show()

C:\Users\main\AppData\Local\Temp\ipykernel_38692\1151883421.py:13: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  fig.tight_layout()

Per-link sinuosity¶

StreamRaster.to_vector(...) now carries a sinuosity column (traced length / straight-line distance per link). Straight links have sinuosity exactly 1.0; meandering links exceed 1.0.

In [9]:

links = streams.to_vector(fdir, dem=filled)
print(links[["link_id", "length_m", "drop_m", "mean_slope", "sinuosity"]])

links = streams.to_vector(fdir, dem=filled) print(links[["link_id", "length_m", "drop_m", "mean_slope", "sinuosity"]])

   link_id    length_m  drop_m  mean_slope  sinuosity
0        0  240.000000    16.0    0.066667   1.000000
1        1  282.426407    12.0    0.042489   1.039626
2        2  300.000000    20.0    0.066667   1.000000
3        3  180.000000     6.0    0.033333   1.000000

Main stem¶

Public API for the longest source-to-outlet path. Returns a binary mask.

In [10]:

main = streams.main_stem(fdir)
print(f"Main-stem cells: {int(main.sum())} of {n_stream} stream cells")
# Verify Hack-order 1 == main_stem mask
assert (hack_arr[main] == 1).all()

main = streams.main_stem(fdir) print(f"Main-stem cells: {int(main.sum())} of {n_stream} stream cells") # Verify Hack-order 1 == main_stem mask assert (hack_arr[main] == 1).all()

Main-stem cells: 20 of 35 stream cells

Remove short stream links¶

prune_short drops headwater links below min_length_m. Internal links (between two confluences) stay so the network topology is preserved.

In [11]:

pruned = streams.prune_short(fdir, min_length_m=120.0)
n_after = int(pruned.read_array().astype(bool).sum())
print(f"Stream cells before prune: {n_stream}")
print(f"Stream cells after  prune: {n_after}")

pruned = streams.prune_short(fdir, min_length_m=120.0) n_after = int(pruned.read_array().astype(bool).sum()) print(f"Stream cells before prune: {n_stream}") print(f"Stream cells after prune: {n_after}")

Stream cells before prune: 35
Stream cells after  prune: 35

Visualise: main stem + pre/post-prune comparison¶

In [12]:

import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 2, figsize=(10, 4))
axes[0].imshow(filled.values, cmap="gray", origin="upper", alpha=0.5)
axes[0].imshow(np.where(sm, 1, np.nan), cmap="Blues", origin="upper", alpha=0.7)
axes[0].imshow(np.where(main, 1, np.nan), cmap="autumn", origin="upper")
axes[0].set_title("Main stem (red) on stream raster (blue)"); axes[0].axis("off")

after = pruned.read_array().astype(bool)
axes[1].imshow(filled.values, cmap="gray", origin="upper", alpha=0.5)
axes[1].imshow(np.where(sm & ~after, 1, np.nan), cmap="autumn", origin="upper", alpha=0.9)
axes[1].imshow(np.where(after, 1, np.nan), cmap="Blues", origin="upper", alpha=0.7)
axes[1].set_title("After prune_short — red cells dropped"); axes[1].axis("off")
fig.tight_layout()
plt.show()

import matplotlib.pyplot as plt fig, axes = plt.subplots(1, 2, figsize=(10, 4)) axes[0].imshow(filled.values, cmap="gray", origin="upper", alpha=0.5) axes[0].imshow(np.where(sm, 1, np.nan), cmap="Blues", origin="upper", alpha=0.7) axes[0].imshow(np.where(main, 1, np.nan), cmap="autumn", origin="upper") axes[0].set_title("Main stem (red) on stream raster (blue)"); axes[0].axis("off") after = pruned.read_array().astype(bool) axes[1].imshow(filled.values, cmap="gray", origin="upper", alpha=0.5) axes[1].imshow(np.where(sm & ~after, 1, np.nan), cmap="autumn", origin="upper", alpha=0.9) axes[1].imshow(np.where(after, 1, np.nan), cmap="Blues", origin="upper", alpha=0.7) axes[1].set_title("After prune_short — red cells dropped"); axes[1].axis("off") fig.tight_layout() plt.show()

Euclidean HAND¶

Height Above Nearest Drainage using 2D Euclidean distance to the nearest stream cell. Stream cells are 0; cells far from streams have larger HAND.

In [13]:

hand_eu = filled.hand(streams, method="euclidean")
hand_arr = hand_eu.read_array()
print(f"HAND range on valid cells: [{hand_arr[hand_arr != -9999.0].min():.2f}, "
      f"{hand_arr[hand_arr != -9999.0].max():.2f}] m")
# Stream cells must report HAND = 0
assert (hand_arr[sm] == 0).all()

hand_eu = filled.hand(streams, method="euclidean") hand_arr = hand_eu.read_array() print(f"HAND range on valid cells: [{hand_arr[hand_arr != -9999.0].min():.2f}, " f"{hand_arr[hand_arr != -9999.0].max():.2f}] m") # Stream cells must report HAND = 0 assert (hand_arr[sm] == 0).all()

HAND range on valid cells: [0.00, 98.00] m

Visualise: Euclidean HAND¶

In [14]:

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(6, 5))
hand_eu.plot(band=0, ax=ax, title="Euclidean HAND (m)", cmap="terrain_r",
             exclude_value=-9999.0)
ax.set_xticks([])
ax.set_yticks([])
fig.tight_layout()
plt.show()

import matplotlib.pyplot as plt fig, ax = plt.subplots(figsize=(6, 5)) hand_eu.plot(band=0, ax=ax, title="Euclidean HAND (m)", cmap="terrain_r", exclude_value=-9999.0) ax.set_xticks([]) ax.set_yticks([]) fig.tight_layout() plt.show()

Summary¶

From a single 20×20 synthetic DEM we built:

• A complete hydro pipeline via one full_hydro_pipeline call.
• Two stream-ordering rasters — Hack and Topological.
• A vector network with sinuosity per link.
• A main-stem binary mask.
• A pruned stream network with short headwater links removed.
• A Euclidean HAND raster.