Python LiDAR

Extracting Individual Tree Heights from a CHM

TL;DR: Smooth the CHM slightly, find local maxima with a window scaled to expected crown radius, and treat every result as a candidate until it survives a minimum-height and minimum-separation test — then report the distribution, because individual heights carry a systematic negative bias no algorithm removes.

# Context and Motivation

This guide is part of Canopy Height Models with filters.hag_nn. Having a raster of canopy height invites the obvious next question — how tall is each tree — and the honest answer is that a CHM constrains it rather than answers it.

Two facts set the limits. A laser pulse strikes the upper crown somewhere, rarely the apex, so measured heights run one to two metres low, more in narrow-crowned conifers and at low pulse density. And a local maximum in a raster is a bump, not a tree: a large crown produces several, two adjacent small crowns produce one, and a branch tip near a gap edge produces one that belongs to nothing. Detection is therefore a search with parameters, and the parameters encode what you already believe about the stand.

Three ways a local maximum is not a tree A canopy height profile with four crowns. A fixed detection window finds the apex of the isolated crown correctly, splits one broad crown into two detections, merges two adjacent small crowns into one, and picks a branch tip at a gap edge that belongs to no crown at all. Only one of the four detections is unambiguously right. correct one crown, two detections two crowns, one detection a branch tip the detection window is a statement about crown size — set it from the stand, and expect it to be wrong somewhere in every stand that contains more than one species or age class.

# Prerequisites and Assumptions

Requirement Detail
A CHM from rasterizing a canopy height model
scipy and rasterio for the filtering and raster I/O
Expected crown radius from the species and age of the stand; the single most important input
Field plots if the output will be quoted as absolute heights rather than compared

# Step-by-Step Implementation

# Step 1 — Smooth, but only just

A CHM at one metre carries single-cell spikes from individual branch returns. A small Gaussian removes them without moving the apexes.

python
smoothed = scipy.ndimage.gaussian_filter(chm, sigma=0.7)

# Step 2 — Find local maxima with a crown-scaled window

python
peaks = smoothed == scipy.ndimage.maximum_filter(smoothed, size=window_cells)

The window should be roughly the diameter of the smallest crown you intend to resolve. Too small splits crowns; too large merges them.

# Step 3 — Reject candidates below a minimum height

Anything under two metres is not a tree for most purposes, and rejecting it removes the majority of spurious detections at once.

# Step 4 — Enforce a minimum separation

Where two candidates sit within one crown radius, keep the taller. This is the cheap approximation to crown delineation and gets most of the benefit.

# Step 5 — Report the distribution, not a table of trees

Stem density, height percentiles and the height histogram are defensible. A list of individual trees with heights to two decimal places is not, unless it has been validated against plots.

# Complete Working Example

python
"""Detect candidate tree tops in a CHM and report the stand distribution."""
from __future__ import annotations

import json
from pathlib import Path

import numpy as np
import rasterio
from scipy import ndimage


def detect(chm_path: Path, crown_radius_m: float = 2.5,
           min_height_m: float = 2.0, sigma_cells: float = 0.7) -> dict:
    with rasterio.open(chm_path) as src:
        chm = src.read(1).astype("float32")
        transform = src.transform
        nodata = src.nodata
        cell = abs(transform.a)

    valid = np.isfinite(chm)
    if nodata is not None:
        valid &= chm != nodata
    work = np.where(valid, chm, 0.0)

    smoothed = ndimage.gaussian_filter(work, sigma=sigma_cells)

    window = max(3, int(round(2 * crown_radius_m / cell)) | 1)  # odd, >= 3
    maxima = smoothed == ndimage.maximum_filter(smoothed, size=window)
    candidates = maxima & valid & (work >= min_height_m)

    rows, cols = np.nonzero(candidates)
    heights = work[rows, cols]

    # Enforce minimum separation by keeping the tallest in each labelled clump.
    keep = np.ones(len(rows), dtype=bool)
    order = np.argsort(-heights)
    taken: list[tuple[int, int]] = []
    min_sep_cells = crown_radius_m / cell
    for idx in order:
        r, c = rows[idx], cols[idx]
        if any((r - tr) ** 2 + (c - tc) ** 2 < min_sep_cells ** 2 for tr, tc in taken):
            keep[idx] = False
        else:
            taken.append((r, c))

    kept = heights[keep]
    area_ha = valid.sum() * cell * cell / 10_000.0
    return {
        "cell_size_m": cell,
        "detection_window_cells": window,
        "candidates": int(len(heights)),
        "after_separation": int(len(kept)),
        "stems_per_ha": round(float(len(kept) / max(area_ha, 1e-9)), 1),
        "height_p50": round(float(np.percentile(kept, 50)), 2) if len(kept) else None,
        "height_p95": round(float(np.percentile(kept, 95)), 2) if len(kept) else None,
        "height_max": round(float(kept.max()), 2) if len(kept) else None,
    }


