Watershed Modeling Python for hydrology automation

Stream Network Vectorization and Ordering

A thresholded accumulation raster is a picture of a stream network, not a stream network. It has no segments, no junctions, no direction and no way to answer the question every downstream consumer asks first — what drains into this reach. As part of the flow routing and stream network extraction workflow, vectorization and ordering are the steps that convert that picture into a topological object: a set of directed segments, each knowing its downstream neighbour, each carrying an order that says where it sits in the branching hierarchy. Everything from nested catchment delineation to reach-scale routing depends on getting this conversion right.

The step is routinely done badly because a generic raster-to-vector polygoniser appears to work. It produces lines, the lines look like rivers, and the failure only surfaces two stages later when a graph traversal reports that a river system has forty-seven disconnected components.

Prerequisites and Environment Setup

Vectorization consumes the outputs of the routing stage, so it assumes both a flow direction grid and an accumulation grid derived from the same conditioned DEM.

bash 
conda create -n streamnet python=3.11
conda activate streamnet
conda install -c conda-forge rasterio=1.3 geopandas=0.14 shapely=2.0 \
    scikit-image=0.22 networkx=3.2 numpy scipy pyogrio
Input Requirement Notes
Flow direction grid Same grid as the DEM, D8 encoding The tracing step walks this, not the accumulation raster
Flow accumulation grid Same grid, cells or area units Only used to apply the threshold
Threshold Calibrated, not defaulted See stream threshold tuning
CRS Projected, metres Segment lengths are meaningless in degrees

One constraint deserves emphasis: the direction grid and the accumulation grid must come from the same conditioning run. Mixing an accumulation raster from one fill pass with directions from another produces segments that terminate in the middle of nowhere, and the symptom — isolated one-cell links scattered through the network — looks like a vectorization bug rather than an input mismatch.

Mechanics: From Cells to Segments

The conversion has three distinct phases, and each has a characteristic way of going wrong.

Mask, Skeleton, Segments A thresholded accumulation mask is several cells wide in places. Thinning reduces it to a one-cell skeleton. Tracing along the flow direction grid converts the skeleton into directed line segments that meet exactly at junction nodes. 1. thresholded mask two cells wide in places a polygoniser makes blobs here 2. thinned skeleton one cell wide everywhere but still just pixels 3. traced segments three links meeting at one node directed, connected, orderable The junction node is shared by all three segments — that shared identity is the whole point of the exercise.

Thinning reduces a mask that may be several cells wide to a one-cell skeleton. A D8 network above threshold is usually already one cell wide, but a network derived from D-infinity routing patterns or multiple flow direction methods is not, because those routers spread accumulation across neighbours. Skimping on thinning is what produces the blobs at confluences that later read as topology errors.

Classification labels every skeleton cell by how many skeleton neighbours it has. One neighbour makes it a source or an outlet, two make it an interior link cell, three or more make it a junction. That count is the entire basis of segmentation.

Tracing walks from each source or junction along the flow direction grid until it reaches the next junction or the network edge, emitting one line per link. Tracing along directions rather than along the skeleton is what guarantees segments meet exactly.

Why order matters, and which order

Stream order is the compact answer to “how far into the network am I”, and the two common definitions answer different questions.

Property Strahler Shreve
Rule at a confluence Increment only when two equal orders meet Sum the incoming orders
Source segment 1 1
First-order plus second-order Stays 2 Becomes 3
Range on a large basin Typically 1–8 1 to several thousand
Correlates with Branching hierarchy, geomorphic classification Upstream magnitude, and thus discharge
Stable under threshold change No — a lower threshold adds sources and raises order No, and it moves further
Typical use Selecting reaches for mapping or display Weighting reaches for routing or load estimation

Neither is stable under a change of accumulation threshold, which is the fact most often forgotten. Lowering the threshold adds headwater links, and every added pair of first-order links promotes the segment below them. An order value is therefore only meaningful alongside the threshold that produced it, and comparing orders between two networks extracted at different thresholds is meaningless.

