Forestry GIS Python · Forestry · Ecology

Watershed Crown Segmentation from a CHM in Python

Detection gives you a stem location; segmentation gives that stem an extent. This guide covers the delineation half only: you already have tree-top markers from Local Maxima Tree Detection with a Variable Window Filter, and you need crown polygons whose areas are usable in an allometric equation. The wider workflow is Individual Tree Detection and Crown Segmentation, inside Canopy Height Modeling & Terrain Extraction.

When to use this approach

Marker-controlled watershed is the default because it guarantees one crown per marker: the segment count is decided by detection and cannot drift. Its weakness is boundary placement in stands where crowns interlock without a height dip between them.

Method Segment count Boundary quality in interlocking crowns Cost When it wins
Marker-controlled watershed Fixed by markers Fair — the ridge follows the height valley, which may not exist Low Almost always; the reliable default
Compact watershed Fixed by markers Better — a compactness term keeps crowns from sprawling Low Wide crowns bleeding into their neighbours
Region growing from seeds Fixed by seeds Good with a well-chosen height-drop rule Moderate Conifers with a strong conical profile
Unseeded watershed Discovered Poor — every local dip becomes a crown Low Never, for inventory work

Use the plain watershed first. Reach for the compact variant only when validation shows crown areas that are implausibly large at stand edges, which is the signature of a flood that had nowhere to stop.

Minimal reproducible example

The whole segmentation is one call plus the bookkeeping that makes it usable. The inputs are the smoothed height array, an integer marker raster in which each tree has its own label, and a boolean forest mask.

import numpy as np
import rasterio
from scipy import ndimage as ndi
from skimage.segmentation import watershed


def segment_crowns(chm_path: str, top_rows, top_cols, min_height: float = 2.0):
    """Marker-controlled watershed over a CHM. Returns (labels, heights, profile)."""
    with rasterio.open(chm_path) as src:
        chm = src.read(1, masked=True).astype("float32")
        profile = src.profile

    heights = chm.filled(0.0)
    heights[heights < 0] = 0.0
    forest = heights >= min_height

    markers = np.zeros(heights.shape, dtype="int32")
    markers[np.asarray(top_rows), np.asarray(top_cols)] = np.arange(
        1, len(top_rows) + 1, dtype="int32")

    labels = watershed(-heights, markers=markers, mask=forest)
    return labels.astype("int32"), heights, profile

Two details carry most of the risk. The surface is negated because a watershed floods basins and a crown is a peak. The mask is not optional: without it the flood assigns every clearing, road, and river to whichever crown happens to be nearest, and crown areas at stand edges become meaningless.

Watershed segmentation with and without the forest mask Two plan-view panels of the same stand. On the left, without a mask, the crown regions expand outward until they meet, so a clearing and a road strip are absorbed into the neighbouring crowns and their areas are inflated. On the right, with a two-metre forest mask applied, each crown region stops at the mask boundary and the clearing and road remain unlabelled background. No mask — crowns absorb the clearing Masked at 2 m — crowns stop at forest clearing, still labelled road, still labelled unassigned unassigned The mask is what keeps crown area an ecological quantity rather than a Voronoi artefact.

Parameter reference

Parameter Type Default Recommended range Ecological rationale
mask boolean array forest ≥ 2 m — Confines the flood to canopy; without it crown area is a distance artefact
connectivity int 1 1–2 1 uses 4-neighbours and produces blockier but more conservative crowns
compactness float 0 0–0.01 Penalises sprawl; raise it only when edge crowns come back implausibly large
watershed_line bool False — True leaves a one-pixel gap between crowns, which loses about 5% of crown area
min_area float (m²) 1.5 1.0–4.0 Drops slivers created where three basins meet
max_area float (m²) 400 150–600 Flags merged clumps; set above the largest plausible single crown

Leave watershed_line at False for inventory work. The visual appeal of separated crowns is not worth systematically shrinking every crown by its perimeter.

Clean up and vectorise

Raw watershed output always contains a few slivers where three floods meet, and occasionally one enormous region where a marker was missing. Filter on area, then polygonise.

import geopandas as gpd
from rasterio.features import shapes
from shapely.geometry import shape


def crowns_to_polygons(labels, heights, transform, crs,
                       min_area=1.5, max_area=400.0):
    """Vectorise crown labels, attach metrics, and drop implausible polygons."""
    keep, dropped = [], {"sliver": 0, "clump": 0}
    for geom, value in shapes(labels, mask=labels > 0, transform=transform):
        poly = shape(geom)
        if poly.area < min_area:
            dropped["sliver"] += 1
            continue
        if poly.area > max_area:
            dropped["clump"] += 1
            continue
        sel = labels == int(value)
        keep.append({
            "tree_id": int(value),
            "height_m": float(heights[sel].max()),
            "crown_area_m2": float(poly.area),
            "crown_diam_m": float(2.0 * (poly.area / np.pi) ** 0.5),
            "compactness": float(4 * np.pi * poly.area / (poly.length ** 2)),
            "geometry": poly,
        })
    return gpd.GeoDataFrame(keep, crs=crs), dropped

