Python LiDAR

Water and Bridge Classification in LiDAR

Water and bridges are the two classes that decide whether a LiDAR terrain model is usable for hydrology. Near-infrared lasers are largely absorbed by water, so lakes and rivers come back as sparse, noisy returns — or no returns at all — and a DTM interpolated across them ripples with artefacts. Bridges are the opposite problem: they are solid, flat and connected to the road at both ends, so ground filters routinely classify decks as ground, and the resulting DTM dams every river at every crossing. Delivery specifications such as the USGS Lidar Base Specification require water classified as class 9, bridge decks as class 17, and water bodies above a size threshold flattened to a single elevation in the DEM. This topic, part of the classification and feature extraction section, covers all three.

Two classes, two opposite failures A profile across a river valley with a road bridge. Over the water, points are sparse and scattered around the surface, with gaps where the laser was absorbed. Over the bridge, points form a solid flat deck at road height. Beneath, a dashed line shows the DTM an unmodified ground filter produces: noisy over the water and raised across the channel where the deck was taken for ground, blocking flow. water surface bridge deck: dense, flat, road height naive DTM: dam across the channel water fails by having too few points, bridges by having the right points in the wrong class

# Prerequisites

  • PDAL 2.5+ with filters.overlay, filters.hag_nn, filters.covariancefeatures, filters.smrf and writers.gdal, plus the Python bindings.
  • Python 3.10+ with NumPy, GeoPandas, Shapely 2.x and rasterio.
  • Ground classified, with the understanding that some decks and some water returns will be in class 2 and must be corrected.
  • Water body and road polygons, if available. National hydrography layers, OpenStreetMap water and bridge ways, or a client breakline set make both classes far more reliable than detection from the points alone.
  • Point source IDs populated, because per-flightline behaviour over water (specular returns at nadir) matters for diagnostics.
  • A projected CRS in metres, the same one the breaklines use, or reprojection of the polygons before overlay.

# Core Workflow Architecture

  1. Find water candidates. Rasterize a return-count grid; cells inside large areas of very low or zero return density, adjacent to ground, are water candidates. Low intensity and flatness strengthen the evidence.
  2. Build water polygons. Polygonize the candidate cells, clean them morphologically, and merge with any supplied hydrography; drop polygons below the specification’s minimum size.
  3. Classify water. Burn polygons into the points with filters.overlay and set class 9 for ground and unclassified points inside them.
  4. Find bridge decks. Within road-over-water or road-over-road crossings, points classified as ground that sit well above the interpolated terrain on either side are deck candidates.
  5. Classify bridges. Set class 17 on deck points, removing them from the ground used for the DTM.
  6. Hydro-flatten. Build the DTM from class 2 only, then burn each water polygon at a single elevation derived from its shoreline.

# Full Implementation

python
"""Water (class 9) and bridge deck (class 17) classification, with hydro-flattening."""
from __future__ import annotations

import json
import logging
from pathlib import Path

import geopandas as gpd
import numpy as np
import pdal
import rasterio
import rasterio.features
from rasterio.transform import from_origin
from scipy import ndimage as ndi
from shapely.geometry import shape

log = logging.getLogger("hydro")


