Hierarchical Pfafstetter Basin Coding in Python

A Pfafstetter code is a single integer string that encodes both the identity and the drainage topology of a catchment: read the digits left to right and you walk down the containment hierarchy, compare two codes and you know instantly which basin lies upstream. This guide implements that numbering scheme from a stream network in Python, as a focused reference within the Basin Partitioning Strategies topic, which itself sits under the broader Watershed Delineation & Catchment Synchronization coverage on this site.

The scheme was designed by Otto Pfafstetter to be topologically self-consistent, storage-efficient, and recursively subdividable to arbitrary depth. The global HydroSHEDS project popularised it through the HydroBASINS dataset, whose PFAF_ID field carries exactly the codes this page reconstructs. Understanding the assignment rules lets you extend those codes into custom, high-resolution catchments that the global product does not resolve.

The numbering rules in one pass

At every level, a basin is divided into nine units by examining the tributaries that join its main stem. The four tributaries contributing the largest drainage area receive even codes; the segments of the main stem valley between and around them receive odd codes.

The odd codes number the interbasins: interbasin 1 is the drainage downstream of tributary 2, interbasin 3 lies between tributaries 2 and 4, and so on up to interbasin 9, which is the headwater area above tributary 8. Because both series are assigned from the outlet upstream, a larger digit always denotes a unit further up the drainage.

The rule is scale-independent, which is the whole point. Once basin 8 exists, you subdivide it by finding its four largest tributaries and appending a second digit — 81 through 89 — and the string 84 now means “tributary basin 4 inside interbasin-headwater 8 of the parent.” That prefix property is what makes the code so useful for spatial joins and hierarchical rollups. The same recursive subdivision motivation drives the DEM-based approach in partitioning large watersheds into sub-basins with RichDEM; Pfafstetter differs in that it derives its hierarchy from vector stream topology rather than from a flow-accumulation raster.

Prerequisites

Before coding you need a stream network with resolved connectivity:

• A LineString GeoDataFrame of stream reaches in a projected CRS with metres, so lengths and areas are meaningful. If your reaches arrive in mixed projections, resolve them first — see coordinate reference system alignment.
• A downstream pointer per reach (a next_down id), which you can derive from a flow-direction grid or from snapping reach endpoints.
• A per-reach contributing drainage area, ideally carried over from flow accumulation extraction so that tributaries can be ranked by the area they deliver at their confluence.

conda create -n pfaf python=3.11 geopandas shapely numpy networkx -c conda-forge
conda activate pfaf

Building the confluence hierarchy

The core algorithm walks the main stem from outlet to source, records each confluence, and ranks the entering tributaries. The networkx directed graph makes “which reaches drain into this node” a one-line query.

import logging
import geopandas as gpd
import networkx as nx

logger = logging.getLogger(__name__)


def build_stream_graph(reaches: gpd.GeoDataFrame) -> nx.DiGraph:
    """Construct a downstream-directed graph from a reach GeoDataFrame.

    Each reach must carry a unique 'reach_id', a 'next_down' id (or -1 at the
    outlet), and an 'up_area' contributing drainage area in square kilometres.
    """
    graph = nx.DiGraph()
    for _, reach in reaches.iterrows():
        graph.add_node(reach.reach_id, up_area=reach.up_area)
    for _, reach in reaches.iterrows():
        if reach.next_down != -1:
            graph.add_edge(reach.reach_id, reach.next_down)
    logger.info(
        "Built stream graph: %d reaches, %d downstream edges",
        graph.number_of_nodes(), graph.number_of_edges()
    )
    return graph

Assigning codes recursively

The heart of the implementation ranks tributaries by contributing area, hands the top four the even codes, and fills the interbasins with the odd codes. It then recurses into every non-trivial unit, appending a digit each time.

import logging
import geopandas as gpd

logger = logging.getLogger(__name__)


