Individual Tree Detection and Crown Segmentation from LiDAR in Python

A stand-level canopy height raster tells you how tall the forest is; an inventory tells you how many trees there are, how big each one is, and how they are arranged. Bridging that gap is the job of individual tree detection (ITD) and crown segmentation: taking a continuous height surface or a normalized point cloud and resolving it into discrete tree objects with a stem location, a height, and a crown boundary. The concrete problem this page solves: you have a validated 0.5 m canopy height model over a mixed-conifer or broadleaf stand, and you need a point layer of tree tops plus a polygon layer of crowns whose counts hold up against a field plot tally. This work sits inside the wider Canopy Height Modeling & Terrain Extraction pipeline and consumes the surface produced by Canopy Height Model Creation directly.

Tree detection is where forestry LiDAR stops being remote sensing and starts being inventory. Every downstream number — stems per hectare, diameter distribution, basal area, individual-tree biomass, the crown metrics that feed Fuel Load Mapping from LiDAR — is an aggregate over the objects this step produces. Under-segment and you merge three suppressed stems into one dominant crown; over-segment and a single wide oak becomes four trees. Neither error is visible in the raster, and both propagate silently.

Prerequisites

Every item below has to be true before detection numbers mean anything — a CHM with pits or a stale terrain baseline produces confident, wrong tree counts:

• A validated canopy height model, pit-filled and at 0.25–1.0 m resolution, produced as described in Canopy Height Model Creation.
• Python ≥ 3.10 with rasterio ≥ 1.3, numpy ≥ 1.24, scipy ≥ 1.10, scikit-image ≥ 0.22, geopandas ≥ 0.14, and shapely ≥ 2.0.
• A projected, metre-based CRS on both the raster and any field plot layer — a search window expressed in metres is meaningless in degrees.
• Pulse density of at least 8 returns/m² if you intend to segment in the point cloud rather than the raster; below that, raster methods are more stable.
• A field plot layer with stem coordinates and measured heights for at least a few plots, in the same CRS, for validation.
• A minimum-height threshold decided in advance (commonly 2 m) so shrubs and slash are excluded from the candidate set.

If the height surface is still noisy, resolve that first — pit artifacts create spurious local maxima that no detection parameter can filter out. The dedicated fix is documented under Removing CHM Pits with the Spike-Free Algorithm.

Concept: Why a Fixed Search Window Fails

The workhorse detector is the local maximum filter: a pixel is a candidate tree top if it is the highest value within a circular window centred on it. The window radius is the entire difficulty. Too small and every branch tip in a wide crown becomes its own “tree”; too large and neighbouring suppressed stems disappear into their dominant neighbour.

Crown radius scales with tree height, so the window must scale too. The standard remedy is a variable window filter: the search radius is a function of the pixel’s own height value. A linear form covers most temperate stands,

$$r(h)=a+b\cdot h$$

where $h$ is canopy height in metres and $r$ the search radius in metres. Typical conifer coefficients are $a\approx 1.2$, $b\approx 0.06$; broadleaf crowns are wider for their height, so $b$ climbs toward 0.10–0.14. Where the height–crown relationship is distinctly non-linear — open-grown oaks, for instance — a saturating form fits better:

$$r(h)=r_{\min }+(r_{\max }-r_{\min })\left(1-e^{-kh}\right)$$

The two coefficients are not universal constants; they are a local allometry, and the honest way to obtain them is to regress measured crown radius on measured height from your own plots. Borrowing published coefficients from a different forest type is the single most common cause of a systematic count bias.

The second half of the problem is extent. A detected top is a point; a crown is an area. Marker-controlled watershed solves it by treating the negated height surface as a topographic basin map: each detected top is a marker, water floods outward from every marker simultaneously, and the ridge where two floods meet becomes the crown boundary. Because the markers are supplied rather than discovered, the number of crowns is fixed by detection, and segmentation only decides where the edges fall.

