Watershed Modeling Python for hydrology automation

Computing Strahler Stream Order in Python

Strahler order compresses a network’s branching structure into a single integer per segment, and almost every downstream use — selecting reaches to map, weighting them for display, classifying them geomorphically — depends on it being computed correctly. The rule is three lines long, and the implementations that get it wrong nearly all fail for the same reason: they process links in the wrong sequence. This guide belongs to the stream network vectorization and ordering topic within flow routing and stream network extraction.

Prerequisites

  • A vector stream network with valid topology: every segment carrying a from_node and a to_node, and every node shared exactly by the segments that meet there.
  • networkx (optional but convenient) or nothing beyond the standard library for the relaxation approach.
  • The accumulation threshold that produced the network, because the result is meaningless without it.

If the topology is not already present, produce it first — a network whose segments merely touch geometrically will not order correctly, and the failure is silent.

Core Technique: The Rule and the Sequence

The Strahler rule at any junction is:

  • A segment with no upstream segments is order 1.
  • Otherwise, take the two highest incoming orders. If they are equal, the segment is that value plus one. If not, it is the maximum.

The rule is local. What makes an implementation correct is the sequence: a segment can only be assigned once every segment upstream of it has been assigned. Processing in arbitrary order produces a network where order propagates part-way and stops, and the symptom is a large river carrying order 2.

The Rule at Three Junction Types Two order-2 links meeting produce order 3. An order-1 joining an order-3 leaves it order 3. Three order-2 links meeting produce order 3, because the rule considers only the two highest incoming orders. 2 2 3 equal orders meet → increment 3 1 3 unequal orders meet → take the maximum 2 2 2 3 three equal orders meet → still only +1 the rule reads the two highest incoming orders and nothing else

Annotated Code Example

The implementation below uses a relaxation loop rather than a graph library, so it works on any table carrying from_node and to_node — including one read straight from a database.

python 
import logging
from collections import defaultdict

import geopandas as gpd

log = logging.getLogger(__name__)
logging.basicConfig(level=logging.INFO, format="%(levelname)s %(message)s")


def assign_strahler(
    links: gpd.GeoDataFrame,
    from_field: str = "from_node",
    to_field: str = "to_node",
    out_field: str = "strahler",
) -> gpd.GeoDataFrame:
    """
    Assign Strahler order to every link in a topologically valid network.

    Parameters
    ----------
    links      : One row per stream segment, carrying node identifiers.
    from_field : Column holding the upstream node id.
    to_field   : Column holding the downstream node id.
    out_field  : Column to write the order into.

    Returns
    -------
    A copy of `links` with the order column added.
    """
    out = links.copy()

    # --- Index links by the node they leave, and by the node they enter.
    # A link's upstream neighbours are exactly the links that ENTER its
    # from_node — this is the lookup the whole algorithm rests on. ---
    entering = defaultdict(list)
    for idx, row in out.iterrows():
        entering[row[to_field]].append(idx)

    upstream_of = {
        idx: entering.get(row[from_field], [])
        for idx, row in out.iterrows()
    }

    order: dict = {}
    sources = [i for i, ups in upstream_of.items() if not ups]
    log.info("Network has %d links and %d source links", len(out), len(sources))
    if not sources:
        raise ValueError("no source links found — the topology is cyclic or the "
                         "from/to fields are reversed")

    # --- Relaxation: repeatedly resolve every link whose upstream neighbours
    # are all resolved. Converges in as many passes as the maximum order. ---
    pending = set(out.index)
    passes = 0
    while pending:
        passes += 1
        progressed = False
        for idx in list(pending):
            ups = upstream_of[idx]
            if any(u not in order for u in ups):
                continue
            if not ups:
                order[idx] = 1
            else:
                incoming = sorted((order[u] for u in ups), reverse=True)
                # The rule reads the top TWO only. Three equal orders meeting
                # still increments by one, not by two.
                order[idx] = (incoming[0] + 1
                              if len(incoming) > 1 and incoming[0] == incoming[1]
                              else incoming[0])
            pending.discard(idx)
            progressed = True
        if not progressed:
            # Every remaining link waits on another remaining link: a cycle.
            log.error("Ordering stalled after %d pass(es) with %d link(s) "
                      "unresolved — the network contains a directed cycle",
                      passes, len(pending))
            for idx in pending:
                order[idx] = -1
            break

    out[out_field] = [order[i] for i in out.index]
    resolved = out[out[out_field] > 0]
    log.info("Resolved in %d pass(es); max order %d over %d links",
             passes, int(resolved[out_field].max()), len(resolved))
    for o, grp in resolved.groupby(out_field):
        log.info("  order %d: %4d links, %8.1f km",
                 o, len(grp), grp.geometry.length.sum() / 1000)
    return out