def water_candidates(src: Path, cell: float = 2.0, min_area_m2: float = 8000.0,
                     crs: str = "EPSG:6341") -> gpd.GeoDataFrame:
    """Polygons of near-empty cells large enough to count as water bodies."""
    p = pdal.Pipeline(json.dumps({"pipeline": [str(src)]}))
    p.execute()
    a = p.arrays[0]
    x0, y1 = np.floor(a["X"].min()), np.ceil(a["Y"].max())
    cols = int(np.ceil((a["X"].max() - x0) / cell))
    rows = int(np.ceil((y1 - a["Y"].min()) / cell))
    c = np.clip(((a["X"] - x0) / cell).astype(int), 0, cols - 1)
    r = np.clip(((y1 - a["Y"]) / cell).astype(int), 0, rows - 1)
    counts = np.zeros((rows, cols), dtype=np.int32)
    np.add.at(counts, (r, c), 1)

    expected = np.median(counts[counts > 0])
    empty = counts < max(1, 0.1 * expected)          # under 10 % of typical density
    empty = ndi.binary_opening(empty, iterations=2)  # drop isolated shadows
    empty = ndi.binary_closing(empty, iterations=3)  # join returns scattered on the surface
    transform = from_origin(x0, y1, cell, cell)
    polys = [shape(g) for g, v in rasterio.features.shapes(
        empty.astype(np.uint8), mask=empty, transform=transform) if v == 1]
    gdf = gpd.GeoDataFrame({"geometry": polys}, crs=crs)
    gdf = gdf[gdf.area >= min_area_m2].reset_index(drop=True)
    gdf["WaterId"] = np.arange(1, len(gdf) + 1)
    log.info("%d water candidates >= %.0f m²", len(gdf), min_area_m2)
    return gdf


def classify(src: Path, dst: Path, water: gpd.GeoDataFrame, crs: str = "EPSG:6341") -> None:
    water_gpkg = dst.with_suffix(".water.gpkg")
    water.to_file(water_gpkg, layer="water", driver="GPKG")
    stages = [
        {"type": "readers.las", "filename": str(src), "override_srs": crs},
        {"type": "filters.ferry", "dimensions": "=>WaterId"},
        {"type": "filters.overlay", "dimension": "WaterId",
         "datasource": str(water_gpkg), "layer": "water", "column": "WaterId"},
        {"type": "filters.assign", "value": [
            "Classification = 9 WHERE WaterId > 0 && (Classification == 1 || Classification == 2)"
        ]},
        # Bridge decks: ground-classified points far above the surrounding terrain.
        {"type": "filters.hag_nn", "count": 6, "max_distance": 60.0,
         "where": "Classification == 2"},
        {"type": "filters.assign", "value": [
            "Classification = 17 WHERE Classification == 2 && HeightAboveGround > 2.5"
        ]},
        {"type": "writers.las", "filename": str(dst), "minor_version": 4,
         "dataformat_id": 6, "forward": "all"},
    ]
    n = pdal.Pipeline(json.dumps({"pipeline": stages})).execute()
    log.info("%s: %d points written", dst.name, n)


def hydro_flatten(dtm_path: Path, water: gpd.GeoDataFrame, out_path: Path) -> None:
    """Burn each water polygon at the 5th percentile of its shoreline elevation."""
    with rasterio.open(dtm_path) as ds:
        dtm = ds.read(1, masked=True).filled(np.nan)
        profile, transform = ds.profile, ds.transform
    for poly in water.geometry:
        ring = poly.exterior.buffer(3.0).difference(poly)
        shore = rasterio.features.geometry_mask([ring], dtm.shape, transform, invert=True)
        inside = rasterio.features.geometry_mask([poly], dtm.shape, transform, invert=True)
        z = np.nanpercentile(dtm[shore], 5)
        dtm[inside] = z
    profile.update(nodata=-9999.0)
    with rasterio.open(out_path, "w", **profile) as out:
        out.write(np.nan_to_num(dtm, nan=-9999.0).astype(profile["dtype"]), 1)


if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO, format="%(levelname)s %(message)s")
    src = Path("valley_0310.laz")
    water = water_candidates(src)
    classify(src, Path("valley_0310_hydro.laz"), water)

# Code Breakdown

Return count as the water signal. A grid of return counts is the most reliable single piece of evidence for water in near-infrared data: open water returns a small fraction of the density of the land around it. Ten percent of the tile’s median cell count is a conservative threshold; wind-roughened water returns more, calm water less.

Opening then closing. Opening removes small empty patches — shadows behind buildings, dark roofs, missing returns under dense canopy — that are not water. Closing then fills the scattered surface returns inside a lake so the polygon is solid. The order matters: closing first would join shadows into lakes.

