🔧 Fix Hydrology & watersheds Terrain & elevation Hydrology & Watersheds · step 17 Notebook By · Updated

Fixing Streams That Run in Straight Parallel Lines Across Flat Ground

On this page

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 209,953 cells with no relief within their 3 × 3 neighbourhood, in 67 flat areas; the largest was a 15.19 km² water surface at 177.30 m.
  • WhiteboxTools' fill_depressions with fix_flats=True drew 17 separate channels of at least 1 km² across the middle of that flat, on 4,139 channel cells.
  • pysheds' resolve_flats drew one, on 600 cells — but conditioning took 19.6 s against 2.4 s.
  • Burning the mapped NHDPlus flow paths 1 m into the flat before filling cut the 17 channels to 7. 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.

Derived channels crossing the middle of the largest flat in the Esopus DEM by flat-handling method. channels crossing a 15.19 km2 flatWhiteboxTools fix_flats17 channels, 2.4 spaths burned 1 m7 channels, 2.8 spysheds resolve_flats1 channel, 19.6 s018Channels of at least 1 km2 crossing the flat's middle row, 3DEP 10 m spatialworkflow.io
Same DEM, same flat: the flat-resolution method decided whether one channel crossed it or fourteen.

Find the flats, count crossings, pick a converging fix

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 17 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 17.

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,651 cells of the Esopus flat and reduced the crossing channels from 17 to 7 (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 19.6 s on the 12.2-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.

Channels crossing the floor of a synthetic reservoir under three flat-resolution setups. synthetic reservoir, 300 x 201 cellssetupchannels crossing the floorWhitebox, perfectly flat floor14Whitebox, 0.1% slope on the floor1pysheds, perfectly flat floor1channels of at least 200 cells crossing row 150 of a 101-cell-wide floorA tenth of a percent of slope is the difference between one channel and fourteen spatialworkflow.io
The artefact needs a genuinely flat surface, which is why it shows up on water.

A flat-floored reservoir, synthetic and real

The data used on this page

The setup builds every input this page reads, with the same code as DEM hydrology explained: the basin of Esopus Creek above Coldbrook, New York (USGS 01362500), from the USGS NLDI; the same 3DEP data at 10 m over the basin plus 4 km (esopus_3dep13_utm.tif); NHDPlus HR flowlines from the USGS National Map.

import os
import time
from pathlib import Path

import geopandas as gpd
import numpy as np
import rasterio
import requests
from rasterio.merge import merge
from rasterio.warp import Resampling, reproject, transform_bounds


def get(url, **kw):
    """GET with patient retries: the USGS services answer 504 when busy, and recover."""
    for attempt in range(6):
        try:
            r = requests.get(url, timeout=600, **kw)
            r.raise_for_status()
            return r
        except requests.RequestException:
            if attempt == 5:
                raise
            time.sleep(30 * (attempt + 1))


# basins from the USGS NLDI
for site in ('01362500',):
    if not os.path.exists(f"nldi_basin_{site}.geojson"):
        Path(f"nldi_basin_{site}.geojson").write_text(get(f"https://api.water.usgs.gov/nldi/linked-data/nwissite/USGS-{site}/basin?f=json").text)
basin = gpd.read_file("nldi_basin_01362500.geojson")     # Esopus Creek at Coldbrook
W, S, E, N = basin.total_bounds + np.array([-0.03, -0.03, 0.03, 0.03])

TILES = {"3dep13": ([f"/vsicurl/https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/13/TIFF/current/{t}/USGS_13_{t}.tif"
                      for t in ("n43w075", "n42w075")], 10),
         "glo30": ([f"/vsicurl/https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_{t}_00_W075_00_DEM/"
                    f"Copernicus_DSM_COG_10_{t}_00_W075_00_DEM.tif" for t in ("N42", "N41")], 30)}


def dem_box(name, buffer_m, out):
    """UTM 18N grid of whole cells around the basin plus buffer_m, bilinear from the source tiles, no nodata."""
    if os.path.exists(out):
        return
    paths, res = TILES[name]
    srcs = [rasterio.open(p) for p in paths]
    mosaic, tf = merge(srcs, bounds=(W - 0.1, S - 0.1, E + 0.1, N + 0.1))
    bx = basin.to_crs(32618).total_bounds
    x0, y0 = np.floor((bx[0] - buffer_m) / res) * res, np.floor((bx[1] - buffer_m) / res) * res
    x1, y1 = np.ceil((bx[2] + buffer_m) / res) * res, np.ceil((bx[3] + buffer_m) / res) * res
    w, h = int((x1 - x0) / res), int((y1 - y0) / res)
    grid_tf = rasterio.transform.from_origin(x0, y1, res, res)
    grid = np.full((h, w), -9999, "float32")
    reproject(mosaic[0], grid, src_transform=tf, src_crs=srcs[0].crs, dst_transform=grid_tf, dst_crs="EPSG:32618",
              resampling=Resampling.bilinear, src_nodata=srcs[0].nodata, dst_nodata=-9999)
    with rasterio.open(out, "w", driver="GTiff", height=h, width=w, count=1, dtype="float32", crs="EPSG:32618",
                       transform=grid_tf, nodata=-9999, compress="deflate") as dst:
        dst.write(grid, 1)


dem_box("3dep13", 4000, "esopus_3dep13_utm.tif")         # the same, basin + 4 km
if not os.path.exists("nhdplushr_flowlines_esopus.gpkg"):  # NHDPlus HR flowlines, paged
    feats, off = [], 0
    while True:
        page = get("https://hydro.nationalmap.gov/arcgis/rest/services/NHDPlus_HR/MapServer/3/query", params={
            "geometry": f"{W},{S},{E},{N}", "geometryType": "esriGeometryEnvelope", "inSR": 4326,
            "spatialRel": "esriSpatialRelIntersects", "outFields": "*", "outSR": 4326, "f": "geojson",
            "resultOffset": off, "resultRecordCount": 2000}).json()
        feats += page["features"]
        if len(page["features"]) < 2000:
            break
        off += len(page["features"])
    flowlines = gpd.GeoDataFrame.from_features(feats, crs=4326)
    flowlines.columns = [c.lower() if c != "geometry" else c for c in flowlines.columns]
    flowlines.to_file("nhdplushr_flowlines_esopus.gpkg")

for f in ['nldi_basin_01362500.geojson', 'esopus_3dep13_utm.tif', 'nhdplushr_flowlines_esopus.gpkg']:
    if f.endswith(".tif"):
        with rasterio.open(f) as r:
            print(f"{f}: {r.width:,} x {r.height:,} cells of {r.res[0]:.6g} {'m' if r.crs.is_projected else 'degrees'}")
    elif f.endswith(".zip"):
        print(f"{f}: {os.path.getsize(f) / 1e6:.1f} MB")
    else:
        print(f"{f}: {len(gpd.read_file(f)):,} features")
nldi_basin_01362500.geojson: 1 features
esopus_3dep13_utm.tif: 4,071 x 3,008 cells of 10 m
nhdplushr_flowlines_esopus.gpkg: 4,739 features

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_max_procs(1)                      # one thread: parallel runs give slightly different results
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}")
209,953 cells with no relief in their 3x3 window, in 67 areas; largest 151,908 cells (15.19 km2) at 177.30 m, rows 2579-3007, cols 3046-3957
centre of the largest flat: x 570934, y 4648121

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  4,139; separate channel pieces   4; channels crossing row 2764: 17
pysheds fill + resolve_flats        19.6 s; channel cells on the flat    600; separate channel pieces   1; channels crossing row 2764: 1
mapped flow paths cross 2,651 cells of the flat
mapped paths burned 1 m, then fill   2.8 s; channel cells on the flat  2,379; separate channel pieces   4; channels crossing row 2764: 7

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

Why flats have no flow direction

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.

A flat resolved towards its outlet edge only produces parallel strips; a gradient from the shores as well converges flow. two ways to make a flat drainparallel strips✕imposed gradient has no sideways part✕neighbours share one direction✕synthetic floor: 14 channelsconverging flow✓gradient from shores and to outlet✓strips drawn together✓pysheds resolve_flats: 1 channelCounts from the synthetic reservoir and the Esopus flat spatialworkflow.io
Parallel lines are what a flat looks like when the imposed gradient has no sideways component.

Flattened lakes, floodplains and integer DEMs

  • 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 17 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 17.

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 17 crossing channels to 7; 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 67 such areas, the largest 15.19 km².