# --- Example usage ---
# net = gpd.read_file("basin_streams.gpkg", layer="streams")
# net = assign_strahler(net)
# net.to_file("basin_streams_ordered.gpkg", layer="streams", driver="GPKG")

The relaxation loop converges in at most max_order passes, which for a real basin is under ten. On a network of a hundred thousand links it runs in a couple of seconds, and the loop is easier to reason about than a recursive traversal that can exhaust the stack on a long main stem.

Parameter Reference

Consideration Value Effect
Accumulation threshold 100–2 000 cells The single largest control on the maximum order — halving it typically adds one to two orders
Junction arity 2 for DEM networks Higher arity appears after snapping or simplification; the rule handles it, but check it is intentional
Node tolerance Exact integer ids Node identity must be exact. Coordinate-matched nodes with a tolerance produce silent fragmentation
Cycle handling Fail loudly A cycle means the direction grid was not fully conditioned; ordering cannot be salvaged

Worked Example: Reading the Distribution

A 10 m DEM over a 1 480 km² basin at a 500-cell threshold:

Order Links Total length (km) Share of length
1 2 841 1 964 52 %
2 704 918 24 %
3 178 449 12 %
4 41 236 6 %
5 11 148 4 %
6 3 62 2 %

The link count falls by roughly a factor of four per order, which is the bifurcation ratio and the signature of a healthy dendritic network. Ratios far from that range are worth investigating: a ratio near two suggests a strongly trellised or structurally controlled network, and a ratio above eight usually means the extraction produced spurious short first-order stubs.

The Bifurcation Ratio Is a Quality Check Link counts by order plotted on a logarithmic scale fall in a near-straight line, giving a bifurcation ratio of about 4.0. Annotations mark that a ratio near 2 suggests structural control and a ratio above 8 suggests spurious headwater stubs. 2 841 704 178 41 11 3 1 2 3 4 5 6 Strahler order links, log scale ratio ≈ 4.0 healthy near 2: trellised above 8: stubs

The relaxation converges in as many passes as the maximum order, which is why it is fast on any real network however large. The pass count is also a free diagnostic: it should equal the maximum order exactly.

The Relaxation Converges in max-order Passes Of 3778 links, 2841 resolve on the first pass, 704 on the second, 178 on the third, 41 on the fourth, 11 on the fifth and 3 on the sixth, after which nothing remains. The pass count equals the maximum Strahler order of six. 2 841 704 178 41 11 3 pass 1 2 3 4 5 6 More passes than the maximum order means the loop is stalling — look for a cycle. links resolved

Gotchas and Edge Cases

  • Processing in file order. Order propagates partially and stops. Symptom: the main stem carries a low order while its tributaries carry higher ones — an impossible configuration that is diagnostic of exactly this bug.
  • Reversed from/to fields. The relaxation finds no sources and raises immediately, which is the good case. Where the network is symmetric enough to have apparent sources at both ends, it produces an order field that is quietly wrong.
  • Coordinate-matched nodes. Matching endpoints by coordinate with a tolerance creates near-duplicate nodes at junctions, splitting one confluence into two, which suppresses the increment. Use integer node ids assigned during tracing.
  • Comparing orders across thresholds. Two networks at different thresholds have incomparable orders. Store the threshold with the layer.
  • Assuming order tracks discharge. It does not, and it is not meant to. An order-4 segment in a wet basin and an order-4 segment in an arid one differ in discharge by orders of magnitude. Shreve magnitude is the ordering that correlates with discharge.
  • Order on a network with distributaries. Braided or deltaic networks have segments with two downstream neighbours, which the Strahler rule does not define. Either resolve the braid to a single path first or use a scheme designed for it.