Python LiDAR

Classifying Water from Intensity and Returns

TL;DR: Grid the tile at 2 m and compute four per-cell cues — return density relative to the tile median, median intensity, share of single returns, and elevation range — combine them into a water score, keep large low-score-free regions adjacent to ground, polygonize them, and set class 9 on ground and unclassified points inside with filters.overlay and filters.assign.

# Context and Motivation

This guide is part of Water and Bridge Classification in LiDAR. The parent workflow uses return density alone, which works for large calm lakes but struggles with small ponds, rivers flown at low sun angle, and wind-roughened water that scatters more returns. Combining several weak cues produces a far more reliable detector, because each false positive fails a different test: dark asphalt has low intensity but normal density; shadow behind a building has low density but is not flat; a flat roof is flat and single-return but is elevated.

This matters most where no reliable hydrography exists — new reservoirs, gravel pits, seasonal ponds, or projects in regions where national water layers are coarse or out of date.

Four weak cues, one strong detector A table-like diagram with four cue columns: low return density, low intensity, single returns, and flat surface. Rows show water passing all four, dark asphalt failing density, building shadow failing flatness, and a flat roof failing because it is elevated above the local ground. Pass and fail marks fill the cells. low densitylow intensitysingle returnsflat, at ground open waterdark asphaltbuilding shadowflat roof yesyesyesyes noyesyesyes yesyesmixedno noyesyeselevated

# Prerequisites and Assumptions

  • A tile with ground classified and noise removed; intensity populated and roughly consistent across flightlines (see normalizing intensity across flightlines).
  • Python with NumPy, SciPy, rasterio, GeoPandas and PDAL bindings.
  • A minimum water body size from the project specification.

# Step-by-Step Implementation

# Step 1 — Grid the cues

For every 2 m cell, compute the return count, median intensity, share of single returns (NumberOfReturns == 1), and the elevation range of ground-or-unclassified points.

# Step 2 — Score each cell

Convert each cue to a 0–1 score and average them. Empty cells — the strongest water evidence — get a density score of 1 and neutral scores for the cues they cannot measure.

# Step 3 — Threshold and clean

Keep cells with a score above 0.7, then open to remove specks and close to fill scattered surface returns.

# Step 4 — Keep only regions at ground level

Require that each region’s boundary cells sit within about 1 m of the ground surface around it, which rejects flat roofs.

# Step 5 — Burn into points

Polygonize the regions, drop those below the minimum area, and assign class 9 to ground and unclassified points inside.

# Complete Working Example

python
"""Multi-cue water detection on a 2 m grid, burned back into points as class 9."""
from __future__ import annotations

import json
from pathlib import Path

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

CELL = 2.0


def grid_cues(a: np.ndarray) -> tuple[dict[str, np.ndarray], tuple[float, float, int, int]]:
    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)
    idx = r * cols + c
    n = np.bincount(idx, minlength=rows * cols)
    single = np.bincount(idx, weights=(a["NumberOfReturns"] == 1), minlength=rows * cols)
    inten = np.bincount(idx, weights=a["Intensity"], minlength=rows * cols)
    low = np.isin(a["Classification"], [1, 2])
    zmin = np.full(rows * cols, np.inf); zmax = np.full(rows * cols, -np.inf)
    np.minimum.at(zmin, idx[low], a["Z"][low]); np.maximum.at(zmax, idx[low], a["Z"][low])
    with np.errstate(invalid="ignore", divide="ignore"):
        cues = {"n": n, "single": single / n, "intensity": inten / n, "zrange": zmax - zmin,
                "zmin": zmin}
    return {k: v.reshape(rows, cols) for k, v in cues.items()}, (x0, y1, rows, cols)


def water_score(cues: dict[str, np.ndarray]) -> np.ndarray:
    med_n = np.median(cues["n"][cues["n"] > 0])
    med_i = np.nanmedian(cues["intensity"])
    s_density = np.clip(1.0 - cues["n"] / (0.3 * med_n), 0, 1)
    s_int = np.where(np.isnan(cues["intensity"]), 0.5, np.clip(1 - cues["intensity"] / med_i, 0, 1))
    s_single = np.where(np.isnan(cues["single"]), 0.5, cues["single"])
    s_flat = np.where(np.isfinite(cues["zrange"]), np.clip(1 - cues["zrange"] / 0.5, 0, 1), 0.5)
    return (2 * s_density + s_int + s_single + s_flat) / 5.0     # density weighted double


