Stream Order Explained: Strahler, Shreve and What They Measure

On this page

Problem statement

"A fifth-order stream" sounds like a property of a river. It is a property of a network drawn at a particular level of detail, counted with a particular rule. Strahler order rises only where two equal orders meet; Shreve magnitude adds up every headwater; Horton order relabels the main stem with the highest order. Computed on a DEM network, all three change with the accumulation threshold, and a published order from a mapped network may not match any of them.

Measured on the Esopus Creek catchment above Coldbrook, New York, with a stream network derived from the 10 m 3DEP DEM at a 0.5 km² threshold, and against the USGS NHDPlus High Resolution network:

At the outlet, Strahler order was 5, Shreve magnitude 250 and Horton order 5.
Inside the basin polygon the network had 253 first-order links, 114 second-order, 69 third-order, 47 fourth-order and 22 fifth-order.
The same river's Strahler order at the outlet ranged from 8 to 4 as the threshold went from 0.05 km² to 10 km².
NHDPlus HR also labels Esopus Creek at Coldbrook as order 5; first-order streams made up 443.6 km of its 805.7 km in the basin (extracting a stream network compares the two).

Quick answer

rule                 at a junction of orders a and b                 at the outlet it counts
Strahler             max(a, b), plus 1 if a == b                     levels of equal-order joining
Shreve               a + b                                           headwater sources upstream
Horton               Strahler, then the main stem gets the maximum   the trunk as one ordered path

With WhiteboxTools, given a D8 pointer and a stream raster:

import whitebox

wbt = whitebox.WhiteboxTools()
wbt.set_max_procs(1)                      # one thread: parallel runs give slightly different results
wbt.set_verbose_mode(False)
wbt.strahler_stream_order("d8.tif", "streams.tif", "strahler.tif", zero_background=True)
wbt.shreve_stream_magnitude("d8.tif", "streams.tif", "shreve.tif", zero_background=True)
wbt.horton_stream_order("d8.tif", "streams.tif", "horton.tif", zero_background=True)

State the threshold and the rule whenever you report an order.

The same network, three numbering rules: Strahler counts levels, Shreve counts sources, Horton follows the trunk.

Links and junctions, then Strahler, Shreve and Horton

1. Build the network as links and junctions

Stream order is defined on a network of links — stretches between a source, a junction or the outlet. On a raster, WhiteboxTools' stream_link_identifier numbers each link; a D8 pointer defines which link drains into which. The network here came from a 0.5 km² accumulation threshold on the least-cost breached 10 m DEM: 499 segments and 442.6 km of channel, as extracting a stream network measured it.

2. Strahler order: levels of branching

Every source link is order 1. Where two links of the same order n meet, the link below is order n + 1; where different orders meet, the link below takes the higher. The Esopus network reached order 5 at the outlet. Strahler order grows slowly — it takes at least two order-4 streams joining to make an order 5 — so it compresses a river's size into a handful of classes.

3. Shreve magnitude: sources upstream

Every source link is magnitude 1, and every junction adds the magnitudes that meet. At the gauge the Shreve magnitude was 250: the number of first-order links upstream of it. Counting links inside the NLDI basin polygon found 253, because the polygon's boundary does not follow the DEM's divides exactly. Magnitude therefore scales with the number of headwaters and, roughly, with drainage area, which makes it more useful than Strahler order for anything proportional to size.

4. Horton order: the trunk carries one label

Horton's scheme starts from Strahler order, then traces the main stream upstream from each junction — following the longer or larger branch — and gives that whole path the higher order. At the outlet it was 5, as for Strahler; upstream, the main stem of Esopus Creek keeps order 5 all the way to its source instead of dropping through 4, 3, 2 and 1.

5. Count links per order

A DEM-derived network is not a natural one, and its ratios do not have to match.

Strahler's rule produces a characteristic geometric series: many first-order links, fewer second-order, and so on.

order   links   cells
1        254   19,816
2        118    9,456
3         64    4,417
4         50    3,364
5         21    1,059

The ratio of links from one order to the next — the bifurcation ratio — was 2.2, 1.7, 1.5 and 2.1. Natural networks typically show ratios between 3 and 5; a DEM network with a single threshold, a clipped headwater area and one basin outlet need not.

6. Expect order to change with the threshold

A lower threshold adds first-order links at the tips of the network, which creates more junctions of equal order downstream. Maximum Strahler order in the basin was 8 at 0.05 km², 7 at 0.1 km², 6 at 0.25 km² and 0.5 km², 5 at 1 and 2 km², and 4 at 5 and 10 km². pysheds' own ordering at 0.5 km² put a single cell at order 6 where WhiteboxTools' maximum was 5 — a reminder that small differences in how tools resolve junctions can shift an order at the top of the network.

7. Compare with mapped networks carefully

