Building a Pit-Free Canopy Height Model

TL;DR: Normalize heights, then build several partial CHMs: one from all first returns and one each from first returns above 2, 5, 10, 15, 20 m … . Interpolate each partial CHM with a small search radius so gaps between high returns are not filled with low values, and take the per-cell maximum across the stack. Pits — low values punched into crowns by pulses that penetrated gaps — disappear because each crown is also represented by a layer that never saw the lower returns.

# Context and Motivation

This guide is part of Canopy Height Models. A CHM built as the maximum height above ground per cell has pits: cells inside crowns whose highest return came from a branch or the ground below, because the laser passed through a gap. Pits fragment crowns and break individual tree segmentation, which relies on smooth crowns with a single peak. Khosravipour, Skidmore, Isenburg and colleagues (2014) proposed a simple fix: compute partial CHMs using only returns above a series of height thresholds, so that within a tall crown the lower returns that caused pits are excluded, then combine the layers by taking the highest value per cell. The approach is easy to implement with PDAL and NumPy, and it improves tree detection noticeably in open-canopy forests.

# Prerequisites and Assumptions

• A ground-classified tile for height normalization.
• PDAL with writers.gdal, and NumPy and rasterio.
• Density of about 5 pts/m² or more; very sparse data produces sparse upper layers.
• A target CHM resolution, commonly 0.5 m.

# Step-by-Step Implementation

# Step 1 — Normalize heights

filters.hag_nn computes HeightAboveGround; ferry it into Z so rasterization works on heights.

# Step 2 — Choose thresholds

0 (all first returns), 2, 5, 10, 15, 20, 25 … up to the canopy’s maximum height, in steps of about 5 m.

# Step 3 — Rasterize each layer with a tight radius

For each threshold, keep first returns above it and rasterize the maximum with a radius of about one cell. The tight radius is what stops a layer from spreading high values into gaps; the original method uses a TIN with an edge-length limit for the same purpose.

# Step 4 — Stack by maximum

Read all layers and take np.fmax across them, ignoring NoData.

# Step 5 — Fill remaining gaps lightly

Cells still empty (true gaps between crowns) take the all-returns value, which is usually near zero there.

# Complete Working Example

"""Pit-free CHM by stacking partial CHMs from successive height thresholds."""
from __future__ import annotations

import json
import math
from pathlib import Path

import numpy as np
import pdal
import rasterio

SRC = "tiles/forest_0822.laz"
RES = 0.5
THRESHOLDS = [0, 2, 5, 10, 15, 20, 25, 30]
OUT = Path("out/chm_layers")


def grid(src: str) -> dict:
    b = pdal.Pipeline(json.dumps({"pipeline": [src]})).quickinfo["readers.las"]["bounds"]
    ox, oy = math.floor(b["minx"]), math.floor(b["miny"])
    return {"origin_x": ox, "origin_y": oy,
            "width": math.ceil((b["maxx"] - ox) / RES), "height": math.ceil((b["maxy"] - oy) / RES)}


def build_layers(src: str) -> list[Path]:
    OUT.mkdir(parents=True, exist_ok=True)
    g = grid(src)
    stages = [
        src,
        {"type": "filters.range", "limits": "Classification![7:7],Classification![18:18]"},
        {"type": "filters.hag_nn", "count": 2},
        {"type": "filters.ferry", "dimensions": "HeightAboveGround=>Z"},
        {"type": "filters.range", "limits": "ReturnNumber[1:1],Z[0:80]", "tag": "first"},
    ]
    paths = []
    for t in THRESHOLDS:
        path = OUT / f"chm_gt{t:02d}.tif"
        stages += [
            {"type": "filters.range", "inputs": ["first"], "limits": f"Z[{t}:80]", "tag": f"gt{t}"},
            {"type": "writers.gdal", "inputs": [f"gt{t}"], "filename": str(path), "resolution": RES,
             "radius": RES * 1.0, "output_type": "max", "data_type": "float32", "nodata": -9999, **g},
        ]
        paths.append(path)
    pdal.Pipeline(json.dumps({"pipeline": stages})).execute()
    return paths


def stack(paths: list[Path], dst: Path) -> np.ndarray:
    layers = []
    for p in paths:
        with rasterio.open(p) as ds:
            layers.append(ds.read(1, masked=True).filled(np.nan))
            profile = ds.profile
    chm = np.nanmax(np.stack(layers), axis=0)          # per-cell maximum across layers
    chm = np.where(np.isnan(chm), 0.0, np.clip(chm, 0, None))
    with rasterio.open(dst, "w", **profile) as out:
        out.write(chm.astype("float32"), 1)
    return chm


if __name__ == "__main__":
    layers = build_layers(SRC)
    chm = stack(layers, Path("out/chm_pitfree.tif"))
    print(f"pit-free CHM: max {chm.max():.1f} m, canopy cells {np.mean(chm > 2):.1%}")

All layers are written in one branched pipeline, so the tile is read and normalized once — see branching a PDAL pipeline with tags.

# Key Parameter Table

# Verification

• Pit count. Count cells more than 2 m below the 3×3 median in the standard and pit-free CHMs; the pit-free count should fall sharply.
• Heights preserved. Tree tops (local maxima) should have the same heights in both CHMs; the method removes pits, it does not raise tops.
• Tree segmentation. Run watershed segmentation on both and compare against plots.

# Gotchas and Edge Cases

Radius too large. A generous radius in upper layers spreads crown-top values sideways into gaps between trees, merging neighbouring crowns. Keep it near one cell.

Sparse upper layers. On low-density data, the highest layers contain only a handful of returns and add speckle. Stop the thresholds where layers become too sparse to form crowns.

Thresholds relative to tree size. Short stands (under 10 m) need finer steps (1–2 m); tall forest can use 5 m steps.

Buildings. In mixed urban areas, roofs appear in the CHM too. Exclude class 6 before building layers if the CHM is meant for vegetation only.

# Frequently Asked Questions

What is a pit-free canopy height model?

A CHM in which the spurious low cells inside crowns, caused by laser pulses passing through gaps, have been removed. The layered method builds partial CHMs from returns above successive heights and combines them with a per-cell maximum.

Why do pits matter?

They fragment crowns into several apparent peaks and holes, which makes individual tree detection count too many trees and distorts crown metrics and canopy cover.

Is a pit-free CHM the same as smoothing?

No. Smoothing lowers tree tops and blurs crown edges. The layered method keeps measured heights and only replaces pit cells with values from higher layers.

How many layers do I need?

Enough to span the canopy height range in steps of about 5 metres for tall forest, or 1 to 2 metres for short vegetation. Beyond the tallest trees, extra layers are empty and add nothing.

# Related

• Canopy Height Models — CHMs in general
• Rasterizing a Canopy Height Model from HAG — the standard CHM
• Segmenting Trees with a Watershed on a CHM — the main consumer of pit-free CHMs
• Removing Pits and Spikes from a DSM — the median-based alternative
• Computing Canopy Cover from LiDAR — another CHM-derived product