⚙️ How-To LiDAR & point clouds Terrain & elevation LiDAR & Point Clouds · step 12 of 14 Updated

How to Create a Canopy Height Model from LiDAR

Problem statement

A canopy height model is the difference between two surfaces — the top of everything and the bare ground — and it inherits every weakness of both.

Derived from a real 3DEP survey at 1 m cells:

DSM  16.0 .. 203.1 m,  7.47% of cells empty
DTM  15.9 .. 194.3 m, 10.93% of cells empty
CHM   0.0 ..  22.7 m, 10.93% undefined

The CHM's holes are exactly the DTM's, because a cell with no ground height has no height above ground. And that survey classifies only ground — there are no vegetation classes at all — so the height has to come from the returns rather than from a class.

Quick answer 

import numpy as np


def canopy_height_model(x, y, z, classification, cell=1.0, ground_class=2):
    """DSM minus DTM, per cell, with the undefined cells left undefined."""
    left, top = x.min(), y.max()
    width = int(np.ceil((x.max() - left) / cell))
    height = int(np.ceil((top - y.min()) / cell))

    col = np.clip(((x - left) / cell).astype(int), 0, width - 1)
    row = np.clip(((top - y) / cell).astype(int), 0, height - 1)
    flat = row * width + col
    ground = classification == ground_class

    dsm = np.full(width * height, -np.inf)
    np.maximum.at(dsm, flat, z)
    dsm[np.isinf(dsm)] = np.nan

    dtm = np.full(width * height, np.inf)
    np.minimum.at(dtm, flat[ground], z[ground])
    dtm[np.isinf(dtm)] = np.nan

    chm = (dsm - dtm).reshape(height, width)
    print(f"  CHM {np.nanmin(chm):.1f}..{np.nanmax(chm):.1f} m, "
          f"{np.isnan(chm).mean():.2%} undefined")
    return chm

Note what this cannot produce: a negative value. The DSM is a maximum over all returns and the DTM a minimum over a subset of them, so DSM >= DTM in every cell by construction — and the measurement confirms zero negative cells.

A canopy height model as the difference between a DSM of maximum returns and a DTM of minimum ground returns, undefined where no ground point exists.
Two reductions and a subtraction. The undefined cells come entirely from the DTM.

Step-by-step solution

1. Reproject first

Everything here is metric — cell size, heights, the eventual crown diameters. The survey used here is stored in Web Mercator, where distances are inflated by 1.92 at its latitude. Heights are unaffected, so a CHM built without reprojecting has correct heights on cells the wrong size.

2. Build both surfaces on the same grid

Not two grids that happen to have the same cell size. Compute the origin and shape once and use them for both, so the subtraction is aligned by construction rather than by luck.

3. Choose a cell size the ground can support 

0.25 m:  78.32% of cells have no ground point
0.50 m:  38.11%
1.00 m:  10.93%
2.00 m:   7.51%

The CHM is limited by the DTM, which is limited by the ground density — 3.9 points per square metre here against 16.9 overall. A 0.5 m CHM would be undefined over more than a third of the area.

4. Decide about the DTM holes before subtracting

Three options, with different consequences:

  • Leave them. The CHM is undefined there. Honest and awkward.
  • Interpolate the DTM, then subtract. Fills the CHM, and introduces the possibility of negative values.
  • Interpolate the CHM afterwards. Worse: it interpolates a difference rather than a surface, and difference fields are rougher.

Interpolate the DTM, and ship a mask.

5. Expect negative values once you interpolate

The per-cell CHM cannot be negative. An interpolated one can, wherever the filled ground height exceeds the highest return in that cell — most often at a break of slope, where a linear fill overshoots.

Negative CHM is a diagnostic about the DTM, not about the canopy. Clipping it to zero throws the diagnostic away; investigate first, then clip if the magnitudes are small.

6. Smooth only if you know why

A CHM built from per-cell maxima is noisy, because the maximum is a single point. Smoothing makes tree crowns look like tree crowns, and it also lowers every peak — which matters if you are measuring tree height.

A per-cell CHM bounded below by zero, and an interpolated version going negative at a break of slope where the filled ground exceeds the highest return.
Negatives appear only after filling. They are a message about the ground surface.

Code examples

Example 1 — CHM with filling, masking and a negative check 

import numpy as np
import rasterio
from rasterio.transform import from_origin
from scipy.ndimage import distance_transform_edt


