Python LiDAR

Rasterizing a Canopy Height Model from HAG

TL;DR: writers.gdal with "dimension": "HeightAboveGround", "output_type": "max" and a cell size your ground density can actually support — then clamp negatives to zero at the raster stage, not in the point cloud, so the evidence survives.

# Context and Motivation

This guide is part of Canopy Height Models with filters.hag_nn. The normalised cloud already knows how tall everything is; this page is about turning that into a raster without losing the information that made the point-domain route worth taking.

Three choices do all the work. The dimension, which decides whether you get a canopy height model or an accidental DSM. The reducer, which decides whether a cell reports its tallest return or a smoothed crown surface. And the cell size, which is not a quality dial but a statement about how much ground data supports the product — a one-metre CHM over a block with 0.4 ground returns per square metre is mostly interpolation wearing a resolution.

Four reducers over one canopy cell One cell containing eight returns spread through a crown. The max reducer reports the tallest at 26.4 metres, which is the canopy top. The mean reports 18.1, a value no physical surface sits at. The 95th percentile reports 25.1, more robust to a single high return. The count reports eight, which is a density layer rather than a height. one 1 m cell max → 26.4 m canopy top; the usual choice p95 → 25.1 m robust to one high return mean → 18.1 m no surface sits at this height count → 8 a density layer, not a height ask for several at once — "output_type": "max,count" writes both bands in a single pass over the cloud

# Prerequisites and Assumptions

Requirement Detail
A normalised cloud carrying HeightAboveGround, from filters.hag_nn
PDAL 2.4+ with writers.gdal
A defensible cell size supported by ground density, not chosen for appearance
GDAL for the post-processing clamp and for validation

# Step-by-Step Implementation

# Step 1 — Point the writer at the right dimension

json
{"type": "writers.gdal", "filename": "chm.tif",
 "dimension": "HeightAboveGround", "output_type": "max",
 "resolution": 1.0, "nodata": -9999, "gdaldriver": "GTiff"}

Omit dimension and you rasterize Z, producing a perfectly plausible DSM that nobody notices is wrong until it is compared against field measurements.

# Step 2 — Fill small gaps, and only small ones

json
{"window_size": 3}

Three cells fills the pinholes between crowns. Larger windows bridge real gaps and manufacture canopy where there is none — the same trade as in filling NoData voids in DTM rasters.

# Step 3 — Clamp negatives at the raster, not in the cloud

Small negative heights are interpolation noise. Clamping them to zero in the raster is honest presentation; deleting the points that produced them destroys the evidence that the ground surface was slightly wrong.

# Step 4 — Write a count band alongside

json
{"output_type": "max,count"}

The count band is the confidence layer: a cell with two returns and a cell with forty both report a height, and only one of them means much.

# Complete Working Example

python
"""Rasterize a canopy height model from a normalised cloud, with a count band."""
from __future__ import annotations

import json
import logging
from pathlib import Path

import numpy as np
import pdal

LOG = logging.getLogger("chm_raster")


def rasterize(src: Path, out: Path, resolution: float = 1.0, window: int = 3) -> int:
    spec = json.dumps({"pipeline": [
        {"type": "readers.las", "filename": str(src)},
        # Non-vegetation returns are heights above ground too; exclude them
        # explicitly rather than hoping the maximum happens to be a tree.
        {"type": "filters.range", "limits": "Classification[3:5]"},
        {"type": "writers.gdal", "filename": str(out),
         "dimension": "HeightAboveGround",
         "output_type": "max,count",
         "resolution": resolution,
         "window_size": window,
         "nodata": -9999,
         "gdaldriver": "GTiff",
         "gdalopts": "COMPRESS=DEFLATE,TILED=YES"},
    ]})
    return pdal.Pipeline(spec).execute()


def summarise(src: Path) -> dict:
    """Report the height distribution from the cloud, which the raster generalises."""
    p = pdal.Pipeline(json.dumps({"pipeline": [
        {"type": "readers.las", "filename": str(src)},
        {"type": "filters.range", "limits": "Classification[3:5]"},
    ]}))
    p.execute()
    hag = p.arrays[0]["HeightAboveGround"]
    return {
        "vegetation_points": int(len(hag)),
        "p50": round(float(np.percentile(hag, 50)), 2),
        "p95": round(float(np.percentile(hag, 95)), 2),
        "max": round(float(hag.max()), 2),
        "below_zero": int((hag < 0).sum()),
    }


