Canopy Height Models with filters.hag_nn

A digital surface model tells you how high something is above sea level. A canopy height model tells you how tall it is — and almost every question about vegetation, from biomass to fuel loading to whether a tree will reach a power line, is a question about the second. Getting there means subtracting the ground from the surface, and PDAL offers two routes: raster arithmetic on a DSM and a DTM, or filters.hag_nn, which computes a height above ground for every point before any rasterization happens. This topic belongs to Ground Filtering and DTM/DSM Generation with PDAL.

The point-domain route is the better one, and the reason is resolution. Raster subtraction inherits the cell size of both inputs, so a canopy gap narrower than a cell disappears before the arithmetic starts. Computing height above ground per point keeps every return’s own height, and you rasterize once, at the end, from a cloud that already knows how tall everything is. That also makes the result reusable: the same normalised cloud answers “how tall” and “how many returns above two metres” without recomputation.

# Prerequisites

# Core Workflow Architecture

1. Classify ground. filters.smrf or filters.pmf, tuned for the vegetation — the whole model rests on this, and an error here propagates to every canopy height.
2. Build the ground index. filters.hag_nn collects the ground-classified points into a nearest-neighbour structure.
3. Interpolate a ground elevation per point. For each non-ground point, the filter averages the elevations of its count nearest ground neighbours.
4. Subtract. HeightAboveGround is written as a new dimension: the point’s Z minus that interpolated ground.
5. Decide about extrapolation. Where a point has no ground neighbours within reach, the filter either extrapolates or writes zero, depending on allow_extrapolation.
6. Rasterize. writers.gdal with output_type: "max" over HeightAboveGround produces the canopy height model.

# Full Implementation

"""Produce a canopy height model through the point domain."""
from __future__ import annotations

import json
import logging
from pathlib import Path

import numpy as np
import pdal

LOG = logging.getLogger("chm")


def build_chm(src: Path, chm: Path, normalised: Path, resolution: float = 1.0) -> dict:
    spec = json.dumps({"pipeline": [
        {"type": "readers.las", "filename": str(src)},
        # Blunders first: a single low noise point becomes a negative height
        # and, worse, can be picked up as a ground neighbour.
        {"type": "filters.outlier", "method": "statistical",
         "mean_k": 12, "multiplier": 2.5},
        {"type": "filters.range", "limits": "Classification![7:7]"},
        {"type": "filters.smrf", "window": 33, "slope": 0.2,
         "threshold": 0.6, "cell": 1.0, "returns": "last, only"},
        {"type": "filters.hag_nn", "count": 6, "allow_extrapolation": False},
        # Keep the normalised cloud: it answers more questions than the raster does.
        {"type": "writers.las", "filename": str(normalised),
         "compression": "laszip", "minor_version": 4, "dataformat_id": 6,
         "extra_dims": "HeightAboveGround=float", "forward": "all"},
        {"type": "writers.gdal", "filename": str(chm),
         "dimension": "HeightAboveGround", "output_type": "max",
         "resolution": resolution, "window_size": 3, "nodata": -9999,
         "gdaldriver": "GTiff"},
    ]})
    written = pdal.Pipeline(spec).execute()
    LOG.info("wrote %d points and a %.1f m CHM", written, resolution)
    return {"points": written, "chm": str(chm), "normalised": str(normalised)}


def audit(normalised: Path) -> dict:
    p = pdal.Pipeline(json.dumps({"pipeline": [
        {"type": "readers.las", "filename": str(normalised)}]}))
    p.execute()
    hag = p.arrays[0]["HeightAboveGround"]
    return {
        "min": round(float(hag.min()), 2),
        "p99": round(float(np.percentile(hag, 99)), 2),
        "max": round(float(hag.max()), 2),
        "negative_fraction": round(float((hag < -0.5).mean()), 5),
    }


if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO)
    build_chm(Path("forest.laz"), Path("chm.tif"), Path("forest_hag.laz"))
    print(json.dumps(audit(Path("forest_hag.laz")), indent=2))

# Code Breakdown

Outlier removal comes before classification, not after. A low blunder that survives into the ground class becomes the nearest ground neighbour for everything around it, and every canopy height in that neighbourhood is inflated by the blunder’s depth.

returns: "last, only" on the classifier. Under canopy the ground is reconstructed from the returns that penetrated, which is the point made at length in tuning SMRF for forested terrain.

count: 6 rather than 1. Averaging six ground neighbours smooths the interpolated surface. A count of one makes every canopy height depend on a single ground return, and that return’s own error appears in the result at full amplitude.

allow_extrapolation: False is a deliberate, conservative choice. Where no ground neighbours are within reach — the middle of a large closed-canopy patch — extrapolating invents a ground surface from distant points. Writing zero instead makes the gap visible rather than plausible.

The normalised cloud is written as well as the raster. It is the reusable artefact; the CHM is one view of it.

# Parameter Reference Table

# Validation and Integrity Checks

Ground points have height zero. By construction, every Classification 2 point should have a HeightAboveGround within a few centimetres of zero. A systematic offset means the interpolation is wrong.

ground = arr[arr["Classification"] == 2]["HeightAboveGround"]
assert abs(float(np.median(ground))) < 0.10, "ground is not at zero — check the classifier"

Negative heights are rare. A small number is normal — interpolation noise around breaklines. More than a fraction of a percent below −0.5 m means blunders survived into the ground class.

The maximum is physically plausible. Compare the 99th percentile against the tallest species in the area. A CHM with 60 m trees in a region whose canopy tops out at 35 m is reporting noise.

Open ground is near zero. Sample cells known to be a field or a car park; the CHM should read close to zero there, and a non-zero floor across open ground indicates a ground surface sitting below the true one.

