Python LiDAR

Segmenting Trees with a Watershed on a CHM

TL;DR: Smooth the CHM lightly, find tree tops as local maxima with a search window that grows with height (a 20 m tree gets a wider window than a 6 m one), run skimage.segmentation.watershed on the negated CHM with those tops as markers and a canopy mask, and polygonize the labels into crowns carrying the maximum CHM height as tree height.

# Context and Motivation

This guide is part of Individual Tree Segmentation from LiDAR. A watershed treats the canopy height model as a landscape turned upside down: every tree top becomes the bottom of a basin, and crowns are the areas that drain into each basin. Left alone, a watershed floods from every local minimum and produces hundreds of fragments per hectare. Marker control fixes that by flooding only from the tops you nominate, which turns the problem into choosing good markers.

The refinement worth making over a fixed-window approach is a variable window. Small trees have small crowns and can stand close together; large trees have wide crowns, and two maxima five metres apart on a 30 m tree are usually the same crown. A window whose radius is a function of height captures that without separate runs for each stand type.

Crowns as basins Upper strip: a canopy profile with three peaks. Lower strip: the same profile negated so the peaks become basins. Markers sit at the bottom of each basin. Water rising from each marker meets its neighbour at the ridges, which become crown boundaries. CHM −CHM markers at tops; boundaries where floods meet

# Prerequisites and Assumptions

  • A canopy height model at 0.5 m, ideally pit-free, in a projected CRS in metres.
  • Python with rasterio, NumPy, SciPy, scikit-image 0.19+ and GeoPandas.
  • A minimum tree height of interest (commonly 2–5 m) and a rough crown-width to height relationship for the forest type.

# Step-by-Step Implementation

# Step 1 — Smooth lightly

A Gaussian with sigma of 0.5–1 m suppresses branch-level bumps. Keep the unsmoothed CHM for heights.

# Step 2 — Detect tops with a variable window

For each candidate maximum, require that it is the highest cell within a radius that depends on its own height, for example r = 0.6 + 0.08·h metres. Implement it by testing candidates against a maximum filter per height band.

# Step 3 — Mask the canopy

Build a mask of cells above a fixed height (such as 2 m) or a fraction of each top’s height, so crowns do not grow over gaps and bare ground.

# Step 4 — Run the watershed

Pass the negated smoothed CHM, the marker image and the mask to watershed. Setting compactness above zero produces more rounded crowns, useful in dense conifer stands.

# Step 5 — Polygonize and attribute

Convert labels to polygons, then attach tree height (maximum of the raw CHM inside the crown), crown area and top location.

# Complete Working Example

python
"""Variable-window marker-controlled watershed on a canopy height model."""
from __future__ import annotations

from pathlib import Path

import geopandas as gpd
import numpy as np
import rasterio
import rasterio.features
from scipy import ndimage as ndi
from shapely.geometry import shape
from skimage.segmentation import watershed


def window_radius(h: np.ndarray) -> np.ndarray:
    """Search radius in metres as a function of height (tune per forest type)."""
    return 0.6 + 0.08 * h


def variable_window_tops(smooth: np.ndarray, res: float, min_h: float) -> np.ndarray:
    tops = np.zeros(smooth.shape, dtype=bool)
    bands = np.arange(min_h, np.nanmax(smooth) + 5.0, 5.0)
    for lo in bands:
        hi = lo + 5.0
        r_cells = max(1, int(round(window_radius(np.array(hi)) / res)))
        size = 2 * r_cells + 1
        yy, xx = np.ogrid[-r_cells:r_cells + 1, -r_cells:r_cells + 1]
        disk = (xx ** 2 + yy ** 2) <= r_cells ** 2
        local_max = ndi.maximum_filter(smooth, footprint=disk, mode="nearest")
        in_band = (smooth >= lo) & (smooth < hi)
        tops |= in_band & (smooth == local_max)
    return tops


def segment(chm_path: Path, out_gpkg: Path, min_h: float = 3.0, sigma_m: float = 0.6,
            compactness: float = 0.0) -> gpd.GeoDataFrame:
    with rasterio.open(chm_path) as ds:
        chm = ds.read(1, masked=True).filled(0.0).astype("float32")
        transform, crs, res = ds.transform, ds.crs, ds.res[0]
    chm[chm < 0] = 0.0
    smooth = ndi.gaussian_filter(chm, sigma=sigma_m / res)

    tops = variable_window_tops(smooth, res, min_h)
    markers, n = ndi.label(tops)
    mask = smooth > max(2.0, 0.5 * min_h)
    labels = watershed(-smooth, markers=markers, mask=mask, compactness=compactness)

    ids = np.arange(1, labels.max() + 1)
    heights = ndi.maximum(chm, labels, ids)
    records = []
    for geom, lab in rasterio.features.shapes(labels.astype("int32"), mask=labels > 0,
                                              transform=transform):
        poly = shape(geom)
        h = float(heights[int(lab) - 1])
        if h < min_h or poly.area < 1.0:
            continue
        records.append({"tree_id": int(lab), "height_m": round(h, 2),
                        "crown_area_m2": round(poly.area, 1), "geometry": poly})
    crowns = gpd.GeoDataFrame(records, crs=crs)
    crowns.to_file(out_gpkg, layer="crowns", driver="GPKG")
    print(f"{n} tops, {len(crowns)} crowns written to {out_gpkg.name}")
    return crowns