def chm_with_mask(x, y, z, classification, cell=1.0, crs="EPSG:32605",
                  max_fill_m=5.0, out_path=None):
    """CHM, an observation mask, and a report on any negative values."""
    left, top = x.min(), y.max()
    width = int(np.ceil((x.max() - left) / cell))
    height = int(np.ceil((top - y.min()) / cell))
    transform = from_origin(left, top, cell, cell)

    col = np.clip(((x - left) / cell).astype(int), 0, width - 1)
    row = np.clip(((top - y) / cell).astype(int), 0, height - 1)
    flat = row * width + col
    ground = classification == 2

    dsm = np.full(width * height, -np.inf); np.maximum.at(dsm, flat, z)
    dsm[np.isinf(dsm)] = np.nan
    dtm = np.full(width * height, np.inf)
    np.minimum.at(dtm, flat[ground], z[ground]); dtm[np.isinf(dtm)] = np.nan

    dsm = dsm.reshape(height, width)
    dtm = dtm.reshape(height, width)
    observed = np.isfinite(dtm)

    indices = distance_transform_edt(~observed, return_distances=False,
                                     return_indices=True)
    distance = distance_transform_edt(~observed) * cell
    dtm_filled = dtm[tuple(indices)]
    dtm_filled[distance > max_fill_m] = np.nan

    chm = dsm - dtm_filled
    negative = np.nansum(chm < 0)
    if negative:
        print(f"  {int(negative):,} negative cells "
              f"({negative / np.isfinite(chm).sum():.2%}), "
              f"minimum {np.nanmin(chm):.2f} m — inspect the DTM fill")
    print(f"  CHM {np.nanmin(chm):.1f}..{np.nanmax(chm):.1f} m, "
          f"{np.isnan(chm).mean():.2%} undefined, "
          f"{(~observed).mean():.2%} of ground interpolated")

    if out_path:
        profile = dict(driver="GTiff", height=height, width=width, count=2,
                       dtype="float32", crs=crs, transform=transform,
                       nodata=np.nan, compress="deflate", tiled=True)
        with rasterio.open(out_path, "w", **profile) as dst:
            dst.write(chm.astype("float32"), 1)
            dst.write(observed.astype("float32"), 2)
            dst.set_band_description(1, "canopy height (m)")
            dst.set_band_description(2, "1 = ground observed, 0 = interpolated")
    return chm, observed, transform

Example 2 — canopy metrics from normalised points, not from the CHM 

import numpy as np


def canopy_metrics(x, y, z, dtm, transform, cell_dtm, cell_out=10.0,
                   cover_threshold=2.0, min_points=10):
    """Cover and height percentiles, which a CHM has already thrown away."""
    left, top = transform.c, transform.f
    h, w = dtm.shape
    col = np.clip(((x - left) / cell_dtm).astype(int), 0, w - 1)
    row = np.clip(((top - y) / cell_dtm).astype(int), 0, h - 1)
    above = z - dtm[row, col]

    ok = np.isfinite(above)
    x, y, above = x[ok], y[ok], above[ok]

    width = int(np.ceil((x.max() - x.min()) / cell_out))
    height = int(np.ceil((y.max() - y.min()) / cell_out))
    oc = np.clip(((x - x.min()) / cell_out).astype(int), 0, width - 1)
    orr = np.clip(((y.max() - y) / cell_out).astype(int), 0, height - 1)
    flat = orr * width + oc

    order = np.argsort(flat)
    fs, hs = flat[order], above[order]
    edges = np.searchsorted(fs, np.arange(width * height + 1))

    p95 = np.full(width * height, np.nan)
    cover = np.full(width * height, np.nan)
    for i in range(width * height):
        a, b = edges[i], edges[i + 1]
        if b - a < min_points:
            continue
        cell_heights = hs[a:b]
        p95[i] = np.percentile(cell_heights, 95)
        cover[i] = (cell_heights > cover_threshold).mean()

    print(f"  {cell_out} m cells: median p95 height "
          f"{np.nanmedian(p95):.1f} m, median cover {np.nanmedian(cover):.1%}")
    return p95.reshape(height, width), cover.reshape(height, width)

Canopy cover — the fraction of returns above a threshold — cannot be recovered from a CHM at all, because the CHM has reduced each cell to one number. Working from normalised points keeps the whole vertical distribution.

Example 3 — individual tree detection, with its assumptions stated 

import numpy as np
from scipy.ndimage import maximum_filter, label