8,000 m² minimum. The USGS Lidar Base Specification requires hydro-flattening for inland ponds and lakes of about two acres or more, which is roughly 8,100 m²; it treats rivers above a nominal width of 30 m (100 ft) similarly. Use your contract’s thresholds, and prefer supplied breaklines to derived polygons where you have them.

filters.overlay into a scratch dimension. Overlay writes the polygon’s WaterId onto every point inside it. The dimension has to exist first, which is what filters.ferry with =>WaterId does. Keeping the ID lets you trace any point back to its water body later.

Only classes 1 and 2 become water. Points of overhanging vegetation, docks and boats inside the polygon keep their classes. Specifications differ on whether to relabel above-water points; the conservative default is to leave them.

Bridges from ground above ground. After water is removed from class 2, filters.hag_nn with a where clause computes, for each remaining ground point, its height above its nearest ground neighbours — with max_distance: 60 so the neighbours are drawn from the banks, not only from the deck itself. Deck points sit several metres above those neighbours. This heuristic works because SMRF has already classified most of the terrain correctly; the bridge deck guide shows a more robust version using road polygons.

Shoreline percentile for flattening. The water surface should sit at or just below the lowest bank, so the 5th percentile of DTM values in a 3 m ring around each polygon is used. Taking the minimum would chase one bad cell; the mean would put the water above parts of the bank.

The dropout that marks water A line chart of returns per square metre along a 600 metre transect. Over land the density hovers around 14. Over the lake between 220 and 430 metres it drops to between zero and one, with a brief spike near the centre where the laser hit the water at nadir. A dashed threshold line at ten percent of the land density sits just above the lake values. 10 % of land density lake nadir specular spike 14 0 distance along transect, 0–600 m (illustrative) returns/m²

The flattening step deserves a picture of its own, because the choice of where to sample the water elevation is what makes the shoreline look right or wrong in a hillshade.

One elevation from the shoreline ring Left: a plan view of a lake polygon surrounded by a three metre ring along its shore, where DTM cells are sampled. Right: a histogram of the sampled shoreline elevations ranging from 212.4 to 213.6 metres, with the fifth percentile at 212.5 marked as the chosen water surface, below almost every bank cell so the water never appears to sit above the land. water polygon 3 m shoreline ring sampled 5th percentile = water surface 212.4 m 213.6 m shoreline DTM elevations

# Parameter Reference Table

Parameter Type Default here Typical range Effect
cell (density grid) float, m 2.0 1–5 Coarser is steadier but blunts narrow channels
density threshold fraction 0.10 0.05–0.25 Fraction of median cell count below which a cell is empty
opening iterations int 2 1–4 Removes shadows and dark roofs from candidates
closing iterations int 3 2–6 Fills scattered surface returns inside lakes
min_area_m2 float 8000 per spec Smallest water body classified and flattened
hag_nn count int 6 3–12 Ground neighbours used to judge a deck point
hag_nn max_distance float, m 60 30–150 Must reach the banks across the widest deck
deck height float, m 2.5 1.5–5 Minimum clearance of a deck above the terrain
shoreline ring float, m 3.0 1–10 Width of the band sampled for water elevation

# Validation and Integrity Checks

  • Water polygons against imagery or hydrography. Overlay the derived polygons on orthophotos. False water usually appears as large flat dark roofs, fresh asphalt and deep shadows; false dry land appears where a lake was flown at an angle that returned strong specular points.
  • No ground inside water. After classification, count class 2 points inside each polygon; it should be zero, or only a handful along the edge where the polygon overlaps the bank.
  • Monotonic river surfaces. For flowing water, flattened elevations must decrease downstream. Sample the flattened DTM along the channel centreline and assert that each step is non-increasing. Lakes need only be flat; rivers need to be flat across and descending along.
  • Deck removal lets flow through. Run a quick flow accumulation on the DTM, or simply profile across each bridge: the channel must be continuous beneath every crossing.
