Python LiDAR

Hydro-Flattening Water Bodies in a DTM

TL;DR: For each lake polygon, set every DTM cell inside to a low percentile of the elevations in a thin shoreline ring. For each river polygon, sample bank elevations along the centreline, force the profile to be non-increasing downstream with a running minimum, and interpolate each cell’s water surface from its position along the centreline. Burn both with rasterio masks, then verify that no bank cell sits below the adjacent water.

# Context and Motivation

This guide is part of Water and Bridge Classification in LiDAR. A DTM interpolated across water is a mess: sparse returns from waves and specular stripes produce bumps and pits, and void-filling across empty areas draws arbitrary slopes. Hydro-flattening replaces all of that with a surface a hydrologist would recognize — lakes perfectly flat, rivers descending smoothly, and banks that meet the water cleanly. The USGS Lidar Base Specification requires it for delivered DEMs, and any flood or drainage model run on an un-flattened DTM will route water uphill somewhere.

The two cases need different treatment. A lake has one surface elevation. A river’s surface falls downstream, so flattening it to one value would create a dam at the downstream end or a waterfall at the upstream end.

Flat for lakes, descending for rivers Three long profiles along a water body. The raw DTM surface is jagged with bumps and pits. A lake treatment replaces it with one horizontal line at the shoreline elevation. A river treatment replaces it with a line that steps down smoothly from upstream to downstream and never rises. raw lake one elevation river monotonic downstream upstream

# Prerequisites and Assumptions

  • A bare-earth DTM GeoTIFF built from class 2, e.g. with writers.gdal.
  • Water polygons with a kind attribute (lake or river), from water detection or supplied breaklines.
  • River centrelines as line strings digitized upstream to downstream, one per river polygon.
  • Python with rasterio, NumPy, Shapely 2.x and GeoPandas.

# Step-by-Step Implementation

# Step 1 — Build a shoreline ring per polygon

Buffer each polygon outward by 1.5 to 3 cells and subtract the polygon. The ring samples the banks, not the noisy water surface.

# Step 2 — Flatten lakes

Take the 5th percentile of DTM values in the ring and assign it to every cell inside the polygon.

# Step 3 — Build a river profile

Sample the centreline every 10 m. At each station, take a low percentile of ring cells within a short distance of the station — the local bank elevation.

# Step 4 — Enforce downstream monotonicity

Apply a running minimum from upstream to downstream: each station’s water elevation is the smaller of its own sample and the previous station’s value.

# Step 5 — Burn river cells

For each cell inside the river polygon, find its projected distance along the centreline and interpolate the water elevation from the station profile.

# Step 6 — Verify banks

Check that ring cells are not lower than the water surface next to them by more than a few centimetres.

# Complete Working Example

python
"""Hydro-flatten lakes (single level) and rivers (monotonic profile) in a DTM."""
from __future__ import annotations

from pathlib import Path

import geopandas as gpd
import numpy as np
import rasterio
from rasterio.features import geometry_mask
from shapely.geometry import LineString, Point


def ring_mask(poly, shape, transform, width):
    ring = poly.buffer(width).difference(poly)
    return geometry_mask([ring], shape, transform, invert=True)


def flatten_lake(dtm, poly, transform, width, pct=5.0):
    shore = ring_mask(poly, dtm.shape, transform, width)
    z = float(np.nanpercentile(dtm[shore], pct))
    inside = geometry_mask([poly], dtm.shape, transform, invert=True)
    dtm[inside] = z
    return z


def flatten_river(dtm, poly, centre: LineString, transform, width, step=10.0, pct=10.0):
    shore = ring_mask(poly, dtm.shape, transform, width)
    rows, cols = np.nonzero(shore)
    xs, ys = rasterio.transform.xy(transform, rows, cols)
    shore_xy = np.column_stack([xs, ys])
    shore_z = dtm[rows, cols]

    stations = np.arange(0.0, centre.length + step, step)
    levels = np.full(len(stations), np.nan)
    for i, s in enumerate(stations):
        p = centre.interpolate(s)
        near = np.hypot(shore_xy[:, 0] - p.x, shore_xy[:, 1] - p.y) < 2 * step
        if near.any():
            levels[i] = np.nanpercentile(shore_z[near], pct)
    good = ~np.isnan(levels)
    levels = np.interp(stations, stations[good], levels[good])
    levels = np.minimum.accumulate(levels)          # never rises downstream

    inside = geometry_mask([poly], dtm.shape, transform, invert=True)
    r, c = np.nonzero(inside)
    cx, cy = rasterio.transform.xy(transform, r, c)
    along = np.array([centre.project(Point(x, y)) for x, y in zip(cx, cy)])
    dtm[r, c] = np.interp(along, stations, levels)
    return levels


