Computing Height Above Ground with filters.hag_nn

TL;DR: {"type": "filters.hag_nn", "count": 6, "max_distance": 25.0, "allow_extrapolation": false} after a ground classification, then check that ground points come back at a height of zero — if they do not, the problem is the classifier, not the filter.

# Context and Motivation

This guide is part of Canopy Height Models with filters.hag_nn. The parent covers the whole workflow; this page is about the one stage that does the arithmetic, because nearly every bad canopy height traces back to how it was configured.

What the stage does is simple enough to state exactly. For each point not classified as ground, it finds the count nearest points that are classified as ground, averages their elevations to get an estimated ground level beneath the point, and writes the difference into a new HeightAboveGround dimension. Every part of that sentence is a place where a decision has to be made, and PDAL’s defaults are not the right ones for forestry.

# Prerequisites and Assumptions

# Step-by-Step Implementation

# Step 1 — Confirm ground exists

pdal info tile.laz --stats --dimensions Classification | grep -i counts

If class 2 is absent, hag_nn writes zero everywhere and reports success.

# Step 2 — Choose count

Six is a good default. One makes each height depend on a single ground return. Above about eight the surface stops improving and the search cost keeps rising.

# Step 3 — Bound the search with max_distance

{"type": "filters.hag_nn", "count": 6, "max_distance": 25.0}

Without it, a point over a large ground-free area searches arbitrarily far and produces a height measured against terrain a hundred metres away.

# Step 4 — Decide about extrapolation deliberately

allow_extrapolation: false leaves points beyond the ground hull unestimated. That is usually right: a gap you can see beats a number you cannot trust.

# Step 5 — Persist the dimension

{"type": "writers.las", "extra_dims": "HeightAboveGround=float",
 "minor_version": 4, "dataformat_id": 6}

# Complete Working Example

"""Compute height above ground and verify it against the ground class."""
from __future__ import annotations

import json
import logging
from pathlib import Path

import numpy as np
import pdal

LOG = logging.getLogger("hag")


def normalise(src: Path, dst: Path, count: int = 6, max_distance: float = 25.0) -> int:
    spec = json.dumps({"pipeline": [
        {"type": "readers.las", "filename": str(src)},
        {"type": "filters.hag_nn", "count": count,
         "max_distance": max_distance, "allow_extrapolation": False},
        {"type": "writers.las", "filename": str(dst), "compression": "laszip",
         "minor_version": 4, "dataformat_id": 6,
         "extra_dims": "HeightAboveGround=float", "forward": "all"},
    ]})
    return pdal.Pipeline(spec).execute()


def verify(path: Path) -> dict:
    p = pdal.Pipeline(json.dumps({"pipeline": [
        {"type": "readers.las", "filename": str(path)}]}))
    p.execute()
    arr = p.arrays[0]

    if "HeightAboveGround" not in arr.dtype.names:
        raise AssertionError("dimension absent — extra_dims was not set on the writer")

    hag = arr["HeightAboveGround"]
    ground = hag[arr["Classification"] == 2]
    if len(ground) == 0:
        raise AssertionError("no ground-classified points — classify before running hag_nn")

    ground_median = float(np.median(ground))
    if abs(ground_median) > 0.10:
        raise AssertionError(
            f"ground sits at {ground_median:.2f} m rather than zero — the classifier is suspect"
        )

    return {
        "points": int(len(hag)),
        "ground_median": round(ground_median, 3),
        "negative_fraction": round(float((hag < -0.5).mean()), 5),
        "p99": round(float(np.percentile(hag, 99)), 2),
        "max": round(float(hag.max()), 2),
    }


if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO)
    normalise(Path("classified.laz"), Path("normalised.laz"))
    print(json.dumps(verify(Path("normalised.laz")), indent=2))

# Key Parameter Table

# Verification

Ground is at zero. The strongest single check, asserted above. A systematic offset means the ground class is wrong, not the filter.

The dimension reached the file. Also asserted — the failure mode where everything works in memory and nothing is written.

Negative heights are rare. A handful is normal around breaklines. A percent or more below −0.5 m means blunders in the ground class.

Heights are physically plausible. The 99th percentile against the tallest species locally.

# Gotchas and Edge Cases

All heights zero, no error. No ground class. This is the single most common report, and the fix is upstream.

A ring of tall values around a clearing. Extrapolation at the hull edge. Turn it off.

The stage dominates the runtime. Expected: it is a nearest-neighbour search per non-ground point. Reduce the input before it rather than tuning it.

Buildings become “canopy”. hag_nn measures height above ground regardless of what the point struck. If the product is about vegetation, filter to the vegetation classes after normalising — the codes are in understanding ASPRS classification codes.

# Frequently Asked Questions

Why is every height above ground zero?

Because no points carry Classification 2. The filter measures against ground-classified points, and with none present it has nothing to subtract, so it writes zero and the pipeline succeeds. Classification has to happen before the stage, in the same pipeline or an earlier one.

What does the count parameter actually change?

How many ground returns are averaged to estimate the ground beneath each point. With one, every height inherits that single return’s error at full amplitude. With six the estimate is smooth and still follows real breaks in slope. Above about eight the surface stops improving and starts cutting corners on concave terrain.

Should I set max_distance?

Always. Without it a point over a large ground-free area searches arbitrarily far and ends up measured against terrain a hundred metres away, which produces a confident and meaningless number. Fifteen to thirty metres suits most airborne forestry work.

Why do buildings appear in my canopy height model?

Because height above ground is exactly that — it does not care what the point struck. Filter to the vegetation classes after normalising if the product is about vegetation; the stage itself is deliberately agnostic.

# Related

• Canopy Height Models with filters.hag_nn — the parent workflow this stage sits inside
• Rasterizing a Canopy Height Model from HAG — turning the normalised cloud into a raster
• Extracting Individual Tree Heights from a CHM — what the raster can and cannot tell you about single trees
• SMRF Ground Classification — the classification every height is measured against
• Measuring Ground Point Density Under Canopy — whether enough ground returns exist to measure against at all