Building a NetworkX Graph from a Stream Network

Once a network is a graph, a whole category of question becomes trivial: what drains into this reach, which gauges are upstream of this intake, where does a contaminant released here end up, which confluence separates these two sites. Without the graph each of those is a spatial query repeated until it terminates. This guide covers construction and the traversal patterns that matter, as part of the stream network vectorization and ordering topic within flow routing and stream network extraction.

Prerequisites

• A network with integer from_node and to_node, ideally written as described in vectorizing raster stream networks to GeoPackage.
• networkx 3.x.
• No geometry is needed for the traversals themselves, which is why they stay fast on very large networks.

Core Technique: Two Representations, One Choice

The construction decision that shapes everything after it is what becomes a node.

The junction representation is the direct one, and it is what a traced network already carries. The segment representation, built by networkx.line_graph or by hand, is worth the extra step whenever the question is about accumulating something along segments — a pollutant load, a length, a count of dams — because then “everything upstream” is a single ancestors call returning the segments themselves.

Annotated Code Example

import logging

import geopandas as gpd
import networkx as nx

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


def build_stream_graph(
    links: gpd.GeoDataFrame,
    from_field: str = "from_node",
    to_field: str = "to_node",
    weight_field: str = "length_m",
) -> nx.DiGraph:
    """
    Build a directed junction graph from a stream network.

    Every edge carries the link's attributes, so a traversal can read segment
    properties without a join back to the GeoDataFrame.
    """
    g = nx.DiGraph()
    for idx, row in links.iterrows():
        attrs = {k: row[k] for k in links.columns if k != links.geometry.name}
        attrs["link_index"] = idx
        g.add_edge(int(row[from_field]), int(row[to_field]),
                   weight=float(row.get(weight_field, 0.0)), **attrs)

    log.info("Graph: %d nodes, %d edges", g.number_of_nodes(), g.number_of_edges())

    # --- Validate. These three checks are what a graph makes cheap, and each
    # of them catches a distinct defect in the tracing step upstream. ---
    if not nx.is_directed_acyclic_graph(g):
        cycle = nx.find_cycle(g)
        log.error("Graph contains a cycle (%d edges, e.g. %s) — the flow "
                  "direction grid was not fully conditioned", len(cycle), cycle[:3])

    outlets = [n for n in g.nodes if g.out_degree(n) == 0]
    sources = [n for n in g.nodes if g.in_degree(n) == 0]
    components = nx.number_weakly_connected_components(g)
    log.info("%d outlet(s), %d source(s), %d weakly connected component(s)",
             len(outlets), len(sources), components)

    if components != len(outlets):
        log.warning("Component count (%d) does not match outlet count (%d) — "
                    "the topology is fragmented, most likely from coordinate "
                    "matching instead of integer node ids", components, len(outlets))

    distributaries = [n for n in g.nodes if g.out_degree(n) > 1]
    if distributaries:
        log.warning("%d node(s) have more than one downstream link; D8 routing "
                    "cannot produce these", len(distributaries))

    return g


def upstream_of(g: nx.DiGraph, node: int) -> set:
    """Every node that drains through `node`, excluding it."""
    return nx.ancestors(g, node)


def downstream_path(g: nx.DiGraph, node: int) -> list:
    """The single flow path from `node` to its outlet."""
    path = [node]
    cur = node
    while True:
        succ = list(g.successors(cur))
        if not succ:
            return path
        if len(succ) > 1:
            # A well-formed D8 network never reaches here; if it does, take the
            # highest-order branch so the caller gets the main stem.
            succ.sort(key=lambda n: -g[cur][n].get("strahler", 0))
            log.warning("Node %d has %d downstream links; following the "
                        "highest-order branch", cur, len(succ))
        cur = succ[0]
        path.append(cur)


