🔧 Fix Hydrology & watersheds Terrain & elevation Hydrology & Watersheds · step 17 of 20 Updated

Fixing Streams That Run in Straight Parallel Lines Across Flat Ground

Problem statement

A stream network derived from a DEM looks natural on hillsides and then, across a lake, a reservoir or a wide valley floor, turns into a comb of straight, parallel lines. None of them follows a real channel. The DEM has no slope there — lakes are hydro-flattened to one elevation, floodplains are nearly level — so flow direction is decided by how a tool resolves flat areas, not by the terrain. Different tools resolve flats very differently.

Measured on a synthetic reservoir and on the largest flat area in the 10 m 3DEP DEM around Esopus Creek, New York:

  • The DEM held 220,358 cells with no relief within their 3 × 3 neighbourhood, in 65 flat areas; the largest was a 15.45 km² water surface at 177.30 m.
  • WhiteboxTools' fill_depressions with fix_flats=True drew 14 separate channels of at least 1 km² across the middle of that flat, on 3,948 channel cells.
  • pysheds' resolve_flats drew one, on 606 cells — but conditioning took 20.3 s against 2.4 s.
  • Burning the mapped NHDPlus flow paths 1 m into the flat before filling cut the 14 channels to 6. On a synthetic flat reservoir floor, a cross-slope of 0.1% towards the centre turned 14 parallel channels into one.

Quick answer 

cause                                   fix
flat area resolved towards many exits   resolve flats towards lower terrain as well as away from higher (pysheds resolve_flats)
lake or reservoir surface               burn the mapped flow path, or mask the water body and route around it
nearly level valley floor               add a gentle gradient towards the mapped channel, or burn it
coarse or smoothed DEM                  use a finer bare-earth DEM

Count the channels crossing a line drawn across the flat. One or two is a river; ten is an artefact.

Bar chart of the number of derived channels crossing the middle of the largest flat area in the Esopus DEM after WhiteboxTools fill with flat fixing, pysheds resolve_flats, and burning mapped flow paths before filling.
Same DEM, same flat: the flat-resolution method decided whether one channel crossed it or fourteen.

Step-by-step solution

1. Find the flats

A cell whose 3 × 3 neighbourhood spans less than a centimetre of elevation has no direction of its own. Label connected flat cells and rank the areas (Example 2). Large flats are almost always water bodies flattened in the DEM, reservoirs, or filled depressions behind embankments.

2. Count the channels crossing each large flat

Take a row or column through the middle of the flat and count how many separate runs of channel cells cross it. On the Esopus reservoir surface the middle row was crossed by 14 channels after WhiteboxTools conditioning; a real river crosses once.

3. Understand why straight lines appear

Flat-resolution algorithms impose a tiny artificial gradient so every flat cell drains. If the gradient only points away from higher ground, or towards the nearest outlet along one axis, neighbouring cells get the same direction and flow runs in parallel straight lines until it reaches the edge of the flat. D8 then keeps those lines separate, because flow never merges sideways on a surface with no cross-gradient.

4. Prefer a flat resolution that converges

pysheds' resolve_flats combines a gradient away from higher terrain with a gradient towards lower terrain, which draws flow across the flat towards its outlet. On the synthetic reservoir it produced a single channel across the floor, where WhiteboxTools' fix produced 14 (Example 1). On the real flat: one channel against 14.

5. Or give the flat the gradient it lacks

Where you know where the channel should be, tell the DEM. A cross-slope of just 0.1% towards the centre of the synthetic reservoir floor was enough for WhiteboxTools to draw one channel. On real data, burning mapped flow paths — NHDPlus artificial paths run through water bodies — lowered 2,794 cells of the Esopus flat and reduced the crossing channels from 14 to 6 (Example 3); see burning a river network into a DEM.

6. Mask water bodies when the path across does not matter

For catchment delineation, flow across a reservoir only needs to reach its outlet. Masking the lake and treating it as a single node, or clipping stream products to land, removes the artefact from maps without changing catchment areas.

7. Check the cost of the conditioning you choose

