Orderings, layerings and partitions for D8 flow networks, in Python.
A D8 flow-direction grid fixes the order in which cells must be visited: a cell can be processed only after its upstream neighbours, or, for some kernels, only after its downstream one. FlowTopo computes that order once, from the flow-direction grid alone, and stores it as a reusable structure. Any routine that walks the network in drainage order (drainage area, flow length, stream order, routing) can then reuse the structure instead of rebuilding the traversal each time.
Three kinds of structure are provided:
- serial orderings (three): a topological sort of the cells, walked in one pass on one core;
- parallel layerings (three): cells grouped into layers so that cells in the same layer are independent of each other and can run on several threads;
- spatial partitions (two): the network cut into independent subregions, one per processor.
Four kernels are bundled to exercise the structures: upstream drainage area, distance to outlet, longest upstream path and Strahler stream order. They are test cases, not the purpose of the package. Each kernel runs on every structure and the results are cross-checked.
In the snippets below, topo is a FlowTopo object built once from a
flow-direction raster; Quick start shows how.
This repository holds the Python reference implementation and one worked example basin. The global 90 m products computed with the same method are released separately on Zenodo; see Global products. The companion paper (Jiang et al.) is in preparation and has no link yet.
A narrated video walks through the structures and the three ways of moving a
value from a cell to its receiver: https://youtu.be/tE5K2wM3TTY. The
animations in this README are reduced previews. Each has a full-size MP4
link beneath it, and all files are in docs/media.
An ordering is a topological sort of the valid cells. One pass over it computes any kernel, because every cell is visited after whatever it needs. Which end that is depends on the kernel: drainage area needs a cell's donors first, flow length needs its receiver first. A sequence is stored in one direction and reversed on demand, and the kernels do the reversing.
topological sort from the sourcesordering="topo" | breadth-first from the pitordering="bfs" | depth-first from the pitordering="dfs" |
|---|---|---|
 |  |  |
| a cell is appended once all its donors are done | cells in order of hop count from the pit | one tributary subtree at a time |
| full-size MP4 | full-size MP4 | full-size MP4 |
The three differ in memory access pattern. Depth-first has the lowest simulated L1 miss rate on the example basin: 10.55%, against 21.83% (breadth-first) and 37.02% (topological sort). Those are one alignment of the two arrays the traversal reads, the one where they collide in every cache set. Shift them apart and the numbers fall to about 8.8%, 14.2% and 36.9%. The order does not move, at any alignment tried.
Direction. An ordering is built in one direction and reversed on demand.
Depth-first and breadth-first start at the pit, so they come out downstream to
upstream (d2u, position 0 is a pit); the topological sort starts at the
headwaters, so it comes out upstream to downstream (u2d, position 0 is a
headwater). The two are reverses of each other.
Which one a kernel needs follows from the way its values travel. Drainage area,
longest upstream path and Strahler order accumulate into the receiver and
need u2d; distance to outlet reads from the receiver and needs d2u. The
kernels flip the sequence for you, so topo.upstream_area(ordering="dfs")
walks the depth-first order upstream to downstream even though it was built the
other way. Ask for a direction explicitly with
topo.ordering("dfs", "u2d").
A layering assigns every cell a layer index such that cells in the same layer are independent of each other. Layers run in order; cells within a layer run in parallel. Layer 0 holds the headwaters.
as soon as possiblelayering="asap" | conflict-free downstreamlayering="cfds" | as late as possiblelayering="alap" |
|---|---|---|
 |  |  |