Strahler and Shreve on the Same Network A branching network of seven segments. Strahler orders rise 1, 1, 2, 1, 2, 2, 3 down the network, incrementing only where equal orders meet. Shreve orders rise 1, 1, 2, 1, 3, 1, 4, incrementing at every confluence. S1 / M1 S1 / M1 S2 / M2 S1 / M1 S1 / M1 S2 / M2 S3 / M4 S = Strahler, M = Shreve magnitude Where they diverge Two order-1 links meeting: Strahler → 2, Shreve → 2. They agree. An order-1 joining an order-2: Strahler stays 2, Shreve becomes 3. Strahler is deliberately insensitive to small tributaries; Shreve counts them all. Both are threshold-dependent: lowering the accumulation threshold adds sources and shifts every order downstream.

Step-by-Step Workflow

  1. Apply the threshold. Produce a boolean mask from the accumulation grid using the calibrated value. Do not use a software default here — the reasons are set out in tuning flow accumulation thresholds for ephemeral streams.
  2. Thin to a skeleton. Apply a morphological thinning that preserves connectivity. skimage.morphology.skeletonize is adequate for D8-derived masks; dispersive masks may need a preliminary binary closing to avoid fragmenting.
  3. Classify cells. Count stream neighbours in the eight-connected neighbourhood of every stream cell. One neighbour is a source or outlet, two an interior cell, three or more a junction.
  4. Trace links. From each source and from each cell immediately downstream of a junction, follow the flow direction grid until reaching a junction or leaving the network. Record the cell path.
  5. Build node identities. Assign an integer id to each unique junction, source and outlet cell, and attach from_node and to_node to every link.
  6. Assign order. Process links in upstream-to-downstream order, applying the Strahler or Shreve rule at each junction.
  7. Export. Write to GeoPackage with the topology fields intact, and record the threshold in the layer metadata.

Production-Ready Code

The function below performs the whole conversion and returns a GeoDataFrame with topology and both orders attached. It traces along the flow direction grid rather than along the skeleton, which is what makes the resulting network connected by construction.

python 
import logging
from collections import defaultdict

import geopandas as gpd
import numpy as np
import rasterio
from shapely.geometry import LineString
from skimage.morphology import skeletonize

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

# D8 encoding used by RichDEM / TauDEM: 1=E, 2=NE, 3=N, 4=NW, 5=W, 6=SW, 7=S, 8=SE
D8_OFFSETS = {
    1: (0, 1), 2: (-1, 1), 3: (-1, 0), 4: (-1, -1),
    5: (0, -1), 6: (1, -1), 7: (1, 0), 8: (1, 1),
}