Step-by-Step Python Pipeline

The flow is: read and mask → smooth → detect tops with a variable window → segment crowns → attribute and export. Each step below is independently runnable.

Step 1 — Read the CHM and mask non-forest

Load the height surface, replace NoData with zero, and apply the minimum-height threshold that keeps shrubs out of the candidate pool.

import numpy as np
import rasterio


def load_chm(chm_path: str, min_height: float = 2.0):
    """Return (height array with NoData as 0, transform, crs, forest mask)."""
    with rasterio.open(chm_path) as src:
        chm = src.read(1, masked=True).astype("float32")
        transform = src.transform
        crs = src.crs
        res = abs(src.transform.a)
    heights = chm.filled(0.0)
    heights[heights < 0] = 0.0
    forest = heights >= min_height
    return heights, transform, crs, res, forest

Filling NoData with zero rather than NaN matters: the maximum filter in the next step propagates NaN across the whole window, which would erase every tree adjacent to a data gap.

Step 2 — Smooth just enough to remove branch-tip noise

A light Gaussian blur suppresses the single-pixel spikes that individual branch tips produce without moving the true apex. Keep $\sigma$ below one pixel — heavier smoothing shifts tops downslope in leaning crowns.

from scipy.ndimage import gaussian_filter


def smooth_chm(heights: np.ndarray, sigma: float = 0.6) -> np.ndarray:
    """Light Gaussian smoothing; sigma is in pixels, deliberately sub-pixel-ish."""
    return gaussian_filter(heights, sigma=sigma, mode="nearest")

Step 3 — Detect tree tops with a variable window

Rather than running one maximum filter, run a small bank of them at the discrete radii the allometry implies, then keep each pixel’s result from the bank entry that matches its own height. This is vectorised and far faster than a per-pixel loop.

from scipy.ndimage import maximum_filter


def variable_window_maxima(heights, res, a=1.2, b=0.09, min_height=2.0):
    """Boolean array of tree tops using radius r(h) = a + b*h metres."""
    radii_m = a + b * heights
    radii_px = np.clip(np.round(radii_m / res).astype("int16"), 1, 25)

    tops = np.zeros(heights.shape, dtype=bool)
    for k in np.unique(radii_px):
        size = int(2 * k + 1)
        local_max = maximum_filter(heights, size=size, mode="nearest")
        # A pixel is a top if it equals the local max computed at *its own* radius.
        band = (radii_px == k) & (heights == local_max) & (heights >= min_height)
        tops |= band
    return tops

The heights == local_max test is exact equality on floats, which is correct here because the maximum filter returns a value copied from the array rather than a computed one. Two adjacent pixels that share the same value both survive; step 4 collapses such ties into one marker.

Step 4 — Turn tops into labelled markers

Adjacent top pixels on a plateau are one tree. Label connected groups and take each group’s centroid.

from scipy import ndimage as ndi


def tops_to_markers(tops: np.ndarray):
    """Label connected top pixels; return (marker raster, list of (row, col))."""
    markers, n = ndi.label(tops, structure=np.ones((3, 3), dtype=int))
    if n == 0:
        return markers, []
    centres = ndi.center_of_mass(tops, markers, range(1, n + 1))
    return markers, [(int(round(r)), int(round(c))) for r, c in centres]

Step 5 — Segment crowns with a marker-controlled watershed

Flood the negated height surface from the markers, constrained to the forest mask so crowns cannot leak across clearings.

from skimage.segmentation import watershed


def segment_crowns(heights: np.ndarray, markers: np.ndarray, forest: np.ndarray):
    """Marker-controlled watershed on the inverted CHM; returns a crown label raster."""
    labels = watershed(-heights, markers=markers, mask=forest, watershed_line=False)
    return labels.astype("int32")

Step 6 — Vectorise crowns and attribute each tree

Polygonise the label raster, attach the top height to each crown, and drop implausibly small or large polygons.

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


