⚙️ How-To Hydrology & watersheds Terrain & elevation Raster & rasterio Hydrology & Watersheds · step 4 of 20 Updated

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,406 × 3,470 cells:

  • pysheds computed D8 accumulation in 0.28 s; WhiteboxTools in 0.63 s. Both reached a maximum of 5,134,838 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.58% drained more than 1 km². The mean was 202 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()))
5134838.0 513.4838

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

Bar chart of flow accumulation percentiles inside the Esopus basin on a logarithmic scale: median 10 cells, 90th percentile 66, 99th 2,590 and 99.9th 352,927.
Six orders of magnitude between a typical cell and the river: accumulation needs a log 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 11.95 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.09 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.28 s with the default iterative algorithm; in an earlier run, algorithm="recursive" gave identical values in 0.39 s against 0.35 s. 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,134,838 cells is 513.48 km². WhiteboxTools can do the conversion itself: out_type="catchment area" returned square metres, a maximum of 513,483,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.58% drained at least 1 km². The median was 10 cells, the mean 2,024. 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,924,422 cells, 492.44 km².

Two panels comparing a linear colour ramp of flow accumulation, where almost every cell is background, with a logarithmic ramp that shows tributaries.
The data are the same; the linear ramp spends its whole colour range on a few hundred river cells.

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,924,422 cells = 492.44 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.48 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           10 cells = 0.0010 km2
p90           66 cells = 0.0066 km2
p99        2,590 cells = 0.2590 km2
p99.9    352,927 cells = 35.2927 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 returned a maximum of 5,134,839 cells from the DEM and pysheds 5,134,838 from the same surface. The tools make their own flow-direction choices where neighbouring drops tie and at the edge of the data, so small cell-level differences are expected. On drainage areas of hundreds of square kilometres a one-cell difference is irrelevant; when comparing tools cell by cell, compare flow directions first.

Bar chart of the share of basin cells with an accumulation of exactly 1, at most 10, at most 100, and at least 1 km2.
Nine cells in ten drain less than a hectare; the river is a fraction of a per cent of the grid.

Edge cases or notes

  • NoData handling matters: with the DEM's nodata tag removed, the −999999 border cells were treated as terrain and the gauge's area fell to 448.4 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.

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,406 × 3,470 grid pysheds took 0.28 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,134,838 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?

Almost. On the Esopus DEM the maxima were 5,134,838 and 5,134,839 cells; small differences come from each tool's own flow-direction choices.