NHDPlus HR assigns order 5 to Esopus Creek near Coldbrook, matching the DEM network at 0.5–2 km² thresholds. Its own lengths by order inside the basin were 443.6 km of order 1, 191.1 km of order 2, 87.0 km of order 3, 55.0 km of order 4 and 29.0 km of order 5. A mapped network's order reflects its cartographic detail, not a physical threshold, so agreement at one threshold is a calibration, not a confirmation.

8. Use order for classes, magnitude for size

Strahler order is good for grouping streams — sampling designs, habitat classes, symbology. It is poor as a size measure because one order can span a large range of drainage areas. Where size matters, use accumulation or Shreve magnitude.

253 first-order links and 22 fifth-order ones: most of any network is headwater.

One network ordered three ways

The data used on this page

The setup builds every input this page reads, with the same code as DEM hydrology explained: the basin of Esopus Creek above Coldbrook, New York (USGS 01362500), from the USGS NLDI; USGS 3DEP 1/3 arc-second elevation, resampled bilinearly to UTM 18N at 10 m over the basin plus 1 km (esopus_3dep13_utm_basin.tif); the 10 m DEM conditioned by least-cost breaching in WhiteboxTools (esopus_breached.tif); NHDPlus HR flowlines from the USGS National Map.

import os
import time
from pathlib import Path

import geopandas as gpd
import numpy as np
import rasterio
import requests
from rasterio.merge import merge
from rasterio.warp import Resampling, reproject, transform_bounds


def get(url, **kw):
    """GET with patient retries: the USGS services answer 504 when busy, and recover."""
    for attempt in range(6):
        try:
            r = requests.get(url, timeout=600, **kw)
            r.raise_for_status()
            return r
        except requests.RequestException:
            if attempt == 5:
                raise
            time.sleep(30 * (attempt + 1))


# basins from the USGS NLDI
for site in ('01362500',):
    if not os.path.exists(f"nldi_basin_{site}.geojson"):
        Path(f"nldi_basin_{site}.geojson").write_text(get(f"https://api.water.usgs.gov/nldi/linked-data/nwissite/USGS-{site}/basin?f=json").text)
basin = gpd.read_file("nldi_basin_01362500.geojson")     # Esopus Creek at Coldbrook
W, S, E, N = basin.total_bounds + np.array([-0.03, -0.03, 0.03, 0.03])

TILES = {"3dep13": ([f"/vsicurl/https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/13/TIFF/current/{t}/USGS_13_{t}.tif"
                      for t in ("n43w075", "n42w075")], 10),
         "glo30": ([f"/vsicurl/https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_{t}_00_W075_00_DEM/"
                    f"Copernicus_DSM_COG_10_{t}_00_W075_00_DEM.tif" for t in ("N42", "N41")], 30)}


def dem_box(name, buffer_m, out):
    """UTM 18N grid of whole cells around the basin plus buffer_m, bilinear from the source tiles, no nodata."""
    if os.path.exists(out):
        return
    paths, res = TILES[name]
    srcs = [rasterio.open(p) for p in paths]
    mosaic, tf = merge(srcs, bounds=(W - 0.1, S - 0.1, E + 0.1, N + 0.1))
    bx = basin.to_crs(32618).total_bounds
    x0, y0 = np.floor((bx[0] - buffer_m) / res) * res, np.floor((bx[1] - buffer_m) / res) * res
    x1, y1 = np.ceil((bx[2] + buffer_m) / res) * res, np.ceil((bx[3] + buffer_m) / res) * res
    w, h = int((x1 - x0) / res), int((y1 - y0) / res)
    grid_tf = rasterio.transform.from_origin(x0, y1, res, res)
    grid = np.full((h, w), -9999, "float32")
    reproject(mosaic[0], grid, src_transform=tf, src_crs=srcs[0].crs, dst_transform=grid_tf, dst_crs="EPSG:32618",
              resampling=Resampling.bilinear, src_nodata=srcs[0].nodata, dst_nodata=-9999)
    with rasterio.open(out, "w", driver="GTiff", height=h, width=w, count=1, dtype="float32", crs="EPSG:32618",
                       transform=grid_tf, nodata=-9999, compress="deflate") as dst:
        dst.write(grid, 1)


dem_box("3dep13", 1000, "esopus_3dep13_utm_basin.tif")   # 3DEP 1/3 arc-second at 10 m, basin + 1 km
if not os.path.exists("esopus_breached.tif"):            # least-cost breaching, 20 cells, one thread so it repeats
    import whitebox
    wbt = whitebox.WhiteboxTools()
    wbt.set_verbose_mode(False)
    wbt.set_max_procs(1)
    wbt.breach_depressions_least_cost(os.path.abspath("esopus_3dep13_utm_basin.tif"), os.path.abspath("esopus_breached.tif"),
                                      dist=20, fill=True)