def crowns_to_gdf(labels, heights, transform, crs, min_area=1.5, max_area=400.0):
    """Vectorise crown labels into a GeoDataFrame with height and area attributes."""
    records = []
    for geom, value in shapes(labels, mask=labels > 0, transform=transform):
        poly = shape(geom)
        if not (min_area <= poly.area <= max_area):
            continue
        mask = labels == int(value)
        records.append({
            "tree_id": int(value),
            "height_m": float(heights[mask].max()),
            "crown_area_m2": float(poly.area),
            "crown_diam_m": float(2.0 * (poly.area / np.pi) ** 0.5),
            "geometry": poly,
        })
    return gpd.GeoDataFrame(records, crs=crs)


def tops_to_gdf(centres, heights, transform, crs):
    """Tree-top points in map coordinates with their height."""
    rows = []
    for r, c in centres:
        x, y = rasterio.transform.xy(transform, r, c)
        rows.append({"height_m": float(heights[r, c]), "geometry": Point(x, y)})
    return gpd.GeoDataFrame(rows, crs=crs)

The crown_diam_m column is an area-equivalent circular diameter, which is what most allometric equations expect. Keep the raw area as well: crown shape irregularity is itself a useful species discriminator.

Parameter Reference

Validation & Verification

Detection accuracy is a matching problem, not a counting problem. Two runs can both report 412 trees on a plot that holds 400 and be very differently wrong. Link each detection to at most one field stem within a horizontal tolerance, then report precision, recall, and their harmonic mean:

$$P=\frac{TP}{TP+FP},R=\frac{TP}{TP+FN},F_1=\frac{2PR}{P+R}$$

A greedy nearest-neighbour match, taken in descending height order so dominant trees claim their partner first, is the standard and reproducible choice.

import numpy as np
import geopandas as gpd


def match_detections(detected: gpd.GeoDataFrame, field: gpd.GeoDataFrame,
                     max_dist: float = 3.0, max_dh: float = 5.0) -> dict:
    """Greedy height-ordered matching of detected tops to field stems."""
    det = detected.sort_values("height_m", ascending=False).reset_index(drop=True)
    used = np.zeros(len(field), dtype=bool)
    tp = 0
    for _, d in det.iterrows():
        dist = field.geometry.distance(d.geometry).to_numpy()
        dh = np.abs(field["height_m"].to_numpy() - d["height_m"])
        ok = (dist <= max_dist) & (dh <= max_dh) & (~used)
        if not ok.any():
            continue
        j = int(np.where(ok, dist, np.inf).argmin())
        used[j] = True
        tp += 1
    fp = len(det) - tp
    fn = int((~used).sum())
    precision = tp / (tp + fp) if tp + fp else 0.0
    recall = tp / (tp + fn) if tp + fn else 0.0
    f1 = 2 * precision * recall / (precision + recall) if precision + recall else 0.0
    return {"tp": tp, "fp": fp, "fn": fn,
            "precision": precision, "recall": recall, "f1": f1}

Read the two error types as separate diagnoses rather than one score. Low precision with high recall means the window is too small — branch tips are being promoted to trees, so raise b. High precision with low recall means the window is too large or the stand is multi-layered: suppressed stems under a closed overstory are simply not visible from above, and no raster parameter recovers them. Published F₁ for raster methods runs roughly 0.85–0.95 in even-aged conifer plantations and falls to 0.55–0.75 in multi-layered broadleaf stands; a score outside that envelope is usually a matching-tolerance artifact rather than a genuine result.

The other half of validation is the height itself. Regress detected height on field height and check that the slope is near one and the intercept near zero:

def height_bias(matched_pairs: np.ndarray) -> dict:
    """matched_pairs: (n, 2) array of [detected_h, field_h]."""
    d, f = matched_pairs[:, 0], matched_pairs[:, 1]
    resid = d - f
    slope, intercept = np.polyfit(f, d, 1)
    return {
        "n": int(len(d)),
        "bias_m": float(resid.mean()),
        "rmse_m": float(np.sqrt((resid ** 2).mean())),
        "slope": float(slope),
        "intercept_m": float(intercept),
    }

