How to Fill Sinks and Derive Flow Direction from a DEM

Problem statement

Every hydrological product — flow direction, accumulation, stream networks, watersheds — assumes water can always run downhill to the edge of the grid. The standard advice is to "fill the sinks" first.

Measure your DEM before following it, because the advice assumes a problem you may not have:

import numpy as np
from scipy.ndimage import label, minimum_filter

neighbour_min = minimum_filter(dem, size=3, mode="nearest")
interior = np.zeros(dem.shape, bool)
interior[1:-1, 1:-1] = True

true_pits = (dem < neighbour_min) & interior          # strictly lower than all 8
no_outlet = (dem <= neighbour_min) & interior         # no strictly-lower neighbour

print(f"true pits:            {true_pits.sum():,} in {label(true_pits)[1]:,} regions")
print(f"cells with no outlet: {no_outlet.sum():,} in {label(no_outlet)[1]:,} regions")

true pits:            0 in 0 regions
cells with no outlet: 1,394 in 131 regions

Zero true pits. Copernicus DEM is void-filled and conditioned, so no cell is strictly lower than all its neighbours. What it does have is 1,394 cells with nowhere to drain — cells on flats, where at least one neighbour is at exactly the same height.

That distinction matters, because filling fixes pits and makes flats worse.

Quick answer

from skimage.morphology import reconstruction


def fill_sinks(dem):
    """Priority-flood: raise every depression to the level of its lowest outlet."""
    seed = np.full_like(dem, dem.max())
    seed[0, :], seed[-1, :] = dem[0, :], dem[-1, :]      # boundary stays put
    seed[:, 0], seed[:, -1] = dem[:, 0], dem[:, -1]
    return reconstruction(seed, dem, method="erosion")


filled = fill_sinks(dem)
raised = filled - dem
print(f"{(raised > 0).sum():,} cells raised, median {np.median(raised[raised > 0]):.2f} m, "
      f"max {raised.max():.2f} m")

after = ((filled <= minimum_filter(filled, size=3, mode="nearest")) & interior)
print(f"cells with no outlet after filling: {after.sum():,} in {label(after)[1]:,} regions")

2,058 cells raised, median 1.07 m, max 10.51 m
cells with no outlet after filling: 2,068 in 115 regions

The count went up, from 1,394 to 2,068. Filling did its job — the 131 depressions became 115 — and every depression it raised became a flat surface at its outlet level, where D8 has no steepest descent.

Pits or flats, priority-flood, D8, accumulation

1. Distinguish pits from flats before doing anything

def drainage_report(surface, label_text):
    neighbour_min = minimum_filter(surface, size=3, mode="nearest")
    interior = np.zeros(surface.shape, bool)
    interior[1:-1, 1:-1] = True

    strict = (surface < neighbour_min) & interior
    stuck = (surface <= neighbour_min) & interior
    print(f"{label_text}")
    print(f"   true pits            {strict.sum():>6,} in {label(strict)[1]:>5,} regions")
    print(f"   no downhill neighbour {stuck.sum():>6,} in {label(stuck)[1]:>5,} regions")


drainage_report(dem, "RAW DEM")

RAW DEM
   true pits                 0 in     0 regions
   no downhill neighbour  1,394 in   131 regions

Two very different diagnoses:

• Many true pits — a noisy DEM. Filling is exactly right and will remove almost all of them.
• Zero true pits, many flats — a conditioned product like Copernicus DEM. Filling changes little and flat-resolution is your real problem.

Running the strict test alone would have reported "no pits, nothing to do", and then produced a fragmented stream network anyway.

2. Fill with priority-flood

seed = np.full_like(dem, dem.max())
seed[0, :], seed[-1, :] = dem[0, :], dem[-1, :]
seed[:, 0], seed[:, -1] = dem[:, 0], dem[:, -1]
filled = reconstruction(seed, dem, method="erosion")

Morphological reconstruction by erosion is priority-flood: it lowers a surface starting at the maximum everywhere except the boundary until it rests on the DEM. The result is the lowest surface above the DEM with no interior minima — each depression raised to exactly its lowest outlet, and no higher.

raised = filled - dem
print(f"{(raised > 0).sum():,} raised, median {np.median(raised[raised > 0]):.2f} m, "
      f"max {raised.max():.2f} m")

2,058 raised, median 1.07 m, max 10.51 m

A median of 1.07 m is consistent with noise in a DEM of about 4 m vertical accuracy. The 10.51 m maximum deserves a look at the hillshade — that is either a real depression or a void-fill artefact.

Filling is destructive. Keep the original and never use the filled surface for slope, contours or elevation reporting.

