Watershed Modeling Python for hydrology automation

Vectorizing Raster Stream Networks to GeoPackage

A traced network is only useful if it survives the trip to disk with its topology intact, and the format and schema choices decide whether it does. This guide covers the export step — schema, indexing, metadata and the pitfalls that silently destroy connectivity — as part of the stream network vectorization and ordering topic within flow routing and stream network extraction.

Prerequisites

  • A traced network in memory as a GeoDataFrame with from_node, to_node and geometry.
  • geopandas with the pyogrio engine, which is substantially faster for GeoPackage writes than the older Fiona path.
  • A projected CRS. Writing a network in degrees produces a file whose length column is not a length.

Core Technique: A Schema That Preserves Topology

GeoPackage is the right container here for three reasons: it is a single file, it holds multiple layers, and it is a SQLite database, so the topology can be indexed and queried without loading the geometry at all.

The schema matters more than the format. A network layer needs three groups of columns.

Three Column Groups, Three Different Jobs Identity and topology columns — link id, from node, to node — carry the graph. Measured columns — length, order, magnitude, upstream area — carry the hydrology. Provenance columns — threshold, algorithm, DEM checksum — make the layer reproducible and comparable. identity and topology link_id, from_node, to_node integers, never floats indexed, never simplified without these it is a picture measured attributes length_m, strahler, shreve, upstream_area_km2, slope derived, recomputable recomputing needs the raster provenance threshold_cells, router, dem_sha256, extracted_at a separate metadata table without these, nothing compares The third group is the one that gets dropped, and it is the one that turns two extractions of the same basin from an unresolvable disagreement into an explainable one.

Annotated Code Example

python 
import hashlib
import json
import logging
import sqlite3
from datetime import datetime, timezone

import geopandas as gpd

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


def write_stream_geopackage(
    links: gpd.GeoDataFrame,
    out_path: str,
    dem_path: str,
    threshold_cells: int,
    router: str = "D8",
    conditioning: str = "priority-flood fill, breached road crossings",
    layer: str = "streams",
) -> dict:
    """
    Write a traced network to GeoPackage with topology, indexes and provenance.

    Parameters
    ----------
    links           : Traced segments carrying from_node, to_node and geometry.
    out_path        : Destination .gpkg path.
    dem_path        : Source DEM, hashed into the provenance record.
    threshold_cells : The accumulation threshold used for extraction.
    router          : Routing algorithm name, recorded verbatim.
    conditioning    : Free-text description of the conditioning applied.
    layer           : Layer name inside the GeoPackage.
    """
    required = {"from_node", "to_node"}
    missing = required - set(links.columns)
    if missing:
        raise ValueError(f"network is missing topology column(s): {sorted(missing)}")
    if links.crs is None or links.crs.is_geographic:
        raise ValueError("network must be in a projected CRS — lengths in degrees "
                         "are not lengths")

    out = links.copy()

    # --- Node ids must be integers. Float ids survive a round trip through
    # some drivers as 1.0000000001 and stop matching, which fragments the
    # topology in a way that is very hard to see. ---
    for col in ("from_node", "to_node"):
        out[col] = out[col].astype("int64")
    if "link_id" not in out.columns:
        out.insert(0, "link_id", range(len(out)))
    out["link_id"] = out["link_id"].astype("int64")

    if "length_m" not in out.columns:
        out["length_m"] = out.geometry.length

    # --- Write geometry unsimplified. Any generalisation moves endpoints,
    # and a moved endpoint no longer coincides with its neighbour's. ---
    out.to_file(out_path, layer=layer, driver="GPKG", engine="pyogrio")
    log.info("Wrote %d links to %s (layer '%s')", len(out), out_path, layer)

    # --- Index the topology columns so upstream/downstream queries do not
    # scan the table, and record provenance in a sibling table. ---
    dem_hash = _sha256(dem_path)
    meta = {
        "threshold_cells": threshold_cells,
        "router": router,
        "conditioning": conditioning,
        "dem_path": dem_path,
        "dem_sha256": dem_hash,
        "crs": out.crs.to_string(),
        "link_count": int(len(out)),
        "total_length_km": round(float(out["length_m"].sum()) / 1000.0, 3),
        "extracted_at": datetime.now(timezone.utc).isoformat(timespec="seconds"),
    }

    with sqlite3.connect(out_path) as con:
        con.execute(f'CREATE INDEX IF NOT EXISTS idx_{layer}_from '
                    f'ON "{layer}" (from_node)')
        con.execute(f'CREATE INDEX IF NOT EXISTS idx_{layer}_to '
                    f'ON "{layer}" (to_node)')
        con.execute("CREATE TABLE IF NOT EXISTS extraction_metadata "
                    "(key TEXT PRIMARY KEY, value TEXT)")
        con.executemany("INSERT OR REPLACE INTO extraction_metadata VALUES (?, ?)",
                        [(k, json.dumps(v)) for k, v in meta.items()])
    log.info("Indexed topology columns and stored %d provenance keys", len(meta))
    log.info("Provenance: threshold=%d cells, router=%s, DEM sha256=%s…",
             threshold_cells, router, dem_hash[:12])
    return meta


