Local Maxima Tree Detection with a Variable Window Filter in Python

The window radius is the only parameter that materially changes a local-maxima tree count, and most workflows inherit it from a paper about a different forest. This guide is about the narrow task of choosing that radius from your own data: fitting a crown-radius-on-height relationship from plot measurements, encoding it as a search window, and sweeping the coefficients until precision and recall balance. It is one step inside Individual Tree Detection and Crown Segmentation, which itself belongs to Canopy Height Modeling & Terrain Extraction.

When to use this approach

A local-maxima filter over a raster is the right detector when the trees you care about have an apex visible from above. It is fast, deterministic, and easy to defend. It cannot see a suppressed stem beneath a closed overstory, and no coefficient choice will make it.

Choose the variable window whenever the stand spans more than one size class, which in practice is nearly always. The fixed radius is defensible only in a uniform, even-aged plantation where every crown really is the same width.

Fit the allometry before you pick a radius

Crown radius scales with height, and the scaling is local. Measure crown radius and height on a sample of stems — two perpendicular crown widths averaged is the field standard — and fit the straight line the filter will use:

$$r(h)=a+bh$$

The intercept $a$ is the radius floor: the crown radius of the smallest tree you intend to detect. The slope $b$ is the part that varies by species and stocking, from around 0.05 m per metre in dense conifer plantations to 0.14 in open-grown broadleaf woodland.

import numpy as np
import pandas as pd


def fit_window_coefficients(plots: pd.DataFrame) -> dict:
    """Fit r = a + b*h from field crown widths.

    plots needs columns: height_m, crown_width_m (mean of two perpendicular widths).
    """
    h = plots["height_m"].to_numpy(dtype="float64")
    r = plots["crown_width_m"].to_numpy(dtype="float64") / 2.0
    b, a = np.polyfit(h, r, 1)
    pred = a + b * h
    ss_res = float(((r - pred) ** 2).sum())
    ss_tot = float(((r - r.mean()) ** 2).sum())
    return {
        "a": float(a),
        "b": float(b),
        "r2": 1.0 - ss_res / ss_tot if ss_tot else float("nan"),
        "rmse_m": float(np.sqrt(((r - pred) ** 2).mean())),
        "n": int(len(h)),
    }

An $R^2$ below about 0.4 is a warning that one line does not describe the stand — usually because two species with very different crown architectures are pooled. Fit them separately and pick the coefficients for whichever dominates, or accept that the detector will systematically over-segment one of them.

Minimal reproducible example

The filter itself is short. The one subtlety is that the search radius depends on the pixel being tested, so a single call to maximum_filter is not enough — you need one pass per distinct radius, then a per-pixel selection of the pass that matches that pixel’s own height.

import numpy as np
import rasterio
from scipy.ndimage import gaussian_filter, maximum_filter
from scipy import ndimage as ndi


def detect_tops(chm_path: str, a: float = 1.2, b: float = 0.09,
                min_height: float = 2.0, sigma: float = 0.6,
                radius_step_px: int = 2):
    """Variable-window local maxima detection. Returns (rows, cols, heights)."""
    with rasterio.open(chm_path) as src:
        chm = src.read(1, masked=True).astype("float32")
        res = abs(src.transform.a)
    h = chm.filled(0.0)
    h[h < 0] = 0.0
    h = gaussian_filter(h, sigma=sigma, mode="nearest")

    # Quantise radii so the filter bank stays small: cost scales with distinct radii.
    radii_px = np.round((a + b * h) / res / radius_step_px) * radius_step_px
    radii_px = np.clip(radii_px, 1, 25).astype("int16")

    tops = np.zeros(h.shape, dtype=bool)
    for k in np.unique(radii_px):
        local_max = maximum_filter(h, size=int(2 * k + 1), mode="nearest")
        tops |= (radii_px == k) & (h == local_max) & (h >= min_height)

    labels, n = ndi.label(tops, structure=np.ones((3, 3), dtype=int))
    if n == 0:
        return np.array([]), np.array([]), np.array([])
    centres = ndi.center_of_mass(tops, labels, range(1, n + 1))
    rows = np.array([int(round(r)) for r, _ in centres])
    cols = np.array([int(round(c)) for _, c in centres])
    return rows, cols, h[rows, cols]

Quantising the radius to every second pixel usually leaves five to eight passes over a temperate stand. Doubling the number of passes buys no measurable accuracy because the underlying allometry has an RMSE of half a metre or more anyway.

Parameter reference

Sweep the coefficients against a validation plot

Coefficients fitted from crown widths are a starting point, not an answer: field crown width and the apparent crown extent in a height raster are not the same quantity. Sweep b around the fitted value and pick the peak F₁ against a plot where every stem is mapped.

import numpy as np


def sweep_slope(chm_path, field, evaluate, b_values=np.arange(0.04, 0.16, 0.01),
                a: float = 1.2):
    """Grid-search the window slope; `evaluate` returns a dict with an 'f1' key."""
    results = []
    for b in b_values:
        rows, cols, heights = detect_tops(chm_path, a=a, b=float(b))
        scores = evaluate(rows, cols, heights, field)
        results.append({"b": float(b), **scores})
    best = max(results, key=lambda r: r["f1"])
    return best, results

The shape of that plot is the diagnosis. If precision is already high at the smallest slope, the CHM is smoother than you think and the real problem is recall. If recall never exceeds about 0.7 at any slope, the stand is multi-layered and the missing trees have no apex to detect — a raster method has reached its ceiling and the answer is point-cloud segmentation, not more tuning.

Expected output and verification

detect_tops returns three parallel arrays of equal length. Two cheap assertions catch nearly every wiring mistake before the result reaches a map:

rows, cols, heights = detect_tops("chm_0p5m.tif", a=1.2, b=0.09)

assert len(rows) == len(cols) == len(heights)
assert heights.min() >= 2.0, "a detection below the minimum-height threshold"

with rasterio.open("chm_0p5m.tif") as src:
    area_ha = (src.width * src.height * abs(src.transform.a) ** 2) / 10_000
density = len(rows) / area_ha
print(f"{len(rows)} trees over {area_ha:.1f} ha = {density:.0f} stems/ha")
assert 20 < density < 3000, "implausible stem density — check the window coefficients"

A density in the low hundreds is normal for mature temperate forest; four figures means the window is far too small for the stand, and a number below about 20 usually means the minimum-height threshold is filtering everything because the CHM is in the wrong units or still un-normalised.

Common pitfalls

• Fitting the allometry on dominant trees only. Plot crews often measure crown width on the plot’s largest stems, which biases the slope upward and produces a window that erases the understory. Sample the whole diameter range.
• Quantising the radius too coarsely. A radius_step_px of 4 or more collapses the window law into two or three classes, which is a fixed window in disguise.
• Running the filter before pit-filling. A CHM pit creates a ring of pixels that are all local maxima relative to the hole; every one becomes a tree.
• Comparing sweep results across different matching tolerances. F₁ is only comparable when max_dist and the height tolerance are held constant across the whole sweep.

Frequently Asked Questions

Related

• Individual Tree Detection and Crown Segmentation — the full detection-to-inventory workflow this step belongs to
• Watershed Crown Segmentation from a CHM in Python — turns these points into crown polygons
• Canopy Height Model Creation — the surface this filter runs over
• Removing CHM Pits with the Spike-Free Algorithm — fixes the artifact that produces spurious maxima

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