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, counting only the cells the survey covered (the surfaces guide measures them):

DSM  16.0 .. 203.1 m,  0.00% of cells empty
DTM  15.9 .. 194.3 m,  3.67% of cells empty
CHM   0.0 ..  22.7 m,  3.67% 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


chm = canopy_height_model(x, y, z, classification)
print(f"  inside the survey: {np.isnan(chm[inside]).mean():.2%} undefined, "
      f"{int((chm < 0).sum())} negative cells")

  CHM 0.0..22.7 m, 10.93% undefined
  inside the survey: 3.67% undefined, 0 negative cells

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.

Reproject, one grid for both surfaces, a supported cell size

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

Counted over the surveyed cells (point density explained runs the sweep):

0.25 m:  76.53% of cells have no ground point
0.50 m:  32.99%
1.00 m:   3.67%
2.00 m:   0.13%

The CHM is limited by the DTM, which is limited by the ground density — 4.2 points per square metre here against 18.3 overall. A 0.5 m CHM would be undefined over 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 CHM with filling, masking and checks

The data used on this page

The USGS 3DEP lidar survey of Brooks Camp, Katmai National Park, Alaska (2012), from the USGS Entwine Point Tile archive on AWS, AK_BrooksCamp_2012: every node under octree node 4-3-13-7 — 287 LAZ files, 12,968,770 points — joined into tile.laz, reprojected from Web Mercator to UTM zone 5N and loaded as arrays. The setup also builds the 1 m lowest-ground-point DTM (dtm, transform), and inside marks the grid cells within the survey's footprint — the patch is turned 2.4° against the UTM grid, so the grid's bounding box has corners that were never surveyed.

import json
import os
import urllib.request

import laspy

EPT = "https://s3-us-west-2.amazonaws.com/usgs-lidar-public/AK_BrooksCamp_2012"
ROOT = (4, 3, 13, 7)                                     # depth, x, y, z of the octree node


def under(key):
    d, x, y, z = map(int, key.split("-"))
    step = 2 ** (d - ROOT[0])
    return d >= ROOT[0] and (x // step, y // step, z // step) == ROOT[1:]


hierarchy = json.load(urllib.request.urlopen(f"{EPT}/ept-hierarchy/0-0-0-0.json"))
nodes = sorted(k for k in hierarchy if under(k))
os.makedirs("nodes", exist_ok=True)
for key in nodes:
    path = f"nodes/{key}.laz"
    if not os.path.exists(path):
        urllib.request.urlretrieve(f"{EPT}/ept-data/{key}.laz", path + ".part")
        os.replace(path + ".part", path)

if not os.path.exists("tile.laz"):                       # the whole subtree as one file
    with laspy.open(f"nodes/{nodes[0]}.laz") as first:
        header = first.header
    with laspy.open("tile.laz", mode="w", header=header) as writer:
        for key in nodes:
            writer.write_points(laspy.read(f"nodes/{key}.laz").points)
with laspy.open("tile.laz") as reader:
    print(len(nodes), "nodes,", f"{reader.header.point_count:,} points,",
          f"{sum(os.path.getsize(f'nodes/{k}.laz') for k in nodes) / 1e6:.1f} MB of nodes")

import numpy as np
from pyproj import Transformer

las = laspy.read("tile.laz")
x, y = Transformer.from_crs(3857, 32605, always_xy=True).transform(las.x, las.y)   # UTM 5N metres
z = np.asarray(las.z)
classification = np.asarray(las.classification)
return_number = np.asarray(las.return_number)
number_of_returns = np.asarray(las.number_of_returns)

import shapely
from scipy.spatial import ConvexHull

hull = ConvexHull(np.c_[x, y])
footprint = shapely.Polygon(np.c_[x, y][hull.vertices])     # the ground the survey covers


def surveyed(transform, shape):
    """True for grid cells whose centre lies inside the survey footprint."""
    rows, cols = np.indices(shape)
    cx, cy = transform * (cols + 0.5, rows + 0.5)
    return shapely.contains_xy(footprint, cx, cy)

from rasterio.transform import from_origin

# the 1 m lowest-ground-point DTM and its ground-point counts, over the bounding box
CELL = 1.0
left, top = x.min(), y.max()
width, height = int(np.ceil((x.max() - left) / CELL)), int(np.ceil((top - y.min()) / CELL))
transform = from_origin(left, top, CELL, CELL)
ground = classification == 2
flat = (np.clip(((top - y[ground]) / CELL).astype(int), 0, height - 1) * width
        + np.clip(((x[ground] - left) / CELL).astype(int), 0, width - 1))
dtm = np.full(width * height, np.inf)
np.minimum.at(dtm, flat, z[ground])
dtm[np.isinf(dtm)] = np.nan
dtm = dtm.reshape(height, width)
counts = np.bincount(flat, minlength=width * height).reshape(height, width)
inside = surveyed(transform, dtm.shape)

287 nodes, 12,968,770 points, 82.8 MB of nodes

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


chm_filled, observed, chm_transform = chm_with_mask(x, y, z, classification, out_path="chm.tif")
print(f"  inside the survey: {np.isnan(chm_filled[inside]).mean():.2%} undefined")

  32 negative cells (0.00%), minimum -0.75 m — inspect the DTM fill
  CHM -0.8..22.7 m, 7.47% undefined, 10.93% of ground interpolated
  inside the survey: 0.00% undefined

Filling the ground turned 32 cells negative, down to −0.75 m — out of about 710,000 — exactly the effect step 5 describes, and small enough to inspect rather than clip. The function's own "undefined" and "interpolated" shares count the whole grid, corners included; inside the survey nothing was left undefined.

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)