def assign_pfafstetter(
    reaches: gpd.GeoDataFrame,
    main_stem_ids: list[int],
    tributaries: list[dict],
    prefix: str = "",
    max_level: int = 3,
) -> dict[int, str]:
    """Assign Pfafstetter codes to reaches within one basin.

    Parameters
    ----------
    reaches       : full reach GeoDataFrame indexed by reach_id
    main_stem_ids : reach ids of the main stem, ordered outlet -> source
    tributaries   : list of {'confluence_pos': int, 'up_area': float,
                    'member_ids': list[int]} for each tributary entering the stem
    prefix        : accumulated code from parent levels
    max_level     : recursion depth guard

    Returns
    -------
    mapping of reach_id -> Pfafstetter code string
    """
    codes: dict[int, str] = {}

    # Rank tributaries by the drainage area they contribute at their confluence
    ranked = sorted(tributaries, key=lambda t: t["up_area"], reverse=True)
    top_four = sorted(ranked[:4], key=lambda t: t["confluence_pos"])
    logger.info(
        "Prefix '%s': %d tributaries, selecting %d largest for even codes",
        prefix or "<root>", len(tributaries), len(top_four)
    )

    # Even codes 2,4,6,8 to the four largest tributaries, downstream -> upstream
    even_positions = []
    for even_code, trib in zip((2, 4, 6, 8), top_four):
        for member in trib["member_ids"]:
            codes[member] = f"{prefix}{even_code}"
        even_positions.append(trib["confluence_pos"])
        logger.debug("Tributary at pos %d -> even code %d", trib["confluence_pos"], even_code)

    # Odd codes 1,3,5,7,9 to the interbasins delimited by those confluences
    boundaries = [0] + sorted(even_positions) + [len(main_stem_ids)]
    for odd_index, odd_code in enumerate((1, 3, 5, 7, 9)):
        lo, hi = boundaries[odd_index], boundaries[odd_index + 1]
        for stem_id in main_stem_ids[lo:hi]:
            codes[stem_id] = f"{prefix}{odd_code}"

    # Recurse: every unit becomes a new basin one level deeper
    if len(prefix) + 1 < max_level:
        for unit_code in ("1", "2", "3", "4", "5", "6", "7", "8", "9"):
            child_prefix = f"{prefix}{unit_code}"
            sub = _extract_unit(reaches, codes, child_prefix)
            if sub is not None:
                codes.update(assign_pfafstetter(reaches, *sub, prefix=child_prefix, max_level=max_level))

    return codes


def _extract_unit(reaches, codes, child_prefix):
    """Rebuild the (main_stem_ids, tributaries) inputs for one sub-unit.

    Returns None when the unit holds no further confluences worth subdividing.
    Implementation depends on your reach schema; omitted here for brevity.
    """
    member_ids = [rid for rid, c in codes.items() if c == child_prefix]
    if len(member_ids) < 5:  # too small to hold four ranked tributaries
        logger.debug("Unit %s has %d reaches; leaving as leaf", child_prefix, len(member_ids))
        return None
    # ... derive ordered main stem and tributary list for member_ids ...
    return None


if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO, format="%(levelname)s %(message)s")
    gdf = gpd.read_file("stream_reaches.gpkg")
    logger.info("Loaded %d reaches from stream_reaches.gpkg", len(gdf))
    # main_stem_ids and tributaries built from build_stream_graph(gdf)

The recursion terminates naturally: a unit with fewer than five reaches cannot host four ranked tributaries plus an interbasin, so it stays a leaf. In practice you cap max_level at the digit depth your application needs — HydroBASINS ships levels 1 through 12.

Reading the codes back

The prefix property turns hierarchy queries into string operations. To find every catchment inside basin 84, filter for codes that startswith("84"). To find the immediate downstream neighbour, decrement the last digit.

Because upstream/downstream is decidable from the digits alone, routing a pollutant or aggregating runoff downstream needs no repeated graph traversal once the codes exist. That is the practical payoff of the scheme over an arbitrary integer id.

Validation checks

1. Prefix containment: for every reach, assert its code begins with its parent unit’s code. A violation means a tributary was assigned to the wrong interbasin during ranking.
2. Complete partition: the union of all nine unit geometries at any level must equal the parent geometry with no gaps or overlaps. Verify this the same way you would any coverage, using the techniques in boundary topology validation.
3. Even-code count: each level should assign at most four even codes. More than four means a tie in contributing area was resolved inconsistently — break ties deterministically by reach id.

Gotchas & edge cases

• Braided or distributary channels break the single-main-stem assumption. Pfafstetter presumes one dominant downstream path per confluence. In deltas and anabranching rivers you must first designate a principal channel, or the tributary ranking becomes ambiguous.
• Ties in contributing area produce nondeterministic codes. When two tributaries deliver nearly identical area, floating-point ranking can flip between runs. Sort with a secondary key such as reach id so the output is reproducible.
• Fewer than four tributaries leaves even codes unused. A basin with only two significant tributaries assigns codes 2 and 4, and interbasins 1, 3, 5 absorb the rest. This is valid; do not pad with phantom units.
• Contributing area, not length, drives ranking. Selecting the four longest tributaries instead of the four largest by area produces a numbering that looks plausible but violates the standard, breaking joins against HydroBASINS.

Related Topics

• Basin Partitioning Strategies — parent topic covering how and why large watersheds are split into modelling units
• Partitioning Large Watersheds into Sub-Basins with RichDEM — the raster-based counterpart that derives sub-basins from flow accumulation rather than vector topology
• Boundary Topology Validation — how to confirm that a coded partition covers its parent with no gaps or overlaps
• Watershed Delineation & Catchment Synchronization — the overview tying together delineation, coding, and reconciliation workflows