3. Expect the flats to grow

drainage_report(filled, "FILLED")

FILLED
   true pits                 0 in     0 regions
   no downhill neighbour  2,068 in   115 regions

Regions down from 131 to 115 — depressions merged and resolved. Cells up from 1,394 to 2,068, because every filled depression is now a plateau at its outlet level.

This is not a failure. It is what filling does, and it moves the problem from "water stops here" to "water does not know which way to go here".

4. Compute D8 flow direction

OFFSETS = [(-1, 0), (-1, 1), (0, 1), (1, 1), (1, 0), (1, -1), (0, -1), (-1, -1)]
CODES = [64, 128, 1, 2, 4, 8, 16, 32]      # ESRI convention, clockwise from north


def d8_direction(surface, cell_x, cell_y):
    best_drop = np.full(surface.shape, -np.inf)
    direction = np.zeros(surface.shape, dtype="uint8")

    for (dr, dc), code in zip(OFFSETS, CODES):
        shifted = np.roll(np.roll(surface, -dr, axis=0), -dc, axis=1)
        distance = np.hypot(dr * cell_y, dc * cell_x)       # diagonals are further
        drop = (surface - shifted) / distance
        better = drop > best_drop
        best_drop = np.where(better, drop, best_drop)
        direction = np.where(better, code, direction).astype("uint8")

    direction[best_drop <= 0] = 0                            # flats and edges
    direction[0, :] = direction[-1, :] = 0                   # np.roll wraps — mask it
    direction[:, 0] = direction[:, -1] = 0
    return direction, best_drop

Dividing by the distance is not optional. Compare raw height drops instead and a diagonal neighbour wins whenever its drop is bigger, even though it is √2 times further away:

def diagonal_share(direction):
    assigned = direction[direction > 0]
    return np.isin(assigned, [128, 2, 8, 32]).mean()

by_gradient, _ = d8_direction(filled, 27.9, 30.7)
raw_best = np.full(filled.shape, -np.inf)
raw_dir = np.zeros(filled.shape, dtype="uint8")
for (dr, dc), code in zip(OFFSETS, CODES):               # raw drop: distance ignored
    drop = filled - np.roll(np.roll(filled, -dr, axis=0), -dc, axis=1)
    raw_dir = np.where(drop > raw_best, code, raw_dir).astype("uint8")
    raw_best = np.maximum(raw_best, drop)
raw_dir[raw_best <= 0] = 0
raw_dir[0, :] = raw_dir[-1, :] = raw_dir[:, 0] = raw_dir[:, -1] = 0
changed = (raw_dir != by_gradient) & (by_gradient > 0)
print(f"diagonal directions: {diagonal_share(by_gradient):.0%} by gradient, {diagonal_share(raw_dir):.0%} by raw drop")
print(f"cells whose direction changes: {changed.mean():.0%}")

diagonal directions: 49% by gradient, 90% by raw drop
cells whose direction changes: 39%

The raw-drop version sends far more of the surface down diagonals, and over whole hillsides that is the distinctive herringbone flow pattern.

5. Accumulate from the top down

def d8_accumulation(surface, direction):
    rows, cols = surface.shape
    order = np.argsort(surface.ravel())[::-1]        # highest cells first
    accumulation = np.ones(surface.size)
    lookup = dict(zip(CODES, OFFSETS))
    flat_dir = direction.ravel()

    for idx in order:
        code = flat_dir[idx]
        if code == 0:
            continue
        dr, dc = lookup[code]
        r, c = divmod(int(idx), cols)
        nr, nc = r + dr, c + dc
        if 0 <= nr < rows and 0 <= nc < cols:
            accumulation[nr * cols + nc] += accumulation[idx]

    return accumulation.reshape(surface.shape)

Processing highest-first guarantees every upstream contribution has arrived before a cell is passed on. It works only because D8 gives each cell exactly one downstream neighbour, so the flow graph is a forest.

Fill, flow, streams and a watershed on Snowdon

The data used on this page

The Copernicus DEM (GLO-30) over the Snowdon massif, read from its public bucket. The window crosses the 4°W tile boundary, so two tiles are merged.

import numpy as np
import rasterio
from rasterio.merge import merge

BBOX = (-4.10, 53.03, -4.00, 53.10)
COP = ("https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N53_00_{0}_00_DEM/"
       "Copernicus_DSM_COG_10_N53_00_{0}_00_DEM.tif")