def detect_trees(chm, cell=1.0, min_height=3.0, window_m=5.0):
    """Local maxima as tree tops. A crude method with visible assumptions."""
    size = max(3, int(round(window_m / cell)) | 1)          # odd window
    smoothed = np.where(np.isfinite(chm), chm, 0.0)

    peaks = (maximum_filter(smoothed, size=size) == smoothed) & \
            (chm > min_height) & np.isfinite(chm)
    labels, n = label(peaks)

    rows, cols = np.nonzero(peaks)
    heights = chm[rows, cols]
    print(f"  {n} local maxima above {min_height} m "
          f"with a {size}x{size} window ({window_m} m)")
    print(f"  heights {heights.min():.1f}..{heights.max():.1f} m, "
          f"median {np.median(heights):.1f} m")
    print(f"  density {n / (chm.size * cell ** 2 / 1e4):.1f} per hectare")
    return rows, cols, heights

The window size is the assumption: it asserts a minimum crown spacing. Too small and one tree becomes several; too large and a group becomes one. Report the window with the tree count, and expect the count to move substantially with it.

Explanation

Why the CHM inherits the DTM's holes and not the DSM's

Height above ground needs a ground height. A cell with returns but no ground point has a DSM value and no DTM value, so the difference is undefined.

Measured, that is the gap between 7.47% and 10.93% — 3.46% of cells where the laser reached something but never the ground. Under closed canopy that fraction rises steeply, which is why forest CHMs rest on more interpolation than open-ground ones.

Why a per-cell CHM cannot be negative

DSM = max(z over all returns in the cell), DTM = min(z over ground returns in the cell). The ground returns are a subset of all returns, so the maximum over the superset is at least the minimum over the subset.

The measurement confirms it: zero negative cells out of the 89% that are defined.

Once the DTM is filled or smoothed, the ground height in a cell no longer comes from that cell's own points, and the guarantee is gone. Negatives then mark places where the fill overshot — typically the top of a bank or the edge of a terrace.

Why maxima make a noisy surface

A DSM cell is a single point: the highest return that happened to land there. Neighbouring cells take different points from different pulses, so the surface is speckled at the scale of the point spacing.

Smoothing fixes the appearance and biases the heights downward, because a maximum is at the top of the local distribution and any average of it with neighbours is lower. If tree height is the output, smooth as little as possible; if the CHM is a background layer, smooth freely.

Why normalised points beat a CHM for structure

A CHM answers "how tall is the tallest thing here". Most forestry questions are about the distribution — cover, layering, the height below which 50% of returns fall — and a per-cell maximum has discarded all of it.

Normalising the cloud (z - dtm[row, col]) keeps every return with its height above ground, and every metric becomes a per-cell statistic over that array. It is more work and strictly more informative.

A canopy height model reducing each cell to one number against a normalised point cloud retaining the vertical distribution.
For height, a CHM is enough. For cover and layering, normalise the cloud instead.

Edge cases or notes

  • The CHM's holes are the DTM's holes. Fix the ground surface, not the difference.
  • Build both surfaces on one grid computed once.
  • A per-cell CHM cannot be negative. Negatives mean the DTM was filled.
  • Do not clip negatives to zero without looking at where they are.
  • Interpolate the DTM, not the CHM. A difference field is rougher than a surface.
  • Smoothing lowers every peak. Avoid it if height is the measurement.
  • Cell size is limited by ground density, not overall density.
  • Water has no returns, so both surfaces are undefined there.

FAQ

How do I make a canopy height model in Python?

Build a DSM as the per-cell maximum of all returns and a DTM as the per-cell minimum of ground returns on the same grid, then subtract.

Why is my CHM full of holes?

Because the DTM is. A cell with returns but no ground point has no height above ground — 3.46% of cells here, on top of 7.47% with no returns at all.

Why does my CHM have negative values?

Because the DTM was interpolated or smoothed, so the ground height in a cell is no longer bounded by that cell's own points. A per-cell CHM cannot go negative.

Should I clip negative CHM values to zero?

Only after investigating. They point at where the ground fill overshot, usually at a break of slope, and clipping hides that.

What cell size should a CHM use?

One the ground density supports. On this survey, 1 m left 10.93% undefined and 0.5 m would have left 38%.

Should I smooth the CHM?

Only if it is a background layer. Smoothing lowers every peak, which biases tree heights downward.

How do I get canopy cover rather than height?

Normalise the point cloud against the DTM and compute the fraction of returns above a threshold per cell. A CHM has already reduced each cell to a single number.