How to Calculate Flow Accumulation from a DEM in Python

Problem statement

Flow accumulation is the grid behind stream networks, watershed outlets and drainage areas: for every cell, the number of cells — or the area — that drains through it. Computing it takes a single function call in pysheds or WhiteboxTools once the DEM is conditioned. Using it well takes knowing its units, its extraordinary range, what a weighted version needs, and why a map of it usually looks empty.

Measured on the conditioned USGS 3DEP 10 m DEM of the Esopus Creek catchment in New York, 2,408 × 3,471 cells:

• pysheds computed D8 accumulation in 0.29 s; WhiteboxTools in 0.58 s. Both reached a maximum of 5,135,399 cells, 513.5 km².
• Half the basin's cells had an accumulation of 10 or less, one in ten had exactly 1, and only 0.57% drained more than 1 km². The mean was 184 times the median.
• 99.75% of cells were under 1% of the maximum, which is why a linear colour ramp shows nothing but the main river.
• Passing a plain NumPy array as weights failed with AttributeError: 'numpy.ndarray' object has no attribute 'nodata'; wrapped as a pysheds Raster, cell-area weights gave exactly cells × 100 m² everywhere.

Quick answer

import numpy as np

if not hasattr(np, "in1d"):          # pysheds 0.5 on NumPy 2.4+
    np.in1d = lambda a, b, **kw: np.isin(np.ravel(a), b, **kw)
from pysheds.grid import Grid

grid = Grid.from_raster("esopus_breached.tif")        # a depression-free DEM
dem = grid.read_raster("esopus_breached.tif")
fdir = grid.flowdir(dem)
acc = grid.accumulation(fdir)                         # cells draining through each cell
area_km2 = acc * abs(grid.affine.a * grid.affine.e) / 1e6
print(float(acc.max()), float(area_km2.max()))

5135399.0 513.5399

Accumulation counts cells, including the cell itself. Multiply by cell area for drainage area, and view it on a logarithmic scale.

Step-by-step solution

1. Start from a conditioned DEM

Accumulation follows flow directions, and flow directions stop at pits. On the raw Esopus DEM the maximum accumulation was 13.03 km²; on the breached DEM it was 513.5 km². Condition first — see breaching or filling a DEM.

2. Compute flow direction

grid.flowdir(dem) returns D8 directions in 0.08 s on this grid. The routing method changes accumulation on hillslopes much more than on channels; see flow direction explained.

3. Accumulate

grid.accumulation(fdir) returned a float64 grid in 0.29 s with the default iterative algorithm; algorithm="recursive" gave identical values in 0.32 s. (Best of three runs on a 12-core Ryzen 9 7900, load average 4.7.) Each value counts the cell itself plus every cell upstream, so the minimum is 1.

4. Convert cells to drainage area

On a projected grid, multiply by the cell area: 10 m × 10 m is 0.0001 km², so 5,135,399 cells is 513.54 km². WhiteboxTools can do the conversion itself: out_type="catchment area" returned square metres, a maximum of 513,539,904. On a latitude–longitude DEM cell area varies by row and a constant factor is wrong — see fixing catchment areas on a geographic DEM.

5. Weight accumulation when cells contribute unequally

Weighted accumulation sums a value per cell instead of counting cells: rainfall, runoff, land-cover share, pollutant load. pysheds expects the weights as a Raster sharing the grid's view; a bare array failed with an AttributeError about nodata. With cell area as the weight, the result matched cells × cell area in every cell, confirming the weighting (Example 2).

6. Understand the distribution before using it

Inside the basin, 10.0% of cells had an accumulation of exactly 1, 50.0% had 10 or less, 93.1% had 100 or less, and 0.57% drained at least 1 km². The median was 11 cells, the mean 2,020. Any statistic that averages accumulation, and any threshold chosen without looking at this distribution, will surprise you.

7. Display it on a logarithmic scale

With a linear colour ramp, 99.75% of cells fall in the bottom 1% of the range and render as background; only the trunk stream is visible. np.log10(acc) spreads the range from 0 to 6.7 and makes tributaries and hillslope flow paths visible (Example 3).

8. Read drainage areas at outlets, not at gauge coordinates

The accumulation at a gauge's exact coordinate is often tiny because the derived channel is a cell or two away. Take the largest accumulation within a short distance, or snap the outlet first — see pour points explained. Within 150 m of the Coldbrook gauge the largest value was 4,923,548 cells, 492.35 km².

Code examples

Example 1 — accumulation and the drainage area at a gauge

from pyproj import Transformer

cell_km2 = abs(grid.affine.a * grid.affine.e) / 1e6
x, y = Transformer.from_crs("EPSG:4269", grid.crs.srs, always_xy=True).transform(-74.2701944, 42.0144722)
col, row = ~grid.affine * (x, y)
row, col = int(row), int(col)

A = np.asarray(acc)
window = A[row - 15:row + 16, col - 15:col + 16]              # 150 m either side on a 10 m grid
print(f"accumulation at the gauge cell: {A[row, col]:,.0f} cells")
print(f"largest within 150 m: {window.max():,.0f} cells = {window.max() * cell_km2:.2f} km2")

accumulation at the gauge cell: 1 cells
largest within 150 m: 4,923,548 cells = 492.35 km2

Example 2 — weighted accumulation

from pysheds.view import Raster

cell_area = np.full(dem.shape, cell_km2)
try:
    grid.accumulation(fdir, weights=cell_area)
except AttributeError as error:
    print("plain array:", error)