def upstream_length_km(g: nx.DiGraph, node: int) -> float:
    """Total channel length draining through a node."""
    ups = nx.ancestors(g, node) | {node}
    total = sum(d["weight"] for u, v, d in g.edges(data=True) if u in ups and v in ups)
    return total / 1000.0


def confluence_between(g: nx.DiGraph, a: int, b: int) -> int | None:
    """The first junction where the flow paths from a and b meet."""
    path_a = downstream_path(g, a)
    seen = set(path_a)
    for node in downstream_path(g, b):
        if node in seen:
            return node
    return None


# --- Example usage ---
# net = gpd.read_file("basin_streams.gpkg", layer="streams")
# g = build_stream_graph(net)
# print(len(upstream_of(g, 4211)), "nodes upstream of the intake")
# print(f"{upstream_length_km(g, 4211):.1f} km of channel upstream")
# print("gauges share confluence", confluence_between(g, 1180, 2904))

Parameter Reference

The complexity column is the reason to build the graph at all. An upstream query answered by repeated spatial selection is quadratic in the worst case and involves geometry at every step; the same query on a graph touches only the nodes that are actually upstream.

Graph memory is dominated by whatever is attached to the edges. Keeping geometry out of the graph is the difference between a network that fits comfortably and one that does not.

Propagating an attribute downstream

The pattern that justifies building the graph at all is attribute propagation: a value attached to one segment has to be accumulated over everything downstream of it, or everything upstream. Both are two lines once the graph exists.

To accumulate an upstream quantity — contributing area, an impervious fraction, a pollutant load — walk the graph in topological order and add each node’s own value to the running total of everything that drains into it. networkx.topological_sort gives the order directly, and because the graph is acyclic the traversal touches each edge exactly once. The same loop run on the reversed graph propagates downstream instead, which is what a travel-time or a dilution calculation needs.

The reason to do this on the graph rather than in the GeoDataFrame is that the GeoDataFrame has no notion of order. A groupby cannot express “everything upstream of this row”, so the pandas version of the same operation is a loop that re-scans the frame per segment — quadratic where the graph traversal is linear. On a continental network that is the difference between a query that returns and one that does not.

When the graph should be rebuilt rather than cached

A serialised graph is a snapshot of a topology, and the topology changes whenever the network is re-extracted at a different threshold, re-conditioned, or clipped to a different extent. Caching a pickled graph beside a GeoPackage that has since been regenerated is the graph-layer equivalent of a stale checkpoint: every query answers correctly against a network that no longer exists. Rebuild from the file, which takes seconds even at scale, or key the cache on the same provenance hash the GeoPackage carries.

Gotchas and Edge Cases

• Nodes matched on coordinates. Two endpoints a millimetre apart become two nodes and the graph fragments. Assign integer node ids during tracing and never re-derive them from geometry.
• A cycle in the graph. Always means the flow direction grid contains one, which means conditioning did not complete. It is not a graph problem and cannot be fixed at this layer.
• MultiDiGraph when a DiGraph is meant. Two links between the same pair of nodes are usually a tracing bug, and a DiGraph collapses them silently. Build a DiGraph and check the edge count against the link count.
• Attribute name collisions. Copying every column onto the edge overwrites weight if the network already has such a column. Set weight last, or namespace the copied attributes.
• Assuming one outlet. A raster clipped to a rectangle usually contains several partial basins, each with its own outlet. Multiple components are correct there; the check is whether the count matches the outlet count.
• Holding geometry in the graph. Attaching shapely objects to every edge multiplies memory for no benefit — the traversals never touch geometry. Keep an index back into the GeoDataFrame instead.

Related Topics

• Stream Network Vectorization and Ordering — the parent topic, and where the topology this graph needs comes from
• Vectorizing Raster Stream Networks to GeoPackage — the file format and schema that survive the round trip into this graph
• Computing Strahler Stream Order in Python — an ordering that a graph traversal computes naturally
• Nested Catchment Delineation — the catchment hierarchy this graph is usually built to support