pysheds' fill-and-resolve sequence took 20.3 s on the 12.5-million-cell DEM against 2.4 s for WhiteboxTools, and on this DEM it left residual pits that need a second pass for delineation — see fixing a watershed that comes out as a few pixels. A converging flat resolution is not automatically the best choice for every product.

8. Clip or flag stream lines on flats in the final product

Even a single derived line across a lake is a routing path, not a channel. Flag network segments that lie on flats or water bodies, or replace them with the mapped path, before publishing stream maps or computing stream lengths.

Bar chart of channels crossing the middle of a synthetic flat reservoir floor for WhiteboxTools flat fixing, WhiteboxTools with a 0.1 percent cross-slope, and pysheds resolve_flats.
A perfectly flat floor has no answer to where water goes; a 0.1% slope supplies one.

Code examples

Example 1 — a synthetic reservoir with a perfectly flat floor 

import os
import time

import numpy as np
import rasterio
import whitebox
from rasterio.transform import from_origin
from scipy import ndimage

wbt = whitebox.WhiteboxTools()
wbt.set_verbose_mode(False)
wbt.set_working_dir(os.path.abspath("."))


def reservoir(cross_slope=0.0):
    """300 x 201 cells of 10 m: a floor 101 cells wide behind a dam with a one-cell spillway at column 100."""
    rows, cols = np.mgrid[0:300, 0:201]
    across = np.abs(cols - 100)
    surface = 100.0 + cross_slope * across * 10 + np.where(across > 50, (across - 50) * 0.5, 0.0)
    surface = surface + np.where(rows < 20, (20 - rows) * 0.5, 0.0)
    surface[270:, :] = np.maximum(surface[270:, :], 105.0)
    surface[270:, 100] = 99.9 - np.arange(30) * 0.1
    return surface.astype("float32")


def channels_crossing(acc, row, threshold):
    line = acc[row] >= threshold
    return int(np.count_nonzero(line[1:] & ~line[:-1]) + line[0])


profile = dict(driver="GTiff", width=201, height=300, count=1, dtype="float32", crs="EPSG:32618",
               transform=from_origin(500_000, 4_600_000, 10, 10), nodata=-9999.0)
for label, surface in (("flat floor", reservoir()), ("floor sloping 0.1% to the centre", reservoir(0.001))):
    with rasterio.open("reservoir.tif", "w", **profile) as dst:
        dst.write(surface, 1)
    wbt.fill_depressions("reservoir.tif", "reservoir_filled.tif", fix_flats=True)
    wbt.d8_pointer("reservoir_filled.tif", "reservoir_d8.tif")
    wbt.d8_flow_accumulation("reservoir_d8.tif", "reservoir_acc.tif", out_type="cells", pntr=True)
    with rasterio.open("reservoir_acc.tif") as src:
        acc = src.read(1)
    print(f"WhiteboxTools, {label:32} channels (>= 200 cells) crossing row 150: {channels_crossing(acc, 150, 200):3d}; "
          f"at the spillway {acc[299, 100]:,.0f} cells")

if not hasattr(np, "in1d"):
    np.in1d = lambda a, b, **kw: np.isin(np.ravel(a), b, **kw)
from pyproj import Proj
from pysheds.grid import Grid
from pysheds.view import Raster, ViewFinder

view = ViewFinder(affine=profile["transform"], shape=(300, 201), crs=Proj("EPSG:32618"), nodata=-9999.0)
synthetic = Grid(viewfinder=view)
surface = synthetic.resolve_flats(synthetic.fill_depressions(Raster(reservoir().astype("float64"), viewfinder=view)))
acc = np.asarray(synthetic.accumulation(synthetic.flowdir(surface)), dtype="float64")
print(f"pysheds resolve_flats, {'flat floor':25} channels (>= 200 cells) crossing row 150: {channels_crossing(acc, 150, 200):3d}; "
      f"at the spillway {acc[299, 100]:,.0f} cells")