def vectorize_stream_network(
    fdir_path: str,
    facc_path: str,
    threshold_cells: int,
    out_path: str | None = None,
) -> gpd.GeoDataFrame:
    """
    Convert a thresholded accumulation raster into an ordered vector network.

    Parameters
    ----------
    fdir_path       : D8 flow direction raster (encoding above).
    facc_path       : Flow accumulation raster on the same grid.
    threshold_cells : Calibrated channel-initiation threshold, in cells.
    out_path        : Optional GeoPackage destination.

    Returns
    -------
    GeoDataFrame with columns: link_id, from_node, to_node, length_m,
    strahler, shreve, and the segment geometry.
    """
    with rasterio.open(fdir_path) as src:
        fdir = src.read(1)
        transform = src.transform
        crs = src.crs
    with rasterio.open(facc_path) as src:
        facc = src.read(1)

    if fdir.shape != facc.shape:
        raise ValueError("direction and accumulation grids are not on the same grid")

    # --- Step 1-2: threshold, then thin to a one-cell skeleton ---
    mask = facc >= threshold_cells
    log.info("Threshold %d cells selects %d stream cells", threshold_cells, mask.sum())
    skel = skeletonize(mask)
    log.info("Skeletonised to %d cells (%.1f %% of the mask)",
             skel.sum(), 100.0 * skel.sum() / max(1, mask.sum()))

    rows, cols = np.nonzero(skel)
    stream = set(zip(rows.tolist(), cols.tolist()))

    def downstream(rc):
        """The cell this one drains to, or None at the network edge."""
        off = D8_OFFSETS.get(int(fdir[rc]))
        if off is None:
            return None
        nxt = (rc[0] + off[0], rc[1] + off[1])
        if not (0 <= nxt[0] < skel.shape[0] and 0 <= nxt[1] < skel.shape[1]):
            return None
        return nxt if nxt in stream else None

    # --- Step 3: an upstream count per cell identifies sources and junctions ---
    upstream_count = defaultdict(int)
    for rc in stream:
        d = downstream(rc)
        if d is not None:
            upstream_count[d] += 1

    sources = [rc for rc in stream if upstream_count[rc] == 0]
    junctions = {rc for rc in stream if upstream_count[rc] > 1}
    log.info("Network has %d sources and %d junctions", len(sources), len(junctions))

    # --- Step 4: trace one link from every source and from below every junction ---
    starts = list(sources)
    for j in junctions:
        d = downstream(j)
        if d is not None:
            starts.append(d)

    links, node_ids = [], {}

    def node_id(rc):
        if rc not in node_ids:
            node_ids[rc] = len(node_ids)
        return node_ids[rc]

    for start in starts:
        path, cur = [start], start
        while True:
            nxt = downstream(cur)
            if nxt is None:
                break
            path.append(nxt)
            if nxt in junctions:          # a link ends at the junction it feeds
                break
            cur = nxt
        if len(path) < 2:
            continue
        # cell indices to map coordinates, at cell centres
        coords = [transform * (c + 0.5, r + 0.5) for r, c in path]
        links.append({
            "geometry": LineString(coords),
            "from_node": node_id(path[0]),
            "to_node": node_id(path[-1]),
            "_head": path[0],
            "_tail": path[-1],
        })

    log.info("Traced %d links", len(links))

    # --- Step 5-6: order the links, then apply both ordering rules ---
    by_tail = defaultdict(list)
    for i, ln in enumerate(links):
        by_tail[ln["_head"]].append(i)     # links that begin where another ended

    strahler = [0] * len(links)
    shreve = [0] * len(links)

    # Process in an order where every upstream link is resolved first.
    resolved, pending = set(), list(range(len(links)))
    while pending:
        progressed = False
        still = []
        for i in pending:
            feeders = [j for j in range(len(links))
                       if links[j]["_tail"] == links[i]["_head"]]
            if all(j in resolved for j in feeders):
                if not feeders:
                    strahler[i], shreve[i] = 1, 1
                else:
                    orders = sorted((strahler[j] for j in feeders), reverse=True)
                    # Strahler: increment only when the top two orders are equal
                    strahler[i] = (orders[0] + 1
                                   if len(orders) > 1 and orders[0] == orders[1]
                                   else orders[0])
                    shreve[i] = sum(shreve[j] for j in feeders)
                resolved.add(i)
                progressed = True
            else:
                still.append(i)
        if not progressed:
            # A cycle cannot exist in a correct flow grid; if we get here the
            # direction raster is inconsistent with the accumulation raster.
            log.error("Ordering stalled with %d links unresolved — check that the "
                      "direction and accumulation grids share a conditioning run",
                      len(still))
            for i in still:
                strahler[i], shreve[i] = -1, -1
            break
        pending = still

    for i, ln in enumerate(links):
        ln["strahler"] = strahler[i]
        ln["shreve"] = shreve[i]
        ln.pop("_head"), ln.pop("_tail")

    gdf = gpd.GeoDataFrame(links, crs=crs)
    gdf.insert(0, "link_id", range(len(gdf)))
    gdf["length_m"] = gdf.geometry.length
    log.info("Max Strahler order %d, max Shreve magnitude %d, total length %.1f km",
             gdf["strahler"].max(), gdf["shreve"].max(), gdf["length_m"].sum() / 1000)

    if out_path:
        gdf.to_file(out_path, layer="streams", driver="GPKG")
        log.info("Wrote %s", out_path)
    return gdf


# --- Example usage ---
# net = vectorize_stream_network(
#     fdir_path="basin_d8.tif",
#     facc_path="basin_acc.tif",
#     threshold_cells=500,
#     out_path="basin_streams.gpkg",
# )
# print(net.groupby("strahler")["length_m"].sum() / 1000)

Validation Protocol

A vector network is valid when it is a directed acyclic graph with exactly one outlet per drainage basin, and every check below is a way of testing that claim.

  • Component count. Build a graph from the from_node/to_node pairs and count weakly connected components. It should equal the number of independent basins draining off the raster edge. A count in the dozens means the tracing fell back to skeleton adjacency somewhere.
  • Out-degree. Every node except a basin outlet must have exactly one downstream link. A node with two is a distributary, which D8 cannot produce, so it is a bug.
  • No cycles. networkx.is_directed_acyclic_graph must return true. A cycle means the direction grid contains one, which means conditioning did not complete.
  • Length against the raster. Total vector length should be within a few percent of the stream-cell count times the mean cell traversal distance. A large shortfall means links were dropped at junctions.
  • Order sanity. Maximum Strahler order should be plausible for the basin size — roughly 4 to 6 for a few hundred square kilometres, 7 to 9 for a large river basin. A maximum of 2 over a large basin means junctions were not detected.