weights = Raster(cell_area, viewfinder=grid.viewfinder)
area = np.asarray(grid.accumulation(fdir, weights=weights))
print(f"weighted maximum {area.max():.2f} km2; equals cells x cell area everywhere: "
      f"{np.allclose(area, np.asarray(acc) * cell_km2)}")

plain array: 'numpy.ndarray' object has no attribute 'nodata'
weighted maximum 513.54 km2; equals cells x cell area everywhere: True

Replace the constant with a rainfall or runoff grid on the same cells to accumulate volumes instead of areas.

Example 3 — describe and display the distribution

import geopandas as gpd
import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt
from rasterio import features

basin = gpd.read_file("nldi_basin_01362500.geojson").to_crs(grid.crs.srs)
inside = features.rasterize(basin.geometry, out_shape=dem.shape, transform=grid.affine).astype(bool)
values = A[inside]
for q in (50, 90, 99, 99.9):
    v = np.percentile(values, q)
    print(f"p{q:<4} {v:>10,.0f} cells = {v * cell_km2:.4f} km2")
print(f"share of cells below 1% of the maximum: {(A < 0.01 * A.max()).mean():.2%}")

fig, (left, right) = plt.subplots(1, 2, figsize=(10, 4))
left.imshow(np.where(inside, A, np.nan), cmap="Blues")
left.set_title("linear")
right.imshow(np.where(inside, np.log10(A), np.nan), cmap="Blues")
right.set_title("log10")
fig.savefig("accumulation_linear_vs_log.png", dpi=100)

p50           11 cells = 0.0011 km2
p90           66 cells = 0.0066 km2
p99        2,574 cells = 0.2574 km2
p99.9    352,002 cells = 35.2002 km2
share of cells below 1% of the maximum: 99.75%

Explanation

Why accumulation spans so many orders of magnitude

Drainage networks branch: every cell on a hillslope drains a handful of neighbours, and only the cells along the main stem collect the whole basin. The number of cells draining a given area falls steeply as the area grows, so a grid of accumulation is almost entirely small values with a thin tree of very large ones.

Why the gauge cell has almost no accumulation

The derived channel is exactly one cell wide. A gauge located even 10 m to one side of it sits on the bank, whose accumulation is a few cells. Accumulation is not a smooth field that can be sampled at a point; it has to be read on the channel.

Why pysheds needs weights as a Raster

pysheds keeps grid metadata — affine transform, CRS, nodata — on its Raster objects and uses the nodata value while routing. A bare NumPy array has none of that, so the accumulation routine fails when it looks for the weights' nodata. Wrapping the array with the grid's viewfinder supplies it.

Why two tools can differ slightly

WhiteboxTools and pysheds both returned a maximum of 5,135,399 cells from this DEM. The tools still make their own flow-direction choices where neighbouring drops tie and at the edge of the data, so cell-level differences are possible elsewhere in the grid. On drainage areas of hundreds of square kilometres a one-cell difference is irrelevant; when comparing tools cell by cell, compare flow directions first.

Edge cases or notes

• NoData handling matters: with the NoData tag removed from a reprojected DEM, its −999999 corner cells were treated as terrain and the gauge's area fell from 492.4 to 448.0 km²; see fixing flow accumulation that is all zeros or NoData.
• Voids inside the DEM block accumulation downstream of them.
• Clipped DEMs lose every upstream cell outside the clip.
• Integer overflow is possible if accumulation is cast to 32-bit integers on continental grids; keep it as float64 or int64.
• Specific contributing area divides by contour width; it is not the same as catchment area.
• D∞ and MFD accumulation differ from D8 mainly on hillslopes.
• Stream thresholds are chosen on accumulation; see extracting a stream network.

Internal links

• DEM hydrology explained: from elevation to where water goes — the whole pipeline
• Flow direction explained: D8, D-infinity and MFD compared — the input to accumulation
• How to condition a DEM for hydrology: breaching or filling — before accumulating
• How to extract a stream network from a DEM — thresholds on accumulation
• Pour points explained: why watershed delineation needs a snapped outlet — reading areas at outlets
• Fixing flow accumulation that is all zeros, NoData or stops at the edge — when it goes wrong
• Fixing catchment areas that are wrong on a latitude–longitude DEM — cell area on geographic grids
• How to design a colour ramp in Python — displaying skewed data

FAQ

How do I calculate flow accumulation in Python?

Condition the DEM, compute flow direction with grid.flowdir, then call grid.accumulation(fdir) in pysheds, or run WhiteboxTools' d8_flow_accumulation. On a 2,408 × 3,471 grid pysheds took 0.29 s.

What units does flow accumulation have?

Cells, including the cell itself. Multiply by cell area for drainage area; on a 10 m grid 5,135,399 cells is 513.5 km².

Why is my flow accumulation map almost blank?

Accumulation is extremely skewed: 99.75% of cells were under 1% of the maximum. Display log10(accumulation) instead of the raw values.

How do I compute weighted flow accumulation with pysheds?

Pass the weights as a pysheds Raster with the grid's viewfinder. A plain NumPy array raised an AttributeError about nodata.

Why is the accumulation at my gauge so small?

The gauge is beside the one-cell-wide derived channel. Take the largest accumulation within about 150 m, or snap the outlet to the channel first.

Do pysheds and WhiteboxTools give the same accumulation?

On the Esopus DEM, yes: both gave a maximum of 5,135,399 cells. Elsewhere small differences can come from each tool's own flow-direction choices at ties and edges.