WhiteboxTools, flat floor                       channels (>= 200 cells) crossing row 150:  14; at the spillway 60,300 cells
WhiteboxTools, floor sloping 0.1% to the centre channels (>= 200 cells) crossing row 150:   1; at the spillway 60,300 cells
pysheds resolve_flats, flat floor                channels (>= 200 cells) crossing row 150:   1; at the spillway 57,539 cells

pysheds needs float64 elevations here: a float32 array failed inside its compiled depression filling. Its spillway total was also lower, 57,539 of the grid's 60,300 cells, so not every cell it routed reached the spillway.

Example 2 — find the flat areas in a real DEM 

with rasterio.open("esopus_3dep13_utm.tif") as src:
    profile, dem = src.profile, src.read(1, masked=True)
    transform = src.transform
low = dem.filled(-1e6)
spread = ndimage.maximum_filter(low, size=3) - ndimage.minimum_filter(low, size=3)
flat = (spread < 0.01) & ~dem.mask
labels, count = ndimage.label(flat, structure=np.ones((3, 3)))
sizes = np.bincount(labels.ravel())
sizes[0] = 0
largest = labels == sizes.argmax()
rows, cols = np.nonzero(largest)
print(f"{int(flat.sum()):,} cells with no relief in their 3x3 window, in {count:,} areas; largest {sizes.max():,} cells "
      f"({sizes.max() / 1e4:.2f} km2) at {float(np.median(dem.data[largest])):.2f} m, rows {rows.min()}-{rows.max()}, cols {cols.min()}-{cols.max()}")
x_mid, y_mid = transform * (cols.mean(), rows.mean())
print(f"centre of the largest flat: x {x_mid:.0f}, y {y_mid:.0f}")
220,358 cells with no relief in their 3x3 window, in 65 areas; largest 154,476 cells (15.45 km2) at 177.30 m, rows 2674-3117, cols 2978-3901
centre of the largest flat: x 570871, y 4648077

Example 3 — compare flat handling on the real flat 

import geopandas as gpd
from rasterio.features import rasterize

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


def describe(label, acc, seconds):
    channel = (acc >= 10_000) & largest                       # 1 km2 of 10 m cells, on the flat only
    _, lines = ndimage.label(channel, structure=np.ones((3, 3)))
    middle = int(np.median(rows))
    crossing = channels_crossing(np.where(largest, acc, 0), middle, 10_000)
    print(f"{label:34} {seconds:5.1f} s; channel cells on the flat {int(channel.sum()):6,}; separate channel pieces {lines:3d}; "
          f"channels crossing row {middle}: {crossing}")


def whitebox_route(dem_path, tag):
    start = time.perf_counter()
    wbt.fill_depressions(dem_path, f"{tag}_filled.tif", fix_flats=True)
    wbt.d8_pointer(f"{tag}_filled.tif", f"{tag}_d8.tif")
    wbt.d8_flow_accumulation(f"{tag}_d8.tif", f"{tag}_acc.tif", out_type="cells", pntr=True)
    with rasterio.open(f"{tag}_acc.tif") as src:
        return src.read(1), time.perf_counter() - start


acc, seconds = whitebox_route("esopus_3dep13_utm.tif", "fill")
describe("WhiteboxTools fill, fix_flats", acc, seconds)

start = time.perf_counter()
grid = Grid.from_raster("esopus_3dep13_utm.tif")
raw = grid.read_raster("esopus_3dep13_utm.tif")
surface = grid.resolve_flats(grid.fill_depressions(grid.fill_pits(raw)))
fdir = grid.flowdir(surface)
describe("pysheds fill + resolve_flats", np.asarray(grid.accumulation(fdir)), time.perf_counter() - start)

flowlines = gpd.read_file("nhdplushr_flowlines_esopus.gpkg").to_crs(profile["crs"])
paths = flowlines[flowlines.ftype.isin([460, 558, 334])]
on_lines = rasterize(paths.geometry, out_shape=dem.shape, transform=transform).astype(bool) & ~dem.mask
print(f"mapped flow paths cross {int((on_lines & largest).sum()):,} cells of the flat")
with rasterio.open("burned.tif", "w", **profile) as dst:
    dst.write(np.where(on_lines, dem.data - 1.0, dem.data).astype(profile["dtype"]), 1)