python
def check_downstream(dtm: np.ndarray, rows: np.ndarray, cols: np.ndarray, tol: float = 0.02) -> None:
    profile = dtm[rows, cols]              # centreline cells ordered upstream to downstream
    rises = np.diff(profile) > tol
    assert not rises.any(), f"water surface rises at {int(rises.sum())} steps along the river"

# Performance Tuning

Water detection is raster work on a coarse grid and takes seconds per tile. The cost is in the classification pipeline, where filters.hag_nn with a large max_distance can be slow if it runs on every point — which is why it is restricted to class 2 with a where clause. On tiles with no water and no roads crossing water, skip the bridge stage entirely; a quick intersection of the tile extent with a road-over-water layer decides that in milliseconds.

Water bodies straddle tiles far more often than buildings, and flattening a lake tile by tile produces steps at every tile edge. Derive water polygons on a whole-project mosaic of the density grid, or merge per-tile polygons before computing shoreline elevations, so each lake gets exactly one elevation.

# Common Errors and Troubleshooting

Dark roofs and fresh asphalt classified as water. They absorb near-infrared too. Require that candidate cells are adjacent to ground, not surrounded by buildings, and filter polygons by compactness — roofs are rectilinear, lakes are not.

Rivers broken into pieces. Narrow channels under overhanging trees retain canopy returns, so the density grid does not drop to zero there. Use supplied centrelines or breaklines for rivers; detection from density alone is reliable only for open water.

Steps in the flattened river. Each polygon was flattened independently. For rivers, flatten along the channel with a gradient — interpolate water elevation from upstream to downstream shoreline samples — rather than to one value.

Deck classified but the DTM still dams. The DTM was built with triangulation across the gap, which bridges the deck’s absence with a straight line at deck height from the approaches. Build it with an interpolation that respects breaklines, or burn the channel from the water polygon before hillshading and hydrology.

Bridge approach embankments labelled as deck. Embankments are genuinely above nearby ground. Restrict the deck test to points over water or over another road, using the water polygons or a road layer as a mask.

# Frequently Asked Questions

Why do water bodies have so few LiDAR returns?

Topographic LiDAR uses near-infrared light, which water absorbs strongly. Most pulses produce no detectable return, except where the beam strikes the surface nearly vertically and reflects straight back, or where waves, foam or debris scatter it. Bathymetric sensors use green light to penetrate water, but that is a different instrument.

What is hydro-flattening?

It is the practice of setting each water body in a DEM to a single flat elevation, or for rivers to a surface that descends smoothly downstream, and making the shoreline consistent with the surrounding terrain. It removes the interpolation noise over water that would otherwise appear as bumps and pits.

Should bridges be classified as ground?

No. ASPRS class 17 exists for bridge decks precisely so that terrain models can exclude them. A DTM that includes decks shows a dam at every crossing, which breaks drainage analysis and flood modelling.

Do I need breaklines, or can I derive water from the points?

Large open lakes can be derived reliably from return density. Rivers, especially narrow or tree-lined ones, are much better handled with supplied hydrography or manually digitized breaklines. Most production work combines both: derive, then correct against supplied layers.

Hands-on guides

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:... Hydro-Flattening Water Bodies in a DTM --- title: "Hydro-Flattening Water Bodies in a DTM" description: "Flatten lakes to one elevation and rivers to a smooth downstream gradient...

More in LiDAR Classification and Feature Extraction in Python

Building Extraction from LiDAR Point Clouds --- title: "Building Extraction from LiDAR Point Clouds" description: "Classify roof points to ASPRS class 6 and turn... Machine Learning Point Classification for LiDAR --- title: "Machine Learning Point Classification for LiDAR" description: "Train a per-point classifier on PDAL-computed... Individual Tree Segmentation from LiDAR --- title: "Individual Tree Segmentation from LiDAR" description: "Split a forest canopy into individual trees:... Power Line Detection in LiDAR Point Clouds --- title: "Power Line Detection in LiDAR Point Clouds" description: "Find conductor points with linearity and height,... Point Cloud Segmentation with PDAL --- title: "Point Cloud Segmentation with PDAL" description: "Group points into objects with filters.cluster and...