The compactness column — the isoperimetric ratio, 1.0 for a perfect circle — is worth keeping. Values below about 0.3 mark polygons that snaked around an obstacle, which is nearly always a segmentation error rather than a genuinely ragged crown.

Expected output and verification

A correct segmentation has one polygon per detection (minus the slivers you deliberately dropped), and a crown-area distribution that rises with height in the way local allometry predicts. The strongest single check is that regression, because it needs no field data at all:

import numpy as np


def crown_area_sanity(crowns) -> dict:
    """Crown area should increase with height; a flat or negative slope is a red flag."""
    h = crowns["height_m"].to_numpy()
    r = np.sqrt(crowns["crown_area_m2"].to_numpy() / np.pi)
    slope, intercept = np.polyfit(h, r, 1)
    return {
        "n_crowns": int(len(crowns)),
        "radius_vs_height_slope": float(slope),
        "median_area_m2": float(np.median(crowns["crown_area_m2"])),
        "p95_area_m2": float(np.percentile(crowns["crown_area_m2"], 95)),
        "median_compactness": float(np.median(crowns["compactness"])),
    }


stats = crown_area_sanity(crowns)
assert stats["radius_vs_height_slope"] > 0.02, "crown radius does not grow with height"
assert stats["median_compactness"] > 0.35, "crowns are snaking — check the mask"

A slope near zero means the markers are wrong, not the watershed: if every crown gets roughly the same area regardless of height, the flood is being stopped by marker spacing rather than by the height surface, which happens when the detection window was far too small.

Crown-area distributions from a healthy and a failed segmentation Two histograms side by side. The left one, from a healthy run, is smooth and right-skewed with a mode near twenty square metres and a thin tail. The right one, from a failed run, has a tall spike in the smallest bin representing sliver polygons and a second cluster of bars beyond three hundred square metres representing merged clumps, with very little mass in between. Healthy run Failed run — slivers and clumps 0 30 60 100+ crown area (m²) slivers clumps 0 30 60 100+ crown area (m²) A bimodal area histogram is the fastest visual diagnosis of a bad marker set. What a missing marker does to the crown map Two panels of the same canopy. In the first, both crowns carry a marker and the watershed produces two crowns of about thirty square metres each. In the second, the marker on the right crown is missing, so the flood from the left marker expands across both and produces a single region of about seventy square metres, which exceeds the maximum plausible crown area and is caught by the area filter. two markers one marker missing 31 m² 28 m² boundary at the height valley 69 m² — one region caught by the max-area filter The area filter is a detector for detection failures, not just a cosmetic clean-up.

Common pitfalls

  • Segmenting the unsmoothed CHM while detecting on the smoothed one. The two surfaces have slightly different ridges, so crown boundaries land a pixel or two off every marker. Flood the same array you detected on, then read heights from the raw surface.
  • Markers outside the mask. A detection on a pixel just below the height threshold becomes a label with no basin, and skimage silently returns an empty region for it — the crown count then quietly disagrees with the detection count.
  • Forgetting that label 0 is background. Marker numbering must start at 1; a tree numbered 0 is treated as unlabelled and its crown is absorbed by a neighbour.
  • Filtering slivers before checking how many there are. A handful is normal. Hundreds means the marker set has clusters of duplicates and the fix belongs upstream in detection, not in an area filter.

Frequently Asked Questions

Should I use watershed_line=True so crowns do not touch?

Only for cartography. The line consumes one pixel around every crown, which at 0.5 m resolution removes roughly 5–10% of the area of a typical crown and biases every area-based allometric estimate downward. Keep it off for analysis and dissolve or buffer for display instead.

My crowns look square. Is something wrong?

Probably not. A watershed on a raster produces pixel-edge boundaries, and with connectivity=1 those edges are axis-aligned. If the shapes matter — for a species classifier that uses crown outline, or for a figure — simplify with a tolerance near half a pixel after computing area, never before.

How do I handle crowns that straddle a tile boundary?

Segment on a buffered tile and keep only crowns whose marker falls inside the nominal extent. That assigns each straddling tree to exactly one tile with no duplicates and no losses, and it is the same halo rule used for the classification stages in Tiling and Buffering Large LiDAR Surveys with PDAL.

Up: Individual Tree Detection & Segmentation · Canopy Height Modeling & Terrain Extraction