acc, seconds = whitebox_route("burned.tif", "burn")
describe("mapped paths burned 1 m, then fill", acc, seconds)
WhiteboxTools fill, fix_flats        2.4 s; channel cells on the flat  3,948; separate channel pieces   3; channels crossing row 2860: 14
pysheds fill + resolve_flats        20.3 s; channel cells on the flat    606; separate channel pieces   1; channels crossing row 2860: 1
mapped flow paths cross 2,794 cells of the flat
mapped paths burned 1 m, then fill   2.3 s; channel cells on the flat  1,920; separate channel pieces   2; channels crossing row 2860: 6

Fourteen channels crossed the middle row yet formed only three connected pieces on the flat, so the parallel lines joined somewhere before leaving it.

Explanation

Why flats have no flow direction

D8 sends each cell's flow to its steepest downslope neighbour. On a flat, every neighbour is level, so there is no steepest one. Every tool must invent a direction, and the invention is a rule about distance to higher or lower edges, not a measurement.

Why some flat fixes make parallel lines

A gradient built only from distance to the flat's outlet edge assigns the same direction to long rows of cells: each flows straight towards the edge and none has a reason to turn towards its neighbour. Accumulation then grows in parallel strips, each carrying the flow of its own strip of the flat.

Why a combined gradient converges

Adding a gradient away from higher terrain pushes flow off the flat's shores towards its middle, while the gradient towards lower terrain pushes it towards the outlet. Their sum points inwards and downstream, so strips merge into one path. That is the approach pysheds' resolve_flats takes, following Barnes and colleagues' method.

Why burning only partly helped

Burning lowers the cells under the mapped lines, but the rest of the flat is still level. Cells next to a burned line can drain into it; cells further away still get their direction from the flat fix. Six crossings remained on the Esopus flat after burning.

Two panels contrasting a flat resolved only towards its outlet edge, where flow runs in parallel strips, with a flat resolved towards the outlet and away from higher shores, where flow converges on one path.
Parallel lines are what a flat looks like when the imposed gradient has no sideways component.

Edge cases or notes

  • Hydro-flattened lakes are exactly flat by design; the artefact is guaranteed unless the lake is masked or its path burned.
  • Floodplains are nearly, not exactly, flat; a centimetre of noise can still produce parallel lines.
  • Filled depressions become flats, so aggressive filling creates new artefacts; breaching creates fewer.
  • Integer DEMs round gentle slopes into staircases of flats; use floating-point elevations.
  • D-infinity and MFD spread flow across flats rather than drawing parallel lines, which suits wetness indices more than stream maps.
  • Very large flats may take long to resolve; clip to the catchment first.
  • Water body masks from NHDPlus or national hydrography locate the largest flats directly, without scanning the DEM.

FAQ

Why do my DEM streams run in parallel straight lines?

The streams cross a flat area, and the tool's flat resolution sent neighbouring cells in the same direction. On the Esopus DEM, WhiteboxTools' flat fix drew 14 channels across one reservoir surface.

How do I fix parallel flow lines on flat areas?

Use a flat resolution that converges, such as pysheds' resolve_flats, burn the mapped flow path through the flat, or add a small gradient towards the known channel. pysheds drew one channel where WhiteboxTools drew 14.

Why do lakes create straight stream lines in a DEM?

Lakes are flattened to a single elevation in most DEMs, so flow direction across them is invented by the flat-resolution rule. Mask the lake or burn its mapped flow path.

Does fix_flats in WhiteboxTools remove parallel streams?

It makes every flat cell drain, but not necessarily towards one path: on a synthetic flat reservoir floor it produced 14 parallel channels.

Does stream burning fix flat areas?

Partly. Burning NHDPlus paths 1 m into the Esopus flat reduced 14 crossing channels to 6; cells away from the burned line still follow the flat fix.

How can I find flat areas in a DEM?

Mark cells whose 3 × 3 neighbourhood spans less than about a centimetre and label connected groups. The Esopus 10 m DEM had 65 such areas, the largest 15.45 km².