tiles = [COP.format("W005"), COP.format("W004")]
mosaic, transform = merge(tiles, bounds=BBOX)
with rasterio.open(tiles[0]) as src:
    profile = src.profile | {"driver": "GTiff", "width": mosaic.shape[2], "height": mosaic.shape[1],
                             "transform": transform, "compress": "deflate"}
with rasterio.open("snowdonia_glo30.tif", "w", **profile) as dst:
    dst.write(mosaic)
dem = mosaic[0].astype("float64")
print(f"snowdonia_glo30.tif {dem.shape[1]} × {dem.shape[0]}, {dem.min():.0f} to {dem.max():.0f} m")

snowdonia_glo30.tif 240 × 252, 54 to 1075 m

Example 1 — the full chain with both diagnoses reported

import math

import numpy as np
import rasterio
from scipy.ndimage import label, minimum_filter
from skimage.morphology import reconstruction

OFFSETS = [(-1, 0), (-1, 1), (0, 1), (1, 1), (1, 0), (1, -1), (0, -1), (-1, -1)]
CODES = [64, 128, 1, 2, 4, 8, 16, 32]


def drainage_counts(surface):
    neighbour_min = minimum_filter(surface, size=3, mode="nearest")
    interior = np.zeros(surface.shape, bool)
    interior[1:-1, 1:-1] = True
    strict = (surface < neighbour_min) & interior
    stuck = (surface <= neighbour_min) & interior
    return {"pits": int(strict.sum()), "pit_regions": int(label(strict)[1]),
            "stuck": int(stuck.sum()), "stuck_regions": int(label(stuck)[1])}


def hydrology(path, *, stream_cells=200):
    with rasterio.open(path) as src:
        dem = src.read(1).astype("float64")
        if src.nodata is not None:
            dem = np.where(dem == src.nodata, np.nan, dem)
        if src.crs.is_geographic:
            lat = (src.bounds.bottom + src.bounds.top) / 2
            cell_x = abs(src.transform.a) * 111_320 * math.cos(math.radians(lat))
            cell_y = abs(src.transform.e) * 110_574
        else:
            cell_x, cell_y = abs(src.transform.a), abs(src.transform.e)
        profile = src.profile.copy()

    if np.isnan(dem).any():
        raise ValueError("fill voids before hydrology — NaN breaks flow routing")

    before = drainage_counts(dem)
    print(f"  before: {before['pits']:,} true pits, "
          f"{before['stuck']:,} cells with no downhill neighbour "
          f"({before['stuck_regions']:,} regions)")

    seed = np.full_like(dem, dem.max())
    seed[0, :], seed[-1, :] = dem[0, :], dem[-1, :]
    seed[:, 0], seed[:, -1] = dem[:, 0], dem[:, -1]
    filled = reconstruction(seed, dem, method="erosion")

    after = drainage_counts(filled)
    if after["pits"]:
        raise RuntimeError(f"{after['pits']} true pits survived filling")
    raised = filled - dem
    print(f"  filled: {(raised > 0).sum():,} cells raised, "
          f"median {np.median(raised[raised > 0]):.2f} m, max {raised.max():.2f} m")
    print(f"  after:  {after['pits']:,} true pits, {after['stuck']:,} on flats "
          f"({after['stuck_regions']:,} regions) — filling converts pits into flats")

    best_drop = np.full(filled.shape, -np.inf)
    direction = np.zeros(filled.shape, dtype="uint8")
    for (dr, dc), code in zip(OFFSETS, CODES):
        shifted = np.roll(np.roll(filled, -dr, axis=0), -dc, axis=1)
        drop = (filled - shifted) / np.hypot(dr * cell_y, dc * cell_x)
        better = drop > best_drop
        best_drop = np.where(better, drop, best_drop)
        direction = np.where(better, code, direction).astype("uint8")
    direction[best_drop <= 0] = 0
    direction[0, :] = direction[-1, :] = direction[:, 0] = direction[:, -1] = 0
    print(f"  {(direction == 0).sum() - 2 * (sum(filled.shape) - 2) - 4:,} interior cells "
          f"have no assigned direction (flats)")

    rows, cols = filled.shape
    order = np.argsort(filled.ravel())[::-1]
    accumulation = np.ones(filled.size)
    lookup = dict(zip(CODES, OFFSETS))
    flat_dir = direction.ravel()
    for idx in order:
        code = flat_dir[idx]
        if code == 0:
            continue
        dr, dc = lookup[code]
        r, c = divmod(int(idx), cols)
        nr, nc = r + dr, c + dc
        if 0 <= nr < rows and 0 <= nc < cols:
            accumulation[nr * cols + nc] += accumulation[idx]
    accumulation = accumulation.reshape(filled.shape)

    cell_km2 = cell_x * cell_y / 1e6
    print(f"  accumulation: max {accumulation.max():,.0f} cells "
          f"({accumulation.max() * cell_km2:.1f} km²)")
    return filled, direction, accumulation, (cell_x, cell_y), profile