# Performance Tuning

filters.hag_nn is the expensive stage, and its cost is a spatial search per non-ground point against a structure built from the ground points. Three levers matter.

Reduce the non-ground point count first. Elevation limits and a crop cost almost nothing and shrink the number of searches directly — the ordering argument made in reordering PDAL stages for speed.

Keep count modest. Six neighbours is roughly six times the search work of one. Beyond about eight the surface stops getting smoother and the runtime keeps rising.

Set max_distance. Without it, a point over a large ground-free area searches unboundedly far. With it, the search terminates and the point is left unhandled, which is both faster and more honest.

The stage does not stream, so peak memory holds the whole cloud plus the ground index. On large tiles the split-pass pattern applies: stream the reduction, run HAG on the reduced cloud, stream the write.

# What a Canopy Height Model Is and Is Not

A CHM is a raster of height above the terrain, and the distance between that sentence and the questions people ask of it accounts for most of the disappointment in vegetation work.

It is not a tree height map. A cell holds the tallest return that fell in it, which for a cell straddling two crowns is the taller of the two, and for a cell between crowns is a branch tip or a patch of understory. Individual tree height needs the crowns delineated first, which is a separate operation on top of the CHM rather than a reading of it.

It is not a measure of the true treetop either. A laser pulse strikes somewhere on the upper crown, rarely the apex. Airborne LiDAR under-measures individual tree height by a metre or two as a rule, and the bias grows with narrower crowns and lower pulse density. For relative comparisons across a block this hardly matters; for a stand table checked against field measurements it matters a great deal, and the usual remedy is a locally calibrated offset rather than a change of processing.

It is not comparable across acquisitions unless both were normalised the same way. Two CHMs of the same forest, one built at 0.5 m from a 20 pts/m² flight and one at 2 m from an 8 pts/m² flight, will disagree systematically — the coarser one lower, because a larger cell is more likely to include a gap. Change detection between epochs requires matching cell sizes and comparable densities, and where those differ the difference has to be modelled rather than ignored.

It is a superb relative measure. Where the canopy is tall, where it is short, where a gap has opened since the last flight, how fuel load varies across a slope — the CHM answers all of these well, because the systematic biases affect the whole raster nearly equally and cancel in comparison.

# Products Derived from the Normalised Cloud

Writing the normalised cloud rather than only the raster is what makes the rest of this cheap. Once every return carries a HeightAboveGround, several standard products are one more pass each, with no reclassification and no second ground filter.

Canopy cover is the fraction of first returns above a height threshold, typically two metres, per cell. It is a filters.range on HeightAboveGround and ReturnNumber followed by a count raster, divided by an unfiltered count raster.

Height percentiles — the 25th, 75th and 95th percentile height per cell — are the standard predictors in area-based biomass models. They come from the same normalised cloud with a percentile reducer, and because they are computed from the point distribution rather than from a surface they are far more robust to gaps than the CHM itself.

Understory density is the count of returns between, say, 0.5 and 3 metres, which is directly interesting for fuel modelling and for habitat work, and impossible to recover from a CHM because the CHM kept only the maximum.

Vertical profiles for a stand come from histogramming HeightAboveGround over an area of interest, which is one NumPy call against the normalised cloud.

All four share the same ground classification, so they are internally consistent by construction — a property that is easy to lose when each product is built by a separate pipeline that classifies ground its own way.

# Common Errors and Troubleshooting

Every height is zero. No points carry Classification 2, so the filter had no ground to measure against. Classification must precede HAG in the same pipeline or in an earlier one.

Heights are systematically too tall. A ground surface sitting below the true one, almost always from a low blunder that survived into the ground class. Run outlier removal before classifying.

A ring of extreme heights around a clearing. Extrapolation working as designed at the edge of the ground points’ hull. Set allow_extrapolation: False and accept the gaps.

The raster is a DSM, not a CHM. dimension was not set on the writer, so it rasterized Z. The output looks entirely plausible and is elevations above sea level.

Runtime is dominated by one stage. Expected — HAG is the search. Reduce the input before it rather than tuning it.

# Frequently Asked Questions

Why compute height above ground per point instead of subtracting rasters?

Resolution. Raster subtraction inherits the cell size of both inputs, so anything narrower than a cell is generalised away before the arithmetic starts. Computing height per return keeps every point’s own height and defers generalisation to a single rasterization at the end.

Why is my canopy height model systematically too tall?

Almost always because the ground surface sits below the true ground. A single low blunder that survived into the ground class becomes the nearest ground neighbour for its whole neighbourhood, and every height there is inflated by the blunder’s depth. Run outlier removal before classification.

Should I allow extrapolation?

Usually not. Where a point has no ground neighbours within reach — the middle of a large closed-canopy patch — extrapolating invents a ground surface from distant points and produces a plausible-looking number with no support. Leaving the gap visible is more useful than filling it convincingly.

Why does my CHM look exactly like a DSM?

Because the writer rasterized Z. Setting dimension to HeightAboveGround on writers.gdal is what makes the raster a height model; without it the output is elevations above sea level and looks entirely reasonable.

# Related

• Ground Filtering and DTM/DSM Generation with PDAL — the section this workflow belongs to
• Computing Height Above Ground with filters.hag_nn — the stage itself, its parameters and its failure modes
• Rasterizing a Canopy Height Model from HAG — turning the normalised cloud into a raster product
• Extracting Individual Tree Heights from a CHM — local maxima, smoothing and what the numbers mean
• DTM vs DSM: Which Surface Model — the raster route and when it is enough
• SMRF Ground Classification — the classification every canopy height rests on