How to Calculate Strahler Stream Order in Python
Problem statement
Strahler order is one of the first attributes asked of a river network: for sampling designs, habitat classes, map symbology and regional statistics. In Python there are three common starting points — a DEM routed with pysheds, a DEM routed with WhiteboxTools, or an existing vector network such as NHDPlus — and each needs a slightly different route to an order on every stream segment. Each also has a trap that produces plausible but wrong orders.
Computed for the Esopus Creek basin above Coldbrook, New York, on the breached 10 m 3DEP DEM with channels defined at 0.5 km², and on the USGS NHDPlus High Resolution network:
- pysheds ordered and vectorised the network in 2.9 s: 249 first-order segments up to 18 fifth-order ones. Sampling the order raster with
int()on the segment coordinates first gave order 0 to 395 of 501 segments, because the coordinates lie exactly on cell corners. - WhiteboxTools went from D8 pointer to ordered lines in 1.3 s, writing a shapefile with the order in a
STRM_VALfield and no coordinate reference system. - Ordering NHDPlus HR from its
fromnodeandtonodeattributes matched NHDPlus' ownstreamordeon 99.4% of flowlines and 99.8% of length — once the 23 minor divergence paths were removed. - With the divergences left in, agreement fell to 91.5% and the highest computed order reached 12, against 6 in NHDPlus: every braid counted as a junction of two equal streams.
Quick answer
From a pysheds flow direction grid and a channel mask:
order = grid.stream_order(fdir, acc > threshold_cells) # raster of Strahler ordersFrom a vector network with from/to nodes:
import networkx as nx
graph = nx.MultiDiGraph()
for line_id, a, b in zip(lines.id, lines.from_node, lines.to_node):
graph.add_edge(a, b, key=line_id)
order = {}
for node in nx.topological_sort(graph):
up = sorted((order[k] for _, _, k in graph.in_edges(node, keys=True)), reverse=True)
level = 1 if not up else up[0] + int(len(up) > 1 and up[0] == up[1])
for _, _, k in graph.out_edges(node, keys=True):
order[k] = levelRemove divergent (braided) branches before ordering a vector network, and say which threshold and which rule produced the orders.
Step-by-step solution
1. Decide where the network comes from
A DEM-derived network gives orders consistent with your flow routing and threshold. A mapped network gives orders consistent with the map. They answer different questions; see stream order explained for why the same river can be order 5 in one and 7 in the other.
2. For a DEM, route and choose a threshold
Condition the DEM, compute D8 flow directions and accumulation, and define channels by an accumulation threshold. At 0.5 km² — 5,000 cells of 10 m — the Esopus network had 501 segments. The threshold sets the order: see extracting a stream network.
3. Order the raster with pysheds
grid.stream_order(fdir, channels) returns an integer raster with 0 off the network. On the clipped catchment, cell counts by order were 19,678, 9,456, 4,417, 3,364 and 909 for orders 1 to 5.
4. Attach the order to each vector segment
grid.extract_river_network returns one line per segment, but no order. Sample the order raster at a vertex in the middle of each line. The vertices are cell corners, not centres: the first vertex was at column 906.000000, row 69.000000. int() of a coordinate that comes back as 905.9999999 picks the neighbouring cell, off the channel; floor with a small tolerance instead (Example 1).
5. Or let WhiteboxTools write ordered lines
strahler_stream_order orders a stream raster using a D8 pointer, and raster_streams_to_vector turns it into lines carrying the order in STRM_VAL. The shapefile had no .prj: set its CRS from the DEM before clipping or measuring (Example 2).
6. For a vector network, build a directed graph
Mapped networks carry topology: NHDPlus HR has fromnode and tonode on every flowline. Build a directed multigraph with each flowline as an edge, visit nodes in topological order, and give every outgoing edge the order implied by the incoming ones. Without node attributes, derive nodes from snapped line endpoints, with lines digitised downstream.
7. Remove divergences before ordering
Where a river splits around an island and rejoins, both branches carry the same order, and the rejoining junction looks like two equal streams meeting — so the order rises without any tributary. The Esopus network had 23 minor divergence paths (divergence == 2). Ordered with them, orders climbed to 12; without them, the highest order was 5, and 1,914 of 1,925 flowlines matched NHDPlus.
8. Check against a reference
Compare with a published order where one exists, and look at the disagreements. Of the 11 remaining mismatches with NHDPlus, 7 were order 2 in NHDPlus but order 1 here, and 4 were first-order here but order 4, 5 or 6 in NHDPlus.
Code examples
Example 1 — pysheds: order raster to ordered segments
import time
import numpy as np
if not hasattr(np, "in1d"): # pysheds 0.5 on NumPy 2.4+
np.in1d = lambda a, b, **kw: np.isin(np.ravel(a), b, **kw)
import geopandas as gpd
import pandas as pd
from pysheds.grid import Grid
from shapely.geometry import shape
grid = Grid.from_raster("esopus_breached.tif")
fdir = grid.flowdir(grid.read_raster("esopus_breached.tif"))
acc = grid.accumulation(fdir)
x, y = grid.snap_to_mask(acc > 10_000, (560427.2, 4651640.7)) # the Coldbrook gauge, on a channel > 1 km2
grid.clip_to(grid.catchment(x=x, y=y, fdir=fdir, xytype="coordinate"))
fdir_c, acc_c = grid.view(fdir), grid.view(acc)
channels = acc_c > 5_000 # 0.5 km2 of 10 m cells
start = time.perf_counter()
order = np.asarray(grid.stream_order(fdir_c, channels))
network = grid.extract_river_network(fdir_c, channels)
segments = gpd.GeoDataFrame(geometry=[shape(f["geometry"]) for f in network["features"]], crs=grid.crs.srs)
print(f"ordered and vectorised in {time.perf_counter() - start:.1f} s")
first_col, first_row = ~grid.affine * segments.geometry.iloc[0].coords[0]
print(f"first vertex in cell units: col {first_col:.6f}, row {first_row:.6f}")
def order_at_middle(line):
col, row = ~grid.affine * line.coords[len(line.coords) // 2]
return int(order[int(np.floor(row + 1e-6)), int(np.floor(col + 1e-6))])
segments["strahler"] = segments.geometry.map(order_at_middle)
by_order = segments.assign(km=segments.length / 1000).groupby("strahler").agg(segments=("km", "size"), km=("km", "sum"))
by_order["cells"] = pd.Series(order[order > 0]).value_counts().sort_index()
print(by_order.round(1).to_string())ordered and vectorised in 2.9 s
first vertex in cell units: col 906.000000, row 69.000000
segments km cells
strahler
1 249 230.3 19678
2 118 110.8 9456
3 65 52.0 4417
4 51 39.5 3364
5 18 10.9 909
Example 2 — WhiteboxTools: ordered lines in one pass
import os
import whitebox
os.makedirs("wbt", exist_ok=True)
wbt = whitebox.WhiteboxTools()
wbt.set_verbose_mode(False)
wbt.set_working_dir(os.path.abspath("wbt"))
start = time.perf_counter()
wbt.d8_pointer(os.path.abspath("esopus_breached.tif"), "d8.tif")
wbt.d8_flow_accumulation("d8.tif", "acc.tif", out_type="cells", pntr=True)
wbt.extract_streams("acc.tif", "streams.tif", threshold=5000, zero_background=True)
wbt.strahler_stream_order("d8.tif", "streams.tif", "strahler.tif", zero_background=True)
wbt.raster_streams_to_vector("strahler.tif", "d8.tif", "strahler_lines.shp")
print(f"WhiteboxTools pointer to ordered lines in {time.perf_counter() - start:.1f} s")
lines = gpd.read_file("wbt/strahler_lines.shp")
print("fields:", [c for c in lines.columns if c != "geometry"], "crs:", lines.crs)
lines = lines.set_crs(32618, allow_override=True)
basin = gpd.read_file("nldi_basin_01362500.geojson").to_crs(32618)
lines = gpd.clip(lines, basin)
value_field = [c for c in lines.columns if c not in ("geometry", "FID")][0]
print(lines.assign(km=lines.length / 1000).groupby(value_field).agg(lines=("km", "size"), km=("km", "sum")).round(1).to_string())WhiteboxTools pointer to ordered lines in 1.3 s
fields: ['FID', 'STRM_VAL'] crs: None
lines km
STRM_VAL
1.0 254 231.9
2.0 118 110.8
3.0 64 52.0
4.0 50 39.5
5.0 21 12.6
The counts differ slightly from pysheds because this network was clipped to the USGS basin polygon rather than traced from the gauge, and each tool splits lines at junctions in its own way.
Example 3 — Strahler order on a vector network
import networkx as nx
flowlines = gpd.read_file("nhdplushr_flowlines_esopus.gpkg").to_crs(32618)
rivers = flowlines[flowlines.ftype.isin([460, 558, 334]) & flowlines.intersects(basin.geometry.iloc[0])]
def strahler_from_nodes(lines, id_col="nhdplusid", from_col="fromnode", to_col="tonode"):
graph = nx.MultiDiGraph()
for line_id, a, b in zip(lines[id_col], lines[from_col], lines[to_col]):
graph.add_edge(a, b, key=line_id)
orders = {}
for node in nx.topological_sort(graph):
incoming = sorted((orders[k] for _, _, k in graph.in_edges(node, keys=True)), reverse=True)
level = 1 if not incoming else incoming[0] + int(len(incoming) > 1 and incoming[1] == incoming[0])
for _, _, k in graph.out_edges(node, keys=True):
orders[k] = level
return pd.Series(orders)
start = time.perf_counter()
for label, subset in (("all flowlines", rivers), ("minor divergences removed", rivers[rivers.divergence != 2])):
computed = strahler_from_nodes(subset)
compared = subset.set_index("nhdplusid").assign(computed=computed)
main = compared[compared.divergence != 2]
agree = main.computed == main.streamorde
print(f"{label}: {len(subset)} lines, {time.perf_counter() - start:.2f} s; agrees with NHDPlus streamorde on "
f"{agree.mean():.1%} of lines, {main.length[agree].sum() / main.length.sum():.1%} of length; "
f"max computed {int(main.computed.max())}, max streamorde {int(main.streamorde.max())}")
print(pd.crosstab(main.computed, main.streamorde).to_string())
print("divergence counts", rivers.divergence.value_counts().to_dict())
print("gauge reach", rivers.loc[rivers.gnis_name.eq("Esopus Creek")].sort_values("totdasqkm").tail(1)[["totdasqkm", "streamorde"]].to_dict("records"))all flowlines: 1948 lines, 0.01 s; agrees with NHDPlus streamorde on 91.5% of lines, 91.6% of length; max computed 12, max streamorde 6
minor divergences removed: 1925 lines, 0.03 s; agrees with NHDPlus streamorde on 99.4% of lines, 99.8% of length; max computed 5, max streamorde 6
streamorde 1 2 3 4 5 6
computed
1 1051 0 0 2 1 1
2 7 458 0 0 0 0
3 0 0 227 0 0 0
4 0 0 0 117 0 0
5 0 0 0 0 61 0
divergence counts {0: 1902, 1: 23, 2: 23}
gauge reach [{'totdasqkm': 1200.45849916, 'streamorde': 5}]
divergence is 0 for ordinary flowlines, 1 for the main path at a split and 2 for the minor path. The single order-6 flowline in NHDPlus is a zero-length piece at the basin edge, and the largest Esopus Creek reach touching the basin — 1,200 km² of drainage, downstream of the gauge — is order 5.
Explanation
Why Strahler order needs a topological order
A link's order depends on every link upstream of it. Visiting nodes in topological order — every node after all nodes that drain into it — guarantees that the orders of incoming links are known when a node is reached. A plain loop over lines in file order does not, and gives wrong orders wherever lines are stored downstream-first.
Why braids inflate orders
Strahler's rule adds one when two equal orders meet. A braid splits one stream into two branches of the same order that meet again downstream, which the rule reads as a confluence of equals. Repeated braids stack those increments down a river, which is how an order-5 river reached 12. NHDPlus handles this by giving divergences their own rules; removing minor paths is the simplest equivalent.
Why the raster and vector orders differ at the top
Raster ordering works cell by cell and vectorising splits lines at junctions and direction changes, so segment counts depend on the tool. Orders themselves agree where the network does: both DEM routes found the same 118 second-order lines and 110.8 km of them.
Why cell corners matter
Vector coordinates from pysheds are the affine transform applied to integer row and column numbers, which places them on cell corners. Inverting the transform returns values like 905.9999999 as often as 906.0000001. Truncating those picks the cell to the left or above, which on a one-cell-wide channel is almost always a zero.
Edge cases or notes
- Lines digitised in both directions break topological ordering; check flow direction attributes or derive them from elevation first.
- Cycles from digitising errors make
nx.topological_sortraiseNetworkXUnfeasible; find them withnx.find_cycle. - Disconnected networks order correctly per component; isolated pieces are all order 1.
- Canals and pipelines join basins artificially; remove them before ordering.
- Clipped networks lose headwaters and understate orders near the clip.
- Shreve magnitude follows from the same traversal by summing instead of taking the maximum.
- Coastal deltas are divergence everywhere; Strahler order is poorly defined there.
Internal links
- Stream order explained: Strahler, Shreve and what they measure — the rules and what they mean
- How to extract a stream network from a DEM — the threshold behind the orders
- How to calculate flow accumulation from a DEM in Python — channel masks
- Flow direction explained: D8, D-infinity and MFD compared — why D8 networks have no braids
- How to split a catchment into sub-basins at many outlets — another traversal of the same network
- How to burn a known river network into a DEM — aligning DEM and mapped networks
- How DEM resolution and source change a drainage network — orders at different resolutions
- DEM hydrology explained: from elevation to where water goes — the whole pipeline
FAQ
How do I calculate Strahler stream order in Python?
From a DEM, use pysheds' grid.stream_order or WhiteboxTools' strahler_stream_order on a channel mask. For a vector network, build a directed graph from from/to nodes and propagate orders in topological order.
How do I get Strahler order onto stream line segments?
Sample the order raster at a vertex of each segment, or use WhiteboxTools' raster_streams_to_vector, which writes the order to a STRM_VAL field. pysheds' vertices are cell corners, so floor them with a small tolerance before indexing.
Why are my Strahler orders too high?
Braided or divergent channels count as confluences of equal streams. NHDPlus HR ordered with its 23 minor divergence paths reached order 12; without them the maximum was 5.
How can I calculate stream order for NHDPlus flowlines?
Use fromnode and tonode to build the graph and remove flowlines with divergence equal to 2. The result matched NHDPlus streamorde on 99.4% of Esopus flowlines.
Why does the WhiteboxTools stream shapefile have no CRS?
raster_streams_to_vector wrote no .prj file here. Set the CRS from the DEM with set_crs before measuring or clipping.
Does pysheds or WhiteboxTools give different Strahler orders?
On the same DEM and threshold the orders agreed; segment counts differ because each tool splits lines differently.