filled, direction, accumulation, cell, profile = hydrology("snowdonia_glo30.tif")

  before: 0 true pits, 1,394 cells with no downhill neighbour (131 regions)
  filled: 2,058 cells raised, median 1.07 m, max 10.51 m
  after:  0 true pits, 2,068 on flats (115 regions) — filling converts pits into flats
  2,064 interior cells have no assigned direction (flats)
  accumulation: max 3,866 cells (3.3 km²)

That maximum accumulation is the largest single basin the flow graph found: 3,866 cells, or 3.3 km². On an 8 × 7 km extent that is small, and the reason is exactly the flats — 115 regions where flow has no assigned direction, each truncating the paths that reach it.

This is the number to watch. A conditioned DEM with unresolved flats produces many small basins; a properly flat-resolved one produces a few large ones.

Example 2 — choosing the stream threshold

import pandas as pd
from scipy.ndimage import label as connected


def stream_threshold_sweep(accumulation, cell, thresholds=(50, 100, 200, 500, 1000, 5000)):
    cell_km2 = cell[0] * cell[1] / 1e6
    rows = []
    for t in thresholds:
        streams = accumulation >= t
        rows.append({
            "cells": t,
            "area_km2": round(t * cell_km2, 3),
            "share": f"{streams.mean():.2%}",
            "segments": int(connected(streams)[1]),
        })
    print(pd.DataFrame(rows).to_string(index=False))


stream_threshold_sweep(accumulation, cell)

 cells  area_km2 share  segments
    50     0.043 6.58%      1320
   100     0.086 3.47%       786
   200     0.171 1.83%       416
   500     0.428 0.78%       209
  1000     0.856 0.37%       100
  5000     4.280 0.00%         0

Read the segments column as a warning rather than a result. A real drainage network over this area is a handful of connected systems, not 416 fragments — and at 5,000 cells nothing survives at all, because no basin in the flow graph is that large.

That pattern, many segments at every threshold and nothing above a moderate one, is the signature of unresolved flats. The area_km2 column is still the right way to express the threshold: "a channel begins where 0.17 km² drains into it" is defensible; "200 cells" is meaningless at a different resolution.

Example 3 — delineating a watershed, and checking it

def watershed(direction, outlet):
    """Walk the flow graph upstream from one cell."""
    rows, cols = direction.shape
    basin = np.zeros(direction.shape, bool)
    basin[outlet] = True
    stack = [outlet]

    while stack:
        r, c = stack.pop()
        for (dr, dc), code in zip(OFFSETS, CODES):
            nr, nc = r - dr, c - dc                   # cells that drain INTO (r, c)
            if not (0 <= nr < rows and 0 <= nc < cols) or basin[nr, nc]:
                continue
            if direction[nr, nc] == code:
                basin[nr, nc] = True
                stack.append((nr, nc))
    return basin


outlet = tuple(int(v) for v in np.unravel_index(accumulation.argmax(), accumulation.shape))
basin = watershed(direction, outlet)

cell_km2 = cell[0] * cell[1] / 1e6
print(f"outlet {outlet}, accumulation {accumulation[outlet]:,.0f}")
print(f"basin {basin.sum():,} cells = {basin.sum() * cell_km2:.2f} km²")
print(f"consistent: {abs(basin.sum() - accumulation[outlet]) < 1}")

outlet (216, 124), accumulation 3,866
basin 3,866 cells = 3.31 km²
consistent: True

The last line is a real check, not decoration. Flow accumulation at a cell is by definition the number of cells draining through it, so an independent upstream walk must return exactly that count.

If the two disagree, the direction grid contains a cycle or an inconsistency — which is what an incomplete fill, or a direction assigned into a flat, produces. It is worth asserting.

Filling and flat resolution are separate problems

Why filling and flat-resolution are separate problems

Filling guarantees no cell is strictly lower than all its neighbours. It says nothing about cells that are equal to a neighbour.

On a filled depression every cell is at the outlet level — a plateau. D8 asks "which neighbour gives the steepest descent" and the answer is none, because the descent is zero in every direction. The cell gets direction 0 and flow stops.

The measured effect on this DEM: filling raised 2,058 cells, and the count of cells with nowhere to drain rose from 1,394 to 2,068. The problem did not go away; it changed shape.