| every cell in the earliest layer its donors allow | as soon as possible, plus one rule: no two cells in a layer share a receiver | every cell in the latest layer possible |
| full-size MP4 | full-size MP4 | full-size MP4 |
The minimum layer count is set by the longest flow path. The conflict-free rule may add a few layers; on the example basin it adds none (949 layers for all three).
Layerings are built u2d, layer 0 at the headwaters. The downstream-propagating
kernel needs them the other way round, and flips them itself;
topo.decomposition("cfds", "d2u") gives that view directly.
A layering spreads work across threads that share memory. Splitting the network across processors needs a second cut, along the drainage hierarchy, so that no value crosses a subregion boundary while a kernel runs.
Set n_parts to the number of processors: one subregion each, so a subregion's
working set stays in its own memory. The paper's benchmark uses four, because
the server has four Xeon Platinum 8270 processors, each a NUMA node with 26
cores, and runs about 13 threads inside each subregion.
Two basins mapped onto two subregions.level="basin" keeps each basin
whole, so the 62-cell basin and the 11-cell basin cannot be balanced.
level="subbasin" walks up the dominant basin's mainstem, and a tributary
subtree moves to the lighter subregion.
Whole basins cannot be split, so one large basin leaves the other processors idle. The bundled example is one basin, so it shows this plainly:
topo.partition(n_parts=4, level="basin")[1] # [93432, 0, 0, 0]topo.partition(n_parts=4, level="subbasin")[1] # [23121, 23121, 23121, 23120]Subbasin-level walks upstream from the outlet, taking the larger tributary at
each confluence. That isolates the mainstem; the tributary subtrees hanging off
it are dealt to the lighter subregions, and the mainstem runs in a second stage
once they finish. Its cells are marked flowtopo.MAINSTEM. That is why the
four subregion loads above come to 92,483 rather than 93,432: the 949 mainstem
cells are held back for the second stage and do not count as first-stage work.
Whichever structure carries the traversal, a kernel still has to move a value from a cell to its receiver. There are three ways to do it:
| pull | atomic push | push |
|---|---|---|
 |  |  |