if __name__ == "__main__":
    segment(Path("out/stand_12/chm.tif"), Path("out/stand_12/crowns.gpkg"))
A window that grows with the tree A straight line of window radius against tree height, from 0.84 metres at 3 metres height to 3.0 metres at 30 metres height. Three example trees are drawn beneath the line with circles of the corresponding radius: a small tree with a tight window, a medium tree, and a tall tree with a wide window. 3 m tree: r 0.8 m 15 m tree: r 1.8 m 30 m tree: r 3.0 m tree height (m) window radius r = 0.6 + 0.08 h — a starting point, tuned against plots

# Key Parameter Table

Parameter Type Default Guidance
sigma_m float, m 0.6 0.4 for young dense stands, 1.0+ for broad crowns
window intercept float, m 0.6 Minimum radius; stops adjacent cells both being tops
window slope float 0.08 Radius growth per metre of height; from a crown-width model if you have one
min_h float, m 3.0 Smallest tree reported
mask height float, m max(2, 0.5·min_h) Crown edge; lower lets crowns spread across gaps
compactness float 0.0 0.001–0.01 rounds crowns in dense conifers

# Verification

  • Tops per hectare. Compare against stand records or plots. A number far above expected density means the window is too small or smoothing too weak.
  • Heights against raw points. For a sample of crowns, the tree height should match the maximum HeightAboveGround of points inside the crown within one CHM cell’s error.
  • Crown shapes. Overlay crowns on a hillshade of the CHM. Boundaries should follow the valleys between crowns; straight boundaries through a crown mean two markers in one tree.
python
crowns = gpd.read_file("out/stand_12/crowns.gpkg", layer="crowns")
ha = crowns.total_bounds
area_ha = (ha[2] - ha[0]) * (ha[3] - ha[1]) / 10_000
print(f"{len(crowns) / area_ha:.0f} trees/ha, median height {crowns.height_m.median():.1f} m")

# Gotchas and Edge Cases

Flat-topped crowns produce plateaus. When several adjacent cells share the exact maximum, the equality test marks all of them. ndi.label merges touching top cells into one marker, which is why the code labels the boolean top image rather than using the cells directly.

Border trees are truncated. Crowns cut by the tile edge are smaller and their tops may be missing. Process with a buffer at least one crown diameter wide and keep only trees whose top lies inside the unbuffered tile.

Deciduous leaf-off data. Without leaves, crowns are open branch networks and the CHM is full of holes. Build the CHM with a larger radius, use a pit-free method, or segment from points instead.

What the canopy mask prevents Two plan views of a forest edge next to a field. Without a mask, the outermost crown's watershed region spreads far into the open field, inflating its area. With a mask at two metres, the crown stops at the canopy edge and its area matches the visible crown. no mask mask at 2 m crown leaks into field open field

# Frequently Asked Questions

Why use a variable window instead of a fixed one?

Crown width grows with tree height, so a single window is either too small for large trees, splitting them, or too large for small ones, merging them. A height-dependent window adapts the separation to each tree.

What does compactness do in skimage watershed?

It adds a penalty for distance from the marker, so regions grow more evenly in all directions. Small values give rounder crowns and help in dense stands where intensity valleys between crowns are weak.

Can I use the smoothed CHM for tree heights?

No. Smoothing lowers peaks, often by tens of centimetres. Use the raw CHM maximum within each crown, or the maximum point height, for reported tree height.

Is watershed segmentation good for broadleaf forests?

It works but less well than in conifers, because broadleaf crowns are wide, flat and interlocking. Expect lower detection rates and tune smoothing more aggressively, or consider point-based methods for dense broadleaf stands.

Related in Individual Tree Segmentation from LiDAR

↑ Individual Tree Segmentation from LiDAR overview --- title: "Individual Tree Segmentation from LiDAR" description: "Split a forest canopy into individual trees: local-maximum tree tops,... Segmenting Trees Directly from Points --- title: "Segmenting Trees Directly from Points" description: "Point-based individual tree segmentation with PDAL filters.litree:... Computing Crown Metrics per Tree --- title: "Computing Crown Metrics per Tree" description: "Turn segmented trees into an inventory table: height, crown base height, crown...