The connectivity checks are worth running as a set, because each one catches a different class of defect and a network that passes all four is genuinely usable as a graph.

Four Checks, Four Distinct Defects Component count catches gaps at junctions. Out-degree catches spurious distributaries. Acyclicity catches an unconditioned direction grid. Total length against the raster catches links silently dropped during tracing. components should equal the number of basins catches: gaps at junctions out-degree exactly 1 for every node but the outlet catches: spurious distributaries acyclicity no directed cycle anywhere catches: incomplete conditioning total length within a few % of the cell count catches: links dropped in tracing Only the first three are cheap enough to run on every build; the length check needs the raster, so it belongs in the extraction job rather than in the consumer that reads the GeoPackage. A network passing all four can be traversed, ordered and cut into catchments without further repair.

Common Failure Modes and Optimization

  • Polygonising instead of tracing. rasterio.features.shapes on the stream mask produces polygon boundaries, not centrelines, and the resulting lines are disconnected at every junction. It is the fastest way to get a plausible-looking, unusable network.
  • Thinning a dispersive mask without closing. MFD accumulation above threshold can be patchy; skeletonising it directly fragments the network into hundreds of stubs. A binary closing with a small structuring element before thinning fixes this.
  • Mixed conditioning runs. Direction from one fill, accumulation from another. The ordering loop stalls, which the code above detects and reports rather than silently emitting order −1 values nobody notices.
  • Geographic CRS. Segment lengths in degrees are not a length. Every threshold, every length statistic and every order-versus-length regression becomes nonsense. Reproject first; see coordinate reference system alignment.
  • Quadratic ordering. The feeder lookup in the reference implementation is a linear scan for clarity. On networks above a few thousand links, index the links by their head cell first — the difference is minutes versus hours.
  • Order compared across thresholds. Two networks extracted at 300 and 800 cells cannot have their orders compared. Store the threshold as layer metadata so the comparison is at least visibly wrong rather than invisibly so.

When to Use This vs. Alternatives

Use a published hydrography instead when one exists at adequate resolution and the analysis does not need the network to agree with your DEM. NHD flowlines are cartographically better than anything extracted from a 10 m grid. The catch is that they will not align with your flow directions, which matters as soon as you delineate from them — the reconciliation is covered in snapping stream gauge locations to NHD flowlines.

Stay in raster space when the downstream consumer is itself a raster operation. Computing a distance-to-stream grid or a HAND surface needs the stream mask, not segments, and vectorizing to feed it back into a raster operation adds error for no gain.

Vectorize whenever the question is topological: which reaches are upstream of this point, what is the order distribution, how do these segments connect to a routing model. Those are graph questions, and a raster cannot answer them.

Frequently Asked Questions

What is the difference between Strahler and Shreve stream order?

Strahler order increments only when two segments of equal order meet: two first-order streams make a second-order stream, but a first-order joining a second-order leaves it second-order. Shreve order is additive — it is the count of source links upstream, so every confluence increases it. Strahler is a measure of branching hierarchy and is what most geomorphic literature uses; Shreve is a measure of upstream magnitude and correlates far better with discharge.

Why does my vectorized stream network have gaps at confluences?

Almost always because the segments were traced from a thinned raster mask using cell centres rather than from the flow direction grid. Skeletonisation can leave a diagonal step at a junction that the polygoniser renders as two lines ending one cell apart. Tracing along flow directions instead guarantees each segment terminates on the exact cell its downstream neighbour begins from, which closes the gap by construction.

Does the flow routing algorithm change the extracted stream network?

Yes, in two ways. Dispersive routers such as D-Infinity and MFD spread accumulation across several downslope neighbours, so peak accumulation values are lower and the calibrated threshold must be lower to extract the same network. They also produce accumulation fields that are not strictly single-path, so a thinning step is required before segments can be traced. D8 produces a network that is already one cell wide wherever accumulation exceeds the threshold.