def _sha256(path: str, chunk: int = 1 << 20) -> str:
    h = hashlib.sha256()
    with open(path, "rb") as fh:
        for block in iter(lambda: fh.read(chunk), b""):
            h.update(block)
    return h.hexdigest()


# --- Example usage ---
# write_stream_geopackage(
#     links=net, out_path="basin_streams.gpkg",
#     dem_path="basin_conditioned.tif", threshold_cells=500,
# )

Parameter Reference

Choice Recommended Why
Driver GPKG Single file, multi-layer, SQLite-queryable topology
Engine pyogrio Several times faster than Fiona for large writes
Node dtype int64 Float node ids break equality matching after a round trip
Simplification None, in the topological layer Moving a vertex breaks the shared-endpoint invariant
Indexes from_node, to_node Turns an upstream query from a table scan into a lookup
Metadata Separate table in the same file Travels with the data instead of in a README nobody reads

Worked Example: Verifying the Round Trip

The point of the schema is that the file can be reloaded and traversed without repair. The check is short and worth running on every export.

Check Expected A failure means
Link count after reload Unchanged The driver dropped rows — usually a mixed geometry type
Node id dtype after reload Integer Ids were written as real and will not match
Distinct to_node values with no matching from_node One per basin outlet Extra values are dangling links
Weakly connected components One per basin More means the topology fragmented
Sum of length_m Matches the pre-write total Geometry was altered on write
What the Round-Trip Check Compares Five properties compared before writing and after reloading: link count, node dtype, dangling links, connected components and total length. All five must match, and each mismatch points at a specific cause. property in memory after reload link count 3 778 3 778 node id dtype int64 int64 connected components 1 1 total length (km) 3 776.4 3 776.4 run this on every export — it costs a second and catches a whole class of silent corruption

The indexes are what make the file useful as a database rather than as a picture. An upstream query on an unindexed topology column scans the whole table.

Two Indexes Turn a Scan into a Lookup An upstream query on a 120 000-link network takes 4.1 seconds without indexes on from_node and to_node and 0.03 seconds with them. Creating the indexes takes 0.8 seconds once. no indexes 4.1 s indexed 0.03 s — a 140× speedup one upstream query, 120 000-link network Building both indexes costs 0.8 seconds, once, at write time.

Storing the raster-side identifiers too

A segment’s usefulness often depends on being able to get back to the cells it came from. Two small columns make that possible without storing the whole cell path: the row and column of the segment’s head cell, and of its tail cell. With those, a consumer can re-enter the raster at either end of any segment — to sample slope along it, to re-derive its upstream area from the accumulation grid, or to seed a delineation at its outlet — without re-running the tracing step.

The alternative, matching a segment’s geometry back to cells by coordinate, is both slower and lossy: a line vertex sits at a cell centre only if nothing has moved it, and something usually has.

Layer conventions that make the file self-describing

A GeoPackage holds several layers, and using that deliberately is worth the small extra effort. A streams layer carrying the segments, a nodes layer carrying the junctions as points, and the extraction_metadata table described above together make a file that a consumer can understand without documentation. The nodes layer in particular is cheap — it is the unique endpoints already computed during tracing — and it lets a map client symbolise confluences and outlets without deriving them.

Where the same file will hold several extractions at different thresholds, suffix the layer names with the threshold rather than writing several files. The metadata table then carries one row set per layer, and a comparison between two extractions becomes a join inside one file rather than a join across two.

Gotchas and Edge Cases

  • rasterio.features.shapes for vectorization. Produces bank outlines, not centrelines, and fragments at every junction. It is the fastest route to a network that looks right and cannot be traversed.
  • Shapefile output. Field names truncate to ten characters, the format caps at 2 GB, and the CRS travels in a sidecar file that gets separated. Use GeoPackage.
  • Float node ids. Survive most round trips and then fail an equality join on a handful of links, fragmenting the graph in a way that only shows up in a component count.
  • Simplifying before writing. Moves endpoints, breaks shared-node identity, and the damage is invisible until a traversal fails.
  • Mixed geometry types in one layer. A traced network that accidentally includes a zero-length link produces a Point among the LineStrings, and some drivers drop the whole row rather than the geometry.
  • Metadata in the filename. streams_500cells_d8.gpkg is better than nothing and worse than a metadata table, because the filename does not survive being copied into a project directory with a different convention.