Skip to contents

Network properties: accumulation, ordering, and coding

Source: vignettes/network-properties.Rmd
network-properties.Rmd

A hydrofabric flow network is a directed acyclic graph (DAG): every feature (id) points to exactly one downstream feature (toid), and terminal features point to 0 (or NA). hfutils provides a family of topology functions that operate on this structure, all are character-safe, so identifiers such as "fp-123" or scientific-notation strings round-trip cleanly, and each returns a vector aligned to the input rows.

Downstream accumulation

accumulate_downstream() propagates a per-feature attribute (drainage area, incremental length, anything additive) downstream, summing all upstream contributions at confluences. It runs a single topological sort followed by one O(E) edge pass, so it scales to continental networks.

Consider two headwaters (1, 2) joining at 3, which flows to the outlet 4:

df <- data.frame(
  flowpath_id   = c("1", "2", "3", "4"),
  flowpath_toid = c("3", "3", "4", "0"),
  area          = c(1.0, 2.0, 0.5, 0.0)
)

accumulate_downstream(df, attr = "area")
#> [1] 1.0 2.0 3.5 3.5

Feature 3 accumulates 1.0 + 2.0 + 0.5 = 3.5, and the outlet inherits the full basin total. The column names default to flowpath_id / flowpath_toid but can be overridden via the id and toid arguments to match any network table (including reference-fabric id/toid).

A cyclic network is rejected rather than silently producing wrong totals:

cyclic <- data.frame(
  flowpath_id   = c("1", "2"),
  flowpath_toid = c("2", "1"),
  area          = c(1, 1)
)

accumulate_downstream(cyclic, attr = "area")
#> Error in `accumulate_downstream()`:
#> ! Network contains cycles; cannot accumulate.

Hydrosequence

get_hydroseq() assigns a total topological ordering to the network. Larger values are upstream; the value decreases monotonically as you move toward the outlet, which makes it a convenient sort key for downstream traversal.

df$hydroseq <- get_hydroseq(df)
df[order(-df$hydroseq), c("flowpath_id", "flowpath_toid", "hydroseq")]
#>   flowpath_id flowpath_toid hydroseq
#> 4           4             0        4
#> 2           2             3        3
#> 1           1             3        2
#> 3           3             4        1

Derived network attributes

The rest of the family builds the standard NHDPlus-style routing attributes on top of the same topological pass. Take a small basin, a mainstem 4 -> 3 -> 1 with a tributary 2 -> 1 draining to the outlet 1:

net <- data.frame(
  flowpath_id   = c("1", "2", "3", "4"),
  flowpath_toid = c("0", "1", "1", "3"),
  lengthkm      = c(4, 2, 3, 5)
)

net$arbolate  <- accumulate_downstream(net, attr = "lengthkm")   # arbolate sum
net$order     <- get_streamorder(net)                            # Strahler order
net$levelpath <- get_levelpath(net, weight = "arbolate")         # mainstem grouping
net$pathlen   <- get_pathlength(net, length = "lengthkm")        # distance to outlet
net$level     <- get_streamlevel(net, levelpath = "levelpath")   # stream level
net
#>   flowpath_id flowpath_toid lengthkm arbolate order levelpath pathlen level
#> 1           1             0        4       14     2         4       0     1
#> 2           2             1        2        2     1         1       4     2
#> 3           3             1        3        8     1         4       4     1
#> 4           4             3        5        5     1         4       7     1
  • get_streamorder(), Strahler order: leaves are 1, and order increases only where two equal-order streams meet (so the confluence at 1 is order 2).
  • get_levelpath(), groups reaches into continuous mainstems, following the largest-weight contributor at each confluence (here the arbolate sum, so 1, 3, 4 share a level path and the tributary 2 is its own).
  • get_pathlength(), distance along the network from each reach’s outlet to the network terminus (0 at the outlet, summing downstream lengths going up).
  • get_streamlevel(), how many level-path steps a reach is from the terminus: the mainstem is 1, the tributary 2.

Pfafstetter codes

get_pfafstetter() assigns hierarchical Pfafstetter basin codes. It needs total drainage area, a hydrosequence, and level paths precomputed, then subdivides each basin’s mainstem by its four largest tributaries, recursing max_level digits deep:

# on a real fabric this is accumulated catchment area; `arbolate` stands in here
net$total_da_sqkm <- net$arbolate
net$topo_sort     <- get_hydroseq(net)
net$pfaf          <- get_pfafstetter(net, max_level = 2)

The result is most meaningful on a full basin, where each digit narrows the location within the drainage hierarchy. On a four-reach toy network there are no tributaries large enough to subdivide, so every code comes back NA.

Upstream queries and partitioning

upstream_index() assigns each feature a nested set pair, upstream_id (its depth-first pre-order position) and num_upstreams (how many features sit strictly above it). Together these turn “everything upstream of X” into an integer range test rather than a graph traversal:

idx <- upstream_index(net)
cbind(net["flowpath_id"], idx)
#>   flowpath_id upstream_id num_upstreams
#> 1           1           1             3
#> 2           2           4             0
#> 3           3           2             1
#> 4           4           3             0

Everything upstream of a feature with upstream_id == u and num_upstreams == k is exactly the half-open range (u, u + k]:

u <- idx$upstream_id[net$flowpath_id == "1"]
k <- idx$num_upstreams[net$flowpath_id == "1"]
net$flowpath_id[idx$upstream_id > u & idx$upstream_id <= u + k]
#> [1] "2" "3" "4"

Because the walk expands the largest-upstream branch first, a mainstem stays contiguous. Note that upstream_id is build-specific: it changes whenever the topology changes, so it is an index, not a persistent key.

For a hydrofabric whose flowpaths route through nexuses, hf_upstream_index() resolves the flowpath -> nexus -> flowpath hops first and then indexes the resulting flowpath graph. write_hydrofabric() calls it for you, so a written GeoPackage carries the index already.

merge_groups() builds on the same pre-order to cut the network into contiguous runs, breaking wherever a feature does not flow directly into its predecessor or (with order) wherever the given order column changes. Each group is a contiguous upstream_id range, so a size-budgeted partitioner can merge whole groups and always end up with complete sub-networks:

merge_groups(net, order = "levelpath")
#> [1] 1 2 1 1

The mainstem 4 -> 3 -> 1 forms one group; the tributary 2 is its own.

In practice

Every function above accepts a data frame, tibble, or sf object. A typical pattern reads a layer lazily, filters to a VPU, materializes it, and accumulates:

library(dplyr)
library(sf)

da <- as_ogr("conus_nextgen.gpkg", "flowpaths") |>
  filter(vpuid == "01") |>
  st_as_sf() |>
  accumulate_downstream(attr = "areasqkm")

See vignette("reading-and-writing") for the I/O side of that pipeline.