if not os.path.exists("nhdplushr_flowlines_esopus.gpkg"):  # NHDPlus HR flowlines, paged
    feats, off = [], 0
    while True:
        page = get("https://hydro.nationalmap.gov/arcgis/rest/services/NHDPlus_HR/MapServer/3/query", params={
            "geometry": f"{W},{S},{E},{N}", "geometryType": "esriGeometryEnvelope", "inSR": 4326,
            "spatialRel": "esriSpatialRelIntersects", "outFields": "*", "outSR": 4326, "f": "geojson",
            "resultOffset": off, "resultRecordCount": 2000}).json()
        feats += page["features"]
        if len(page["features"]) < 2000:
            break
        off += len(page["features"])
    flowlines = gpd.GeoDataFrame.from_features(feats, crs=4326)
    flowlines.columns = [c.lower() if c != "geometry" else c for c in flowlines.columns]
    flowlines.to_file("nhdplushr_flowlines_esopus.gpkg")

for f in ['nldi_basin_01362500.geojson', 'esopus_3dep13_utm_basin.tif', 'esopus_breached.tif', 'nhdplushr_flowlines_esopus.gpkg']:
    if f.endswith(".tif"):
        with rasterio.open(f) as r:
            print(f"{f}: {r.width:,} x {r.height:,} cells of {r.res[0]:.6g} {'m' if r.crs.is_projected else 'degrees'}")
    elif f.endswith(".zip"):
        print(f"{f}: {os.path.getsize(f) / 1e6:.1f} MB")
    else:
        print(f"{f}: {len(gpd.read_file(f)):,} features")

nldi_basin_01362500.geojson: 1 features
esopus_3dep13_utm_basin.tif: 3,471 x 2,408 cells of 10 m
esopus_breached.tif: 3,471 x 2,408 cells of 10 m
nhdplushr_flowlines_esopus.gpkg: 4,739 features

Example 1 — derive a network and order it three ways

import os

import numpy as np
import rasterio
import whitebox

os.makedirs("order", exist_ok=True)
wbt = whitebox.WhiteboxTools()
wbt.set_max_procs(1)                      # one thread: parallel runs give slightly different results
wbt.set_verbose_mode(False)
wbt.set_working_dir(os.path.abspath("order"))
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)   # 0.5 km2 on a 10 m grid
wbt.strahler_stream_order("d8.tif", "streams.tif", "strahler.tif", zero_background=True)
wbt.shreve_stream_magnitude("d8.tif", "streams.tif", "shreve.tif", zero_background=True)
wbt.horton_stream_order("d8.tif", "streams.tif", "horton.tif", zero_background=True)
wbt.stream_link_identifier("d8.tif", "streams.tif", "links.tif", zero_background=True)


def read(name):
    with rasterio.open(f"order/{name}.tif") as src:
        return src.read(1)


acc, strahler, shreve, horton, links = (read(n) for n in ("acc", "strahler", "shreve", "horton", "links"))

from pyproj import Transformer

with rasterio.open("order/acc.tif") as src:
    transform, crs = src.transform, src.crs
gx, gy = Transformer.from_crs("EPSG:4269", crs, always_xy=True).transform(-74.2701944, 42.0144722)   # Coldbrook gauge
c0, r0 = ~transform * (gx, gy)
r0, c0 = int(r0), int(c0)
window = np.where(strahler[r0 - 15:r0 + 16, c0 - 15:c0 + 16] > 0, acc[r0 - 15:r0 + 16, c0 - 15:c0 + 16], -1)
i = np.unravel_index(window.argmax(), window.shape)
row, col = r0 - 15 + i[0], c0 - 15 + i[1]                     # the largest network cell within 150 m
print(f"outlet ({row}, {col}): drains {acc[row, col] / 1e4:.2f} km2, Strahler {strahler[row, col]:.0f}, "
      f"Shreve {shreve[row, col]:.0f}, Horton {horton[row, col]:.0f}")

outlet (2139, 2163): drains 492.35 km2, Strahler 5, Shreve 250, Horton 5

The outlet is the largest network cell within 150 m of the Coldbrook gauge. The DEM's own lowest network cell, where Esopus Creek leaves the grid below the gauge, drained 513.54 km² (measured in flow accumulation) and was already Strahler order 6 with Shreve magnitude 260, because another tributary joins in between.

Example 2 — links per order and bifurcation ratios

import pandas as pd

from rasterio import features
import geopandas as gpd

with rasterio.open("order/acc.tif") as src:
    transform, crs = src.transform, src.crs
basin = gpd.read_file("nldi_basin_01362500.geojson").to_crs(crs)
inside = features.rasterize(basin.geometry, out_shape=acc.shape, transform=transform).astype(bool)