A negative bias of 0.3–1.0 m is expected and physical: the laser rarely strikes the exact apex, so LiDAR heights under-estimate slightly. A bias beyond about 1.5 m points at the terrain baseline instead, and belongs back in Digital Terrain Model Generation.

Failure Modes & Gotchas

• Coefficients borrowed from another forest type. A conifer window applied to open-grown broadleaves over-segments every large crown. Fit a and b against your own plot crown radii before trusting a regional number.
• NoData filled with NaN. maximum_filter propagates NaN across the whole kernel, silently deleting every tree within one window of a data gap. Fill with zero and mask afterwards.
• Detecting on an unsmoothed CHM. Raw surfaces carry single-pixel spikes from branch tips and from pit-filling artifacts; both are perfect local maxima and neither is a tree.
• Watershed without a mask. An unmasked flood assigns clearings, roads, and water to the nearest crown, inflating every crown area near a stand edge. Always pass the forest mask.
• Tile-edge duplication. A tree whose crown straddles two tiles is detected twice, once per tile. Detect on buffered tiles and drop detections that fall outside the nominal extent — the same halo discipline used in Tiling and Buffering Large LiDAR Surveys with PDAL.
• Comparing counts against a plot with a different minimum diameter. Field crews usually tally stems above a diameter threshold; the raster counts crowns visible from above. Reconcile the two definitions before declaring a bias.

Performance & Scale Notes

The maximum-filter bank is the expensive part: cost grows with the number of distinct radii, not with the number of trees. Quantise radii_px to a handful of classes — rounding to the nearest 2 px typically leaves fewer than eight filter passes with no measurable accuracy loss. On a 2000 × 2000 tile at 0.5 m the whole pipeline runs in a few seconds per tile, so a regional survey is a matter of parallel tiles rather than algorithmic cleverness.

import concurrent.futures
from pathlib import Path


def detect_tile(chm_path: str, out_dir: str) -> str:
    """Full detection + segmentation for one CHM tile; returns the crown file path."""
    heights, transform, crs, res, forest = load_chm(chm_path)
    smooth = smooth_chm(heights)
    tops = variable_window_maxima(smooth, res)
    markers, centres = tops_to_markers(tops)
    labels = segment_crowns(smooth, markers, forest)
    crowns = crowns_to_gdf(labels, heights, transform, crs)
    out = Path(out_dir) / (Path(chm_path).stem + "_crowns.gpkg")
    crowns.to_file(out, driver="GPKG")
    return str(out)


def detect_survey(chm_tiles: list[str], out_dir: str, workers: int = 6) -> list[str]:
    Path(out_dir).mkdir(parents=True, exist_ok=True)
    with concurrent.futures.ProcessPoolExecutor(max_workers=workers) as pool:
        return list(pool.map(detect_tile, chm_tiles, [out_dir] * len(chm_tiles)))

Segment on a buffered tile — 20 m is ample for temperate crowns — and clip the resulting crown polygons back to the nominal tile before mosaicking, keeping only crowns whose top falls inside the nominal extent. That rule assigns every straddling tree to exactly one tile without duplication or loss.

Frequently Asked Questions

Related

• Canopy Height Model Creation — produces the height surface this page consumes.
• Local Maxima Tree Detection with a Variable Window Filter — the detection step in isolation, with coefficient fitting.
• Watershed Crown Segmentation from a CHM in Python — the segmentation step in isolation.
• Point Cloud Tree Segmentation with the Li 2012 Algorithm — the 3-D alternative for high-density surveys.
• Forest Gap & Understory Analysis — the complement: what happens between the crowns.

Up: Canopy Height Modeling & Terrain Extraction