def water_polygons(src: Path, crs: str, min_area: float = 8000.0) -> gpd.GeoDataFrame:
    p = pdal.Pipeline(json.dumps({"pipeline": [str(src)]}))
    p.execute()
    cues, (x0, y1, rows, cols) = grid_cues(p.arrays[0])
    mask = water_score(cues) > 0.7
    mask = ndi.binary_closing(ndi.binary_opening(mask, iterations=1), iterations=3)
    labels, n = ndi.label(mask)
    ground = np.where(np.isfinite(cues["zmin"]), cues["zmin"], np.nan)
    keep = np.zeros_like(mask)
    for lab in range(1, n + 1):
        region = labels == lab
        ring = ndi.binary_dilation(region, iterations=3) & ~region
        inner_edge = region & ~ndi.binary_erosion(region)
        if np.nanmedian(ground[ring]) - np.nanmedian(ground[inner_edge]) > -1.0:
            keep |= region                                  # sits at ground level
    transform = from_origin(x0, y1, CELL, CELL)
    polys = [shape(g) for g, v in rasterio.features.shapes(keep.astype("uint8"), mask=keep,
                                                           transform=transform)]
    gdf = gpd.GeoDataFrame({"geometry": polys}, crs=crs)
    gdf = gdf[gdf.area >= min_area].reset_index(drop=True)
    gdf["WaterId"] = np.arange(1, len(gdf) + 1)
    return gdf


def burn(src: Path, dst: Path, water: gpd.GeoDataFrame) -> None:
    gpkg = dst.with_suffix(".water.gpkg")
    water.to_file(gpkg, layer="water", driver="GPKG")
    pdal.Pipeline(json.dumps({"pipeline": [
        str(src),
        {"type": "filters.ferry", "dimensions": "=>WaterId"},
        {"type": "filters.overlay", "dimension": "WaterId", "datasource": str(gpkg),
         "layer": "water", "column": "WaterId"},
        {"type": "filters.assign", "value": [
            "Classification = 9 WHERE WaterId > 0 && (Classification == 1 || Classification == 2)"]},
        {"type": "writers.las", "filename": str(dst), "minor_version": 4,
         "dataformat_id": 6, "forward": "all"},
    ]})).execute()


if __name__ == "__main__":
    src = Path("valley_0310.laz")
    water = water_polygons(src, "EPSG:6341")
    print(f"{len(water)} water bodies, {water.area.sum() / 1e4:.1f} ha")
    burn(src, Path("valley_0310_water.laz"), water)

# Key Parameter Table

Parameter Type Default Guidance
CELL float, m 2.0 1 m for narrow channels at high density; 3–5 m for sparse data
density reference fraction 0.3 × median Cells above 30 % of normal density score zero for this cue
flatness range float, m 0.5 Elevation range at which a cell stops counting as flat
density weight int 2 The most reliable cue counts double
score threshold float 0.7 Lower finds more small ponds and more false positives
ground-level tolerance float, m 1.0 Region edge must sit within this of surrounding ground
Where the score separates water from land Two overlapping distributions of per-cell water score. Land cells cluster between 0.1 and 0.5. Water cells cluster between 0.75 and 0.95. A small overlap between 0.55 and 0.75 contains shadows and dark roofs. The threshold at 0.7 sits near the right edge of the overlap. 0.7 land cells water cells 0 1 per-cell water score

# Verification

  • Against imagery. Overlay polygons on an orthophoto; every large detected region should be visibly water.
  • Against hydrography. Where national or client water layers exist, compute the share of their area covered by detections, and the share of detected area outside them.
  • No elevated water. Class 9 points should have near-zero height above surrounding ground. Any class 9 point more than a metre above its neighbours’ ground indicates a roof region slipped through.
python
w = gpd.read_file("valley_0310_water.water.gpkg", layer="water")
assert (w.area >= 8000).all()
print(w.area.describe())

# Gotchas and Edge Cases

Specular returns at nadir. Directly below the aircraft, calm water can reflect strongly, producing a stripe of dense, bright returns down the middle of a lake. The closing step fills it when the stripe is narrow; for wide stripes, use per-flightline processing or accept the stripe as water because it is surrounded by water.

Uncalibrated intensity. If intensity differs strongly between flightlines, the intensity cue becomes a flightline detector. Normalize first or drop the cue.

Wetlands and saturated ground. Marshes return weakly and are flat; they will score as water. Whether that is right depends on the specification — decide explicitly and document it.

The nadir stripe and how closing handles it Left: a lake mask with a narrow vertical stripe of non-water cells down the middle, where specular returns made cells look like land. Right: after closing with three iterations, the stripe is filled and the lake is one region. A note says stripes wider than the closing size need per-flightline handling. before closing after closing bright stripe splits the lake one region

Tile edges. A lake cut by a tile edge may fall below the minimum area on each side. Detect on a mosaic of grids, or merge polygons across tiles before applying the area filter.

# Frequently Asked Questions

Why not rely on intensity alone to find water?

Intensity depends on range, incidence angle and sensor calibration, and many dry surfaces such as fresh asphalt and dark roofing are just as dark. It helps as one cue among several but produces too many false positives alone.

What does it mean when water has dense bright returns?

Near nadir, calm water reflects the pulse straight back like a mirror, producing strong returns in a narrow stripe under the flight path. Everywhere else the pulse is reflected away or absorbed, so returns are sparse.

Which points should be set to class 9?

Ground and unclassified returns inside the water polygon. Returns from boats, docks and overhanging vegetation should keep their classes. Specifications differ, so confirm the rule for your project.

Can this detect narrow rivers?

Rivers wider than a few cells at the chosen resolution, if they are open to the sky. Tree-lined streams retain canopy returns and rarely look like water from above; they need breaklines or centreline data.

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