| each receiver reads its donors; needs the upstream table | donors write through atomics; correct, but float sums are not reproducible | donors write directly; deterministic, no locks; requires the conflict-free layering |
| full-size MP4 | full-size MP4 | full-size MP4 |
A push is only safe if no two cells in a layer write to the same receiver. Conflict counts on the example basin (93,432 cells):
| layering | conflicting writes inside a layer |
|---|---|
| as soon as possible | 12,122 |
| conflict-free downstream | 0 |
| as late as possible | 39,130 |
The count is a property of the layering and can be checked before running. A test run is not a reliable check: a race does not always trigger. Strahler order has no atomic form, because its confluence rule is a comparison rather than an addition, so its only parallel push is under the conflict-free layering.
Because no two cells in a layer share a receiver, the sums always happen in the
same order: the conflict-free push returns bit-identical results at any thread
count, which an atomic push cannot promise. If manner is not given, FlowTopo
picks a safe one for the layering.
From the paper's benchmark on the full 90 m network (22.2 billion cells, 65 regions):
- One sweep, one core — the depth-first ordering. Up to 5.1× faster than the slowest serial ordering; subtree contiguity drops the L3 miss rate from 37% to about 6%. It is also the baseline to measure parallel speedup against; a slower baseline overstates the speedup.
- Repeated traversal (calibration, ensembles) — the as-late-as-possible layering with pull. Fastest in parallel: each layer's working set stays in cache. Pull must store the donor table.
- Non-linear kernels, or when RAM is tight — the conflict-free downstream layering with push. Lock-free and deterministic, stores only the receiver pointer, and the only parallel option for Strahler order.
- Across processors — the subbasin partition, run at about 13 threads each, beyond which memory bandwidth rather than the algorithm bounds the speedup.
The structures are computed once from the static D8 field and reused without limit.
Those rankings come from the paper's C run at continental scale. This package
is numba over numpy at a much smaller size, and it does not reproduce them
term for term. On a 16-million-cell grid here, push under cfds is the
fastest parallel form, and pull loses even to the serial sequence: building
and reading the donor table costs more than ten threads save. The structures
are the same either way. Which one wins depends on your machine and your grid,
so measure with benchmark.py before choosing.
pip install -e .# numpy + rasterio
pip install -e ".[speed]"# + numba, for threaded kernelsNeeds Python 3.10 or newer. The tests run on 3.10, 3.11 and 3.12. Without numba everything still runs, in pure Python.
importflowtopotopo=flowtopo.FlowTopo.from_raster("data/dir_example.tif")
upa=topo.upstream_area(ordering="dfs") # serialupa=topo.upstream_area(layering="cfds", manner="push") # a layer at a timeldn=topo.distance_to_outlet(ordering="dfs")
lup=topo.longest_upstream_path(ldn, ordering="dfs")
strord=topo.strahler_order(ordering="dfs",
channel_mask=topo.channel_mask(upa, 10.0))
part, load=topo.partition(n_parts=4, level="subbasin") # across processorstopo.to_2d(upa) # back on the gridWith numba, the same kernels with threads:
fromflowtopoimportparallelupa=parallel.upstream_area(topo, layering="cfds", manner="push")Run everything on the example basin:
python example.py # every structure, kernel and manner; ends PASS or FAIL
python benchmark.py # serial vs threaded at several grid sizesdocs/user-guide.md— choosing an ordering, a layering and a manner; threads; raster I/O.docs/methods.md— each method with its origin and its complexity.examples/quickstart.ipynb— a notebook on the bundled data, stored with its output so it reads without running anything.docs/review-checklist.md— what has been checked and how, the bugs those checks found, and the angles still unattacked.
data/dir_example.tif is the example basin of the paper: 292 rows by 614
columns at 3 arc-seconds, 93,432 valid cells, 731 km², cut from
MERIT Hydro (Yamazaki et al., 2019).
The GeoJSON files are the basin boundary and the outlet.
Any D8 GeoTIFF in the same convention works: codes are powers of two clockwise from east, with 0 and 255 terminal. The file's own nodata value is honoured too, which matters because 255 means a terminal here, not nodata.
python example.py --data my_dir.tifThe structures for the whole 90 m network are on Zenodo, one tile per hydrological region, ready to read. To compute anything on them you also need the flow-direction grid they are indexed against, and that is not redistributed here: get it from its authors at https://global-hydrodynamics.github.io/MERIT_Hydro/, under the CC BY-NC 4.0 terms they set.
Working out which of their files a region needs is the fiddly part, and https://fullhydro.org/fullbasin/regions/ does that for you. Pick a region on the map or in the table and it lists the MERIT Hydro tiles covering it, with the row and column offset to place each one. It does not host MERIT Hydro, only the list. The same page has all 96 region outlines as one GeoJSON, tile names and offsets included, for doing it in a script. Those 96 are 65 continental regions, two that straddle the antimeridian and 29 island groups; the products below cover the 65. Region boxes sit on whole degrees and one degree is exactly 1,200 cells at 3 arc-seconds, so a region cuts out of the global rasters by integer arithmetic, with no resampling.
You do not need a whole region. Clip whatever you like out of one and the
structures clip with it: filter a released sequence to the cells you kept and
it is still a topological sort of them, and a filtered layering keeps its
layers mutually independent, the conflict-free guarantee included. That holds
for any subset, because dropping cells from a valid order cannot put a cell
before something it depends on, and dropping cells from a layer cannot make
two of the survivors depend on each other. Checked on a 12,809-cell subbasin
of the bundled example: all three orderings stayed valid, all three layerings
kept their layers independent, cfds kept zero conflicts, and recomputing
from scratch on the clip reproduced all 12,809 values bit for bit.
Clip however suits you: a basin, a rectangle, a country. Two things behave differently and are worth keeping apart.
The structures describe the network you hand over. Clip it and the
structures of the clip are exact: still a topological sort, still layers of
independent cells, still zero conflicts under cfds. Any clip, no exceptions.
A kernel answers a question about that same network. Drainage area asks how much area drains into a cell, and how much area is however much you supplied. Clip a catchment in half and the answer halves, not because the computation slipped but because you asked about a smaller catchment. The number is exactly right for the data it was given.
So the only question is whether your clip contains the whole catchment of the cells you care about. A basin does, by definition: clipping one out of the bundled example and recomputing reproduces every cell to the last decimal. A rectangle does for most of them. The cells it gets wrong are exactly the ones with catchment outside the cut, so the share depends on how much drainage the cut intercepts rather than on any fixed number: across thirteen rectangles through the bundled example, between 92.7% and 99.6% of the kept cells came back identical, median 97.8%. For the rest, the ones below a channel the cut crossed, the number describes your rectangle rather than the world.
Clip a rectangle when the area you kept is what you are studying. Clip whole basins, or read the values out of MERIT-DrainAttr, when the numbers have to mean what they mean globally.
Any D8 grid in this convention works. Build the eight structures yourself and compute on them:
importflowtopotopo=flowtopo.FlowTopo.from_raster("my_dir.tif")
seq=topo.ordering("dfs") # one of three orderingslayers, n=topo.layering("cfds") # one of three layeringspart, load=topo.partition(n_parts=4) # one of two partitionsupa=topo.upstream_area(ordering="dfs") # or any kernel, on any structureThis is the same code that produced the released structures, so a region you build yourself and a region you download are the same thing.
Work a region or a basin at a time rather than the globe. A region encloses only complete basins and stays within 38° × 38°, which keeps its cell indices inside 32-bit integers, which is why this package uses int32 throughout. The region boundaries are in MERIT-FullBasin, below.
This repository holds the method and one example basin. Applied to the whole 90 m MERIT Hydro network, it produces two Zenodo records, both as per-region GeoTIFFs for 65 regions.
MERIT-FlowTopo (10.5281/zenodo.20653058) holds the traversal structures themselves. Every region carries the depth-first sequence, the conflict-free downstream and as-late-as-possible layerings, and the subbasin partition; Region 43 (South China) carries all eight, so the alternatives can be compared somewhere. 978 GB uncompressed, 49 GB compressed.
All eight structures, drawn over the whole network: three serial orderings, three parallel layerings, two spatial partitions. The release itself carries four of them per region, and all eight for Region 43.
Because the D8 field does not change, these are computed once and reused without limit. That is the point: a cost every tool currently pays on every run becomes a read.
MERIT-DrainAttr (10.5281/zenodo.20686664) is the four kernels run over those structures: a baseline, so the variables almost everyone wants need not be recomputed either.
Flow length downstream, flow length upstream and Strahler stream order for every region; upstream drainage area for Region 43 only, since MERIT Hydro already distributes it globally. 622 GB uncompressed, 60 GB compressed.
Anything else is one pass over a structure you already have: a different accumulation, a routing state, a variable nobody has asked for yet. That is why the release is the structures and not the variables.
MERIT-FullBasin (10.5281/zenodo.20344112) divides the network into the 65 hydrologically independent regions everything above is organised by. It comes from the companion dataset, not from here.
MERIT Hydro (Yamazaki et al., 2019; 10.1029/2019WR024873) provides the flow-direction field everything here traverses: the example basin is cut from it, and the global products are built on its 90 m network. The bundled excerpt keeps its CC BY-NC 4.0 terms.
The representation this package works on comes from pyflwdir (D. Eilander, Deltares and the Institute for Environmental Studies, Vrije Universiteit Amsterdam; MIT licence; 10.5281/zenodo.4287337). What FlowTopo takes from it:
- the flat downstream-pointer array — for every cell, the linear index of the cell it drains into, with a pit pointing at itself. Every structure and every kernel here is derived from that one array;
- the D8 decoding conventions — codes as powers of two clockwise from east, 0 and 255 terminal, 247 nodata, and a cell draining off the grid or into nodata treated as a pit;
- the donor-count array and its sentinel convention;
- the chain-tracing rank and the breadth-first sequence builder, which
appear here in the variants set out in
docs/methods.md.
No pyflwdir source is included. The code here was written against those
conventions rather than copied from them, and FlowTopo.from_d8 takes the same
array pyflwdir.from_array(d8, ftype="d8") takes, so the two read the same
rasters.
The three layerings, the conflict-free downstream rule, the propagation manners, the locality metrics and the benchmark drivers are this project's own.
Claude (Anthropic) helped prepare this repository: the Python implementation, the tests, the documentation and the packaging were drafted with it and checked by the authors, and the commit history records where. The structures, the algorithms and the results they produce are the authors' own work, described in the companion manuscript.
Expected outputs for the example basin are stored with the tests, so a fresh clone verifies itself:
pytest170 tests. Every ordering, layering, partition, kernel and manner is cross-checked on the example basin, write-conflict counts included. The structures are checked against their definitions directly: an ordering is a topological sort, a layer is an antichain, a receiver comes after its donors, and all of that survives clipping a basin out of a region. Degenerate inputs get their own: an empty grid, a lone cell, networks with cycles. So do the API's promises: cached arrays are read-only, repeated calls agree bit for bit, accumulation keeps the precision it was given.
Those tests and example.py run on Python 3.10, 3.11 and 3.12 on every push.
The badge at the top links to the runs.
Questions, problems and suggestions are welcome by email at lulu_jiang@pku.edu.cn, or as a GitHub issue. Email is the surer way to reach us.
The companion manuscript is MERIT-FlowTopo v1.0: a reusable computational foundation for hyperresolution hydrology on the global 90 m drainage network (Jiang et al., in preparation). A citation file will be added once it appears; until then, cite the manuscript.
MIT for the code; see LICENSE. The bundled MERIT Hydro excerpt keeps
its own CC BY-NC 4.0 terms; see DATA_NOTICE.md.