if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO)
    n = rasterize(Path("normalised.laz"), Path("chm.tif"))
    LOG.info("rasterized %d vegetation returns", n)
    print(json.dumps(summarise(Path("normalised.laz")), indent=2))

# Key Parameter Table

Option Value Effect
dimension HeightAboveGround Without it the raster is a DSM
output_type max Canopy top; p95 where a single high return is a risk
resolution 1.0–2.0 m Set by ground density, not by how it looks
window_size 3 Fills pinholes; larger bridges real gaps
nodata −9999 Never 0, which is a legitimate height
gdalopts COMPRESS=DEFLATE,TILED=YES Cheap, and makes the product usable over a network
Ship the confidence band with the height band A two-band product. Band one is the maximum height above ground per cell, which is the canopy height model itself. Band two is the count of returns that produced it, which distinguishes a cell backed by forty returns from one backed by two. Without the second band nothing downstream can tell them apart. band 1 — max height above ground, metres float32, nodata −9999 the product everyone asks for band 2 — count returns contributing to the cell uint16, nodata 0 the one that says how much to trust it "output_type": "max,count" writes both in one pass — the cost is one extra band and the benefit is that a cell backed by two returns can be told apart from one backed by forty.

# Verification

The raster is heights, not elevations. Sample a cell over open ground: it should read near zero, not near the terrain elevation. This single check catches the missing dimension option.

The maximum is plausible. Compare against the tallest local species.

NoData is where you expect it. Over water and hard surfaces if you filtered to vegetation classes; not scattered randomly through the canopy, which would mean the cell size is too fine for the data.

The count band supports the height band. Cells with one or two returns should be treated as provisional, and a product that does not ship the count band cannot express that.

Cell size changes the answer, not just the picture Four cell sizes over the same block. At half a metre, 38 percent of cells are empty and the mean canopy height reads 24.9 metres. At one metre, 11 percent empty and 23.8. At two metres, 1 percent empty and 22.4. At four metres, no empty cells and 20.9. The coarser the grid, the lower the reported canopy, because a larger cell is more likely to include a gap. the same normalised cloud, four rasterizations 0.5 m 38% empty mean height 24.9 m 1.0 m 11% empty mean height 23.8 m 2.0 m 1% empty mean height 22.4 m 4.0 m 0% empty mean height 20.9 m a change-detection study that compares two epochs at different cell sizes measures the cell size

# Gotchas and Edge Cases

Non-vegetation returns are heights too. A building roof is 12 m above ground and will win the max reducer. Filter to Classification 3 to 5 if the product is about vegetation.

Zero is a legitimate height. Setting nodata to 0 makes every bare-earth cell disappear. Use −9999.

A large window_size invents canopy. The fill does not know that the gap it is bridging is a clearing.

Epochs must match. Comparing a 0.5 m CHM against a 2 m one measures the cell size, as the chart above shows.

# Frequently Asked Questions

Why does my canopy height model look like a surface model?

Because the dimension option was not set, so writers.gdal rasterized Z. The result is elevations above sea level and looks entirely reasonable, which is why it usually survives until someone compares it with field measurements.

Should I use max or mean as the reducer?

Max, in almost every case — it reports the tallest return in the cell, which is the canopy top. A mean over returns spread through a crown reports a height at which no physical surface sits. Where one spurious high return is a real risk, a high percentile is the robust compromise.

Does cell size change the numbers or just the picture?

The numbers. A larger cell is more likely to contain a gap, so mean canopy height falls as the grid coarsens — around four metres between a half-metre and a four-metre raster over the same block. Change detection between epochs is only meaningful if the cell sizes match.

What should NoData be?

Anything but zero. Zero is a legitimate canopy height over bare ground, so using it as the NoData value silently deletes every open cell from the product.

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... Extracting Individual Tree Heights from a CHM --- title: "Extracting Individual Tree Heights from a CHM" description: "Local maxima detection with a crown-scaled window, why every...