def hydro_flatten(dtm_in: Path, water_gpkg: Path, centre_gpkg: Path, dtm_out: Path) -> None:
    with rasterio.open(dtm_in) as ds:
        profile, transform = ds.profile, ds.transform
        dtm = ds.read(1, masked=True).filled(np.nan).astype("float64")
        width = 2.0 * ds.res[0]
    water = gpd.read_file(water_gpkg)
    centres = gpd.read_file(centre_gpkg).set_index("water_id")
    for row in water.itertuples():
        if row.kind == "lake":
            z = flatten_lake(dtm, row.geometry, transform, width)
            print(f"lake {row.water_id}: {z:.2f} m")
        else:
            lv = flatten_river(dtm, row.geometry, centres.loc[row.water_id].geometry,
                               transform, width)
            print(f"river {row.water_id}: {lv[0]:.2f} → {lv[-1]:.2f} m")
    profile.update(dtype="float32", nodata=-9999.0, compress="deflate")
    with rasterio.open(dtm_out, "w", **profile) as out:
        out.write(np.nan_to_num(dtm, nan=-9999.0).astype("float32"), 1)


if __name__ == "__main__":
    hydro_flatten(Path("dtm_valley.tif"), Path("water.gpkg"), Path("centrelines.gpkg"),
                  Path("dtm_valley_hydro.tif"))
The running minimum removes uphill steps Dots show bank elevation samples at stations along a river, generally falling downstream but with two upward bumps where a bank sample caught a levee or vegetation. A solid line shows the running minimum, which follows the samples down and stays flat across each bump instead of rising. levee sample ignored vegetation sample ignored distance downstream along centreline elevation

# Key Parameter Table

Parameter Type Default Guidance
ring width cells 2 1–3 cells; wider picks up bank slopes, narrower picks up water noise
lake percentile float 5 Low enough to sit at or below every bank cell
river station spacing float, m 10 5–25 m depending on river size and gradient
river bank percentile float 10 Slightly higher than lakes, since banks vary more along a river
station search radius float, m 2 × spacing Must reach both banks on the widest reach

# Verification

  • Banks above water. For every ring cell, the bank should not be more than about 5 cm below the flattened water next to it. More than that creates a “moat” in hillshades.
  • Monotonic rivers. Sample the burned DTM along each centreline and assert the profile never increases.
  • Lakes exactly flat. The standard deviation of cells inside each lake polygon is zero.
python
with rasterio.open("dtm_valley_hydro.tif") as ds:
    dtm = ds.read(1, masked=True)
    for row in water[water.kind == "lake"].itertuples():
        inside = geometry_mask([row.geometry], dtm.shape, ds.transform, invert=True)
        assert float(dtm[inside].std()) == 0.0, f"lake {row.water_id} not flat"

# Gotchas and Edge Cases

Centreline direction. The running minimum assumes the centreline is digitized upstream to downstream. Reversed lines produce a profile that is flat from the source to the mouth at the mouth’s elevation. Check direction against the DTM before flattening.

Lakes inside rivers. A reservoir behind a dam is a lake; the river above and below it is not. Split the water polygon at the dam and treat each part appropriately.

Tile boundaries. Flattening tile by tile gives each tile its own lake level and its own river profile. Always flatten on a mosaic or on polygons processed as whole features across tiles.

One lake, two tiles, two levels A lake crosses a vertical tile boundary. Flattened per tile, the left half gets 212.4 metres and the right half 212.9 metres, producing a visible 0.5 metre step at the boundary. Flattened as one feature on a mosaic, both halves share 212.4 metres. 212.4 m 212.9 m 0.5 m step at the tile edge — flatten whole features instead

Islands. Water polygons with holes need the holes excluded from both the ring (they are banks) and the fill. geometry_mask honours interior rings, so pass the polygon with its holes intact.

# Frequently Asked Questions

What elevation should a flattened lake have?

At or just below the lowest part of its shoreline, so water never appears to sit above the land. A low percentile of bank cells achieves that while ignoring a few anomalously low cells.

Why can’t rivers be flattened to a single elevation?

Because rivers flow downhill. A single level would either flood the upstream reach or leave a waterfall at the downstream end. Rivers need a surface that is flat across the channel and descends along it.

What size of water body needs hydro-flattening?

It depends on the specification. The USGS Lidar Base Specification applies it to ponds and lakes of about two acres or more and to streams about 100 feet wide or more. Contracts often set their own thresholds.

Does hydro-flattening change the point cloud?

No. It is applied to the raster DEM. The points are classified as water but keep their measured elevations.

Related in Water and Bridge Classification in LiDAR

↑ Water and Bridge Classification in LiDAR overview --- title: "Water and Bridge Classification in LiDAR" description: "Classify water to ASPRS class 9 and bridge decks to class 17, then... Classifying Bridge Decks --- title: "Classifying Bridge Decks" description: "Find bridge decks that ground filters mistook for terrain and reclassify them to ASPRS... Classifying Water from Intensity and Returns --- title: "Classifying Water from Intensity and Returns" description: "Detect water in topographic LiDAR without supplied polygons:...