p95, cover = canopy_metrics(x, y, z, dtm, transform, 1.0)

  10.0 m cells: median p95 height 10.3 m, median cover 49.6%

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")
    defined_ha = np.isfinite(chm).sum() * cell ** 2 / 1e4     # not the whole grid
    print(f"  density {n / defined_ha:.1f} per hectare")
    return rows, cols, heights


for window in (3.0, 5.0, 9.0):
    tops = detect_trees(chm, window_m=window)

  38397 local maxima above 3.0 m with a 3x3 window (3.0 m)
  heights 3.0..22.7 m, median 10.8 m
  density 561.1 per hectare
  15242 local maxima above 3.0 m with a 5x5 window (5.0 m)
  heights 3.0..22.7 m, median 11.4 m
  density 222.7 per hectare
  5576 local maxima above 3.0 m with a 9x9 window (9.0 m)
  heights 3.4..22.7 m, median 12.2 m
  density 81.5 per hectare

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. Here a 3 m window found 561 tops per hectare, 5 m found 223 and 9 m found 82 — a sevenfold spread from one parameter. Report the window with the tree count, and treat any count without it as meaningless.

The CHM inherits the DTM's holes

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, every surveyed 1 m cell had a return and 3.67% had no ground point — 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 96.3% of surveyed cells 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.

Shared grids and negative values

• 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 often has no returns, so both surfaces are undefined there — though this survey had no such cells.

Internal links

• DSM, DTM and CHM from LiDAR: how each surface is derived — the three surfaces compared
• How to create a DTM from LiDAR ground points in Python — the surface a CHM depends on
• Point density explained — choosing a cell size the ground supports
• How to filter a point cloud by class and return number — selecting returns when classes are missing
• My LiDAR DTM has holes, spikes or terraces — where negative CHM values come from
• How to rasterise a point cloud to a grid in Python — the reduction underneath
• How to measure building heights from LiDAR in Python — the same normalisation, different target
• How to calculate zonal statistics in Python — summarising a CHM by stand or parcel

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.67% of the surveyed cells here. Count only cells inside the survey, or a grid's empty corners look like more holes.

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 3.67% of the surveyed cells undefined and 0.5 m would have left 33%.

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.