cells = pd.DataFrame({"link": links[inside & (links > 0)], "order": strahler[inside & (links > 0)]})
per_link = cells.groupby("link").order.max()
table = pd.DataFrame({"links": per_link.value_counts().sort_index(), "cells": cells.groupby("order").size()})
table["bifurcation_ratio"] = (table.links / table.links.shift(-1)).round(2)
print(table.to_string())
print("first-order links:", int(table.links.iloc[0]))

       links  cells  bifurcation_ratio
order                                 
1.0      253  19751               2.22
2.0      114   8956               1.65
3.0       69   4971               1.47
4.0       47   3340               2.14
5.0       22   1057                NaN
first-order links: 253

Example 3 — orders in a mapped network

mapped = gpd.read_file("nhdplushr_flowlines_esopus.gpkg").to_crs(crs)
mapped = gpd.clip(mapped[mapped.ftype.isin([460, 558, 334])], basin)
by_order = mapped.assign(km=mapped.length / 1000).groupby("streamorde").agg(flowlines=("km", "size"), km=("km", "sum"))
print(by_order.round(1).to_string())

            flowlines     km
streamorde                  
1                1066  443.6
2                 458  191.1
3                 232   87.0
4                 125   55.0
5                  66   29.0
6                   1    0.0

streamorde is NHDPlus' Strahler order attribute, computed on its own mapped network.

Why Strahler order rises so slowly

To reach order n + 1, two order-n streams must meet. Each of those needed two order n − 1 streams, and so on. An order-5 outlet therefore needs at least 16 first-order sources arranged just so, and in practice hundreds; the Esopus had 250. Orders grow with the logarithm of network size, which is why rivers the size of the Mississippi are only order 10 to 12 on typical maps.

Why Shreve magnitude equals the number of sources

Each source contributes 1, and magnitudes only ever add at junctions. Nothing is lost or created along the way, so the value at any link is the count of sources upstream of it. That makes magnitude additive, like drainage area, where Strahler order is not.

Why thresholds move orders

Adding tiny headwater links does not change the main river, but each one can turn an unequal junction into an equal one and push orders up by one all the way down the network. The network's order is as much a statement about the smallest streams included as about the largest river.

Why tools can disagree at the top

Order depends on which links meet at each junction. Where D8 routes two streams into the same cell from different sides, or where a tool treats a cell as a junction slightly differently, the order of one short link can change, and the change propagates downstream. A single cell of order 6 in pysheds against a maximum of 5 in WhiteboxTools is such a case.

Order is a property of where you decided a stream begins.

Braids, clipped networks and main-stem rules

Braided channels and canals create loops that ordering rules do not handle; D8 networks have none, mapped networks may.
Clipped networks lose headwaters outside the clip and understate orders.
Horton order depends on the rule for the main stem — longest path or largest area — which differs between tools.
Vector ordering on mapped lines needs topology: every junction snapped and every line directed downstream.
Hack order numbers streams from the outlet upstream, the reverse of Strahler.
Order classes for sampling should be defined with a fixed threshold across all basins compared.
Stream order is dimensionless; do not average it or use it in regressions as if it were a size.

Internal links

How to calculate Strahler stream order in Python — ordering in practice
How to extract a stream network from a DEM — the thresholds that set the order
How to calculate flow accumulation from a DEM in Python — size measured directly
How to split a catchment into sub-basins at many outlets — another way to structure a network
How DEM resolution and source change a drainage network — resolution effects on networks
DEM hydrology explained: from elevation to where water goes — the pipeline
Flow direction explained: D8, D-infinity and MFD compared — why D8 networks have no loops
How to burn a known river network into a DEM — aligning DEM and mapped networks

FAQ

What is Strahler stream order?

A classification in which headwater streams are order 1 and the order increases by one only where two streams of equal order meet. The Esopus network at a 0.5 km² threshold reached order 5 at Coldbrook.

What is the difference between Strahler order and Shreve magnitude?

Strahler order takes the maximum at junctions and rises only on equal orders; Shreve magnitude adds the values. Shreve magnitude at the Esopus outlet was 250, the number of headwater links.

Why does my river's stream order change?

Because the network's level of detail changed. The Esopus outlet was order 8 with a 0.05 km² threshold and order 4 with 10 km².

Is stream order a good measure of river size?

Only roughly. One Strahler order covers a wide range of drainage areas; use accumulation or Shreve magnitude when size matters.

What is a bifurcation ratio?

The number of links of one order divided by the number of the next order. For the Esopus DEM network the ratios were 2.2, 1.7, 1.5 and 2.1.

Do DEM-derived orders match NHDPlus stream orders?

They can at a suitable threshold: NHDPlus HR gives Esopus Creek at Coldbrook order 5, as did the DEM network at 0.5–2 km². Agreement depends on both networks' level of detail.