if __name__ == "__main__":
    print(json.dumps(detect(Path("chm.tif"), crown_radius_m=2.5), indent=2))
Four stages, and how many candidates survive each Detection as a funnel. The raw local-maximum search on a smoothed CHM returns 41,200 candidates over the block. Rejecting anything under two metres removes 22,800 of them. Enforcing a minimum separation of one crown radius removes another 12,500. What remains is 5,900 candidate trees, which is 780 per hectare. local maxima 41,200 candidates after min height 2 m 18,400 remain after min separation 5,900 remain per hectare 780 stems/ha 7.6 ha block, 1 m CHM, crown radius 2.5 m most of what a raw maximum filter returns is not a tree, which is why the two rejection stages are the method

# Key Parameter Table

Parameter Typical Effect
crown_radius_m 1.5–5 Sets the detection window; the dominant parameter
min_height_m 2.0 Rejects understory and ground noise
sigma_cells 0.5–1.0 Removes single-cell spikes; larger flattens real apexes
CHM cell size 0.5–1.0 m Finer than crown radius, or crowns cannot be separated
minimum separation ≈ crown radius Cheap stand-in for crown delineation

# Verification

Stem density is plausible. A managed conifer plantation runs 800–2,000 stems per hectare; open oak woodland 50–200. An answer of 12,000 means the window is too small.

The height distribution matches the stand. An even-aged plantation should produce a narrow mode. A broad, flat distribution in a plantation means the detection is picking branches.

Sensitivity is bounded. Re-run with the crown radius halved and doubled. If stem density moves by an order of magnitude, the number is a parameter choice rather than a measurement — say so when reporting it.

How much of the answer is the parameter Detected stems per hectare plotted against the assumed crown radius, for the same CHM. At one metre the detector reports 4,100 stems per hectare; at two metres 1,450; at three metres 780; at five metres 310. The field-measured value of about 900 is marked, and it corresponds to a crown radius of roughly 2.7 metres. field plots: about 900 stems/ha 4,100 at r = 1 m 1,450 at r = 2 m 780 at r = 3 m 0 2,000 4,200 1 m 3 m 5 m assumed crown radius

# Gotchas and Edge Cases

Heights are biased low and the bias is not constant. It grows with narrower crowns and lower pulse density, so it does not cancel between a dense and a sparse acquisition.

A CHM coarser than the crowns cannot separate them. At a two-metre cell in a stand with three-metre crowns, detection is measuring the raster.

Buildings and poles detect beautifully. Filter to vegetation classes before rasterizing, as in rasterizing a canopy height model.

Validate before quoting absolutes. Relative comparisons across one block are robust; absolute stem counts and heights need plots.

# Frequently Asked Questions

Why are LiDAR tree heights lower than field measurements?

Because a pulse strikes the upper crown somewhere, rarely the apex. The resulting bias is typically one to two metres and grows with narrower crowns and lower pulse density, so it does not cancel between a dense and a sparse acquisition. Relative comparisons within one block are robust; absolute heights need local calibration.

How do I choose the detection window?

From the stand, not from the raster. The window should be roughly the diameter of the smallest crown you intend to resolve — too small and one broad crown becomes several detections, too large and two adjacent crowns become one. Expect it to be wrong somewhere in any mixed-species or mixed-age stand.

Is a local maximum a tree?

No, it is a bump. A large crown can produce several, two small adjacent crowns can produce one, and a branch tip at a gap edge produces one belonging to nothing. Treating detections as candidates that must survive a height and separation test is what makes the output defensible.

How do I know whether my stem count is real?

Run the sensitivity test. If halving and doubling the crown radius moves stem density by an order of magnitude — which it usually does — the number is a parameter choice, and it should be reported as one unless field plots have pinned the parameter down.

Related in Canopy Height Models with filters.hag_nn

↑ Canopy Height Models with filters.hag_nn overview --- title: "Canopy Height Models with filters.hagnn" description: "Producing a canopy height model in the point domain: classify ground,... Computing Height Above Ground with filters.hag_nn --- title: "Computing Height Above Ground with filters.hagnn" description: "What filters.hagnn actually computes, why count and maxdistance... Rasterizing a Canopy Height Model from HAG --- title: "Rasterizing a Canopy Height Model from HAG" description: "The three writers.gdal choices that decide whether you get a canopy...