Proper flat-resolution imposes a small artificial gradient across each flat — away from higher ground and toward the outlet — so that D8 has something to follow. That is the step hand-rolled implementations usually omit, and it is why pysheds, richdem and whitebox produce connected networks where a fifty-line NumPy version produces fragments.

Why zero true pits is a clue about your DEM

A raw photogrammetric or radar DEM is noisy and has thousands of genuine single-cell pits. Copernicus DEM has none, because it is a conditioned product — voids filled, noise smoothed, and in places already hydrologically enforced.

So the strict test returning zero is informative rather than reassuring. It tells you the DEM has been processed, that filling will change little, and that flats are where your effort should go.

Run both tests. The strict one characterises the DEM; the <= one tells you what will actually break.

Why D8 is crude and still standard

D8 sends all of a cell's flow to one neighbour. Real water spreads, especially on hillslopes and fans, and D8 cannot represent that — a slope draining south-south-east becomes a zigzag of south and south-east steps.

D-infinity apportions flow between two neighbours by angle; MFD spreads it among all downhill neighbours. Both give smoother, more realistic accumulation.

D8 survives because it gives each cell exactly one downstream neighbour, which makes the flow graph a forest. That is what allows the single-pass accumulation in Example 1 and the simple upstream walk in Example 3. With MFD neither works, and watershed delineation needs a fundamentally different algorithm.

Why the diagonal correction matters

Comparing raw height drops rather than gradients biases every decision toward diagonals: a diagonal neighbour is √2 further away, so it can drop 41% more without being steeper.

On a uniform planar slope this produces alternating diagonal and orthogonal steps in a regular pattern — a herringbone across the whole hillside, visible immediately in the direction raster and inherited by every stream. Dividing by hypot(dr * cell_y, dc * cell_x) fixes it, and also handles the non-square cells that a geographic DEM produces.

Both pit tests, wrapped edges and diagonal distance

• Run both pit tests. Strictly-lower finds true pits; less-than-or-equal finds everything that cannot drain, which is what actually breaks flow.
• Filling increases the flat area. Expect the "no downhill neighbour" count to go up, not down.
• Never use the filled DEM for elevation, slope or contours. It is a hydrological working surface.
• Real closed basins are destroyed by filling. Many single-cell depressions are noise; a few large ones may be real.
• NaN breaks flow routing entirely. Fill voids first.
• np.roll wraps at the edges. Mask the boundary rows and columns.
• Divide diagonal drops by the true distance, or diagonals win and you get a herringbone.
• pysheds, richdem and whitebox implement flat-resolution and D-infinity properly. For real hydrology, use them — this article is for understanding what they do.

Internal links

• Digital elevation models explained — why a conditioned DEM has no true pits
• How to calculate slope and aspect from a DEM in Python — use the unfilled DEM for this
• Slope values are wrong or absurdly steep — NoData contamination, which also breaks routing
• How to make a hillshade from a DEM in Python — for judging whether a depression is real
• How to generate contour lines from a DEM in Python — again, from the unfilled DEM
• The raster data model explained — NoData and dtype
• How to clean line network dangles in Python — the vector equivalent of a broken network
• How to calculate zonal statistics in Python — summarising by delineated watersheds
• How to condition a DEM for hydrology: breaching or filling — the alternative to filling
• How to delineate a watershed from a pour point — the next step after flow direction

FAQ

What is the difference between a pit and a flat?

A pit is strictly lower than all eight neighbours. A flat has at least one neighbour at exactly the same height. Filling removes pits and creates flats.

My DEM has zero pits. Do I still need to fill it?

Probably not, but check the other test. Copernicus DEM has zero true pits and 1,394 cells with no downhill neighbour — all of them flats, which filling will not help with.

Why did filling make the problem worse?

It did not make it worse, it changed its shape. Every filled depression becomes a plateau at the outlet level, so the count of cells with nowhere to drain rises while the number of separate problem regions falls.

Why is my stream network fragmented?

Unresolved flats. Each one truncates the flow paths reaching it. Use a library that implements flat-resolution — pysheds, richdem or whitebox.

What stream threshold should I use?

Express it as a contributing area in km² rather than a cell count, and calibrate against a mapped network. A cell count is meaningless at a different resolution.

Why divide diagonal drops by √2?

Because diagonal neighbours are further away. On the Snowdon DEM, comparing raw drops sends 90% of cells down a diagonal against 49% by gradient, and changes the direction of 39% of cells — a herringbone flow pattern.

Can I use the filled DEM for anything else?

No. Keep the original for slope, contours and any elevation reporting. The filled surface exists only to make flow routing possible.