🧭 Concept LiDAR & point clouds Terrain & elevation LiDAR & Point Clouds · step 9 of 14 Updated

DSM, DTM and CHM from LiDAR: How Each Surface Is Derived

Problem statement

Three rasters come out of a lidar survey and people use the names interchangeably. They are produced by different reductions of different subsets of the returns, and they fail in different ways.

  • DSM — digital surface model. The top of everything: canopy, roofs, the ground where it is bare.
  • DTM — digital terrain model. The bare earth, from ground-classified returns only.
  • CHM — canopy height model. DSM minus DTM: height above ground.

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

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

The DTM has more holes than the DSM, and the CHM inherits the worse of the two. That asymmetry is the central fact about these products.

Quick answer

Each surface is a different per-cell reduction:

import numpy as np


def cell_reduce(x, y, z, cell, how, bounds):
    """Per-cell min/max/mean of the points falling in each cell."""
    left, bottom, right, top = bounds
    width = int(np.ceil((right - left) / cell))
    height = int(np.ceil((top - bottom) / 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

    if how == "min":
        out = np.full(width * height, np.inf)
        np.minimum.at(out, flat, z); out[np.isinf(out)] = np.nan
    elif how == "max":
        out = np.full(width * height, -np.inf)
        np.maximum.at(out, flat, z); out[np.isinf(out)] = np.nan
    return out.reshape(height, width)


dtm = cell_reduce(x[ground], y[ground], z[ground], 1.0, "min", bounds)
dsm = cell_reduce(x, y, z, 1.0, "max", bounds)
chm = dsm - dtm
A cross-section through canopy showing the DSM as the top of returns, the DTM as the ground surface, and the CHM as the difference.
Three surfaces from one cloud. Only the DTM needs a classification, which is why it has the most holes.

Step-by-step solution

1. DSM: the maximum of all returns

Every pulse contributes, so a DSM is limited only by overall density. In the survey measured here that meant 7.47% empty cells at 1 m — and that 7.47% is water, which returns nothing to a near-infrared laser, not a sampling shortfall.

Use the maximum rather than the mean. A mean over a cell containing both canopy and a gap gives a height that is neither.

2. DTM: the minimum of ground-classified returns

Two constraints, not one. The cell must contain a ground point, and the ground classifier must have found it.

That is why the DTM has 10.93% empty cells against the DSM's 7.47%. The extra 3.46% are cells with returns that contain no ground point — closed canopy, building roofs, and places the classifier declined.

3. Choose min or mean for the DTM deliberately 

min vs mean within the cell: mean |difference| 0.078 m, p99 0.326 m, max 0.94 m

The minimum is standard, and it is the choice that makes you vulnerable to low noise: a single point below the ground drags the whole cell down. The mean is more robust to that and slightly less faithful to the true ground on a slope, where the mean sits above the lowest true point.

At 8 cm typical difference, either is defensible here. On steeper terrain or noisier data the gap widens.

4. CHM: subtract, and expect it to inherit the DTM's holes 

CHM 0.0 .. 22.7 m, 10.93% undefined

The undefined fraction is exactly the DTM's, because a cell without a ground height has no height above ground.

Note the CHM is never negative here. That is a consequence of the construction: the DSM is the maximum over all returns including ground, and the DTM is the minimum over the ground subset, so DSM >= DTM in every cell by definition.

Negative CHM values appear as soon as the DTM is smoothed or interpolated, because then the ground surface in a cell is no longer bounded by that cell's own points. See My LiDAR DTM has holes, spikes or terraces.

5. Normalise the points if you need canopy structure

A CHM is a raster of maximum heights. For canopy metrics — cover, height percentiles, vertical profiles — normalise the point cloud instead:

height_above_ground = z - dtm[row, col]

Then any statistic is available: the 95th percentile of return height, the fraction of returns above 2 m, the number of distinct layers.

A DSM with 7.47 percent empty cells from water against a DTM with 10.93 percent, the extra 3.46 percent being cells with returns but no ground point.
7.47% is genuine absence. The extra 3.46% is the ground classifier finding nothing under canopy and roofs.

Code examples

Example 1 — all three surfaces from one pass 

import numpy as np
import rasterio
from rasterio.transform import from_origin


def lidar_surfaces(x, y, z, classification, cell=1.0, crs="EPSG:32605",
                   ground_class=2):
    """DSM, DTM, CHM and the observation counts behind them."""
    left, bottom = x.min(), y.min()
    right, top = x.max(), y.max()
    width = int(np.ceil((right - left) / cell))
    height = int(np.ceil((top - bottom) / 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 == 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

    n_all = np.zeros(width * height)
    np.add.at(n_all, flat, 1)
    n_ground = np.zeros(width * height)
    np.add.at(n_ground, flat[ground], 1)

    dsm = dsm.reshape(height, width)
    dtm = dtm.reshape(height, width)
    chm = dsm - dtm

    print(f"  {width} x {height} at {cell} m")
    print(f"  DSM {np.nanmin(dsm):7.1f}..{np.nanmax(dsm):7.1f} m, "
          f"{np.isnan(dsm).mean():6.2%} empty")
    print(f"  DTM {np.nanmin(dtm):7.1f}..{np.nanmax(dtm):7.1f} m, "
          f"{np.isnan(dtm).mean():6.2%} empty")
    print(f"  CHM {np.nanmin(chm):7.1f}..{np.nanmax(chm):7.1f} m, "
          f"{np.isnan(chm).mean():6.2%} undefined")

    return {"dsm": dsm, "dtm": dtm, "chm": chm,
            "n_all": n_all.reshape(height, width),
            "n_ground": n_ground.reshape(height, width),
            "transform": transform, "crs": crs}

Returning the counts alongside the surfaces is what makes them interpretable. A DTM cell built from one ground point and one built from forty look identical in the raster.

Example 2 — normalising the cloud, which beats a CHM for structure 

import numpy as np


def normalise_heights(x, y, z, dtm, transform, cell):
    """Height above ground per point, from the DTM cell beneath it."""
    left, top = transform.c, transform.f
    height, width = dtm.shape

    col = np.clip(((x - left) / cell).astype(int), 0, width - 1)
    row = np.clip(((top - y) / cell).astype(int), 0, height - 1)
    ground_z = dtm[row, col]

    normalised = z - ground_z
    supported = np.isfinite(normalised)

    print(f"  {supported.mean():.1%} of points sit over an observed ground cell")
    print(f"  height above ground {np.nanmin(normalised):.2f} .. "
          f"{np.nanmax(normalised):.2f} m")
    below = np.nansum(normalised < -0.5)
    if below:
        print(f"  {int(below):,} points more than 0.5 m below ground "
              "— low noise, or a DTM cell dragged down by a bad point")
    return normalised, supported

The "below ground" count is a free diagnostic. With a per-cell minimum DTM, a point can only be below ground if it is in a different cell from the one that set that minimum — which usually means the ground surface is stepping between cells on a slope.

Example 3 — canopy metrics from normalised heights 

import numpy as np


def canopy_metrics(x, y, height_above_ground, cell=10.0, bounds=None,
                   cover_threshold=2.0):
    """Per-cell canopy statistics that a CHM cannot express."""
    left, bottom, right, top = bounds
    width = int(np.ceil((right - left) / cell))
    height = int(np.ceil((top - bottom) / 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

    order = np.argsort(flat)
    flat_sorted = flat[order]
    h_sorted = height_above_ground[order]
    boundaries = np.searchsorted(flat_sorted, 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 = boundaries[i], boundaries[i + 1]
        if b - a < 5:
            continue
        cell_heights = h_sorted[a:b]
        cell_heights = cell_heights[np.isfinite(cell_heights)]
        if cell_heights.size < 5:
            continue
        p95[i] = np.percentile(cell_heights, 95)
        cover[i] = (cell_heights > cover_threshold).mean()

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

Canopy cover — the fraction of returns above a threshold — is the metric a CHM cannot give, because a CHM has already reduced each cell to one number. Working from the normalised points keeps the vertical distribution available.

Explanation

Why the DTM is the hard one

A DSM asks "what is the highest thing here", which every return answers. A DTM asks "what is the ground here", which only ground returns answer — and identifying them is a classification problem that can fail.

The measured gap is the whole story: 7.47% of cells have no return at all, and 10.93% have no ground return. The extra 3.46% is where the ground exists and was not observed or not recognised.

Under dense canopy that fraction rises sharply, which is why forest DTMs are interpolated over larger gaps than open-ground DTMs and are correspondingly less certain.

Why a per-cell CHM cannot be negative but an interpolated one can

With DSM = max(all returns in cell) and DTM = min(ground returns in cell), the ground returns are a subset of all returns, so the maximum is at least the minimum. The CHM is non-negative by construction — and indeed measured zero negative cells.

The moment the DTM is filled or smoothed, that guarantee is gone. An interpolated ground height in a cell can exceed the highest return in that cell, particularly at a break of slope, and the CHM goes negative.

Negative CHM is therefore a signal about the DTM, not about the canopy. Clipping it to zero hides the diagnostic.

Why maximum for the DSM and minimum for the DTM

Both reductions are chosen to be robust in the direction that matters.

A cell containing canopy and a gap should report the canopy for a DSM, because a DSM is a surface over the top of everything — so maximum. A cell containing ground and low vegetation misclassified as ground should report the lowest, because that is the most likely true ground — so minimum.

The asymmetry has an asymmetric weakness. The DSM is vulnerable to high noise (a bird, a cloud return); the DTM is vulnerable to low noise (a multipath return below the surface). Filtering noise classes before reducing addresses both.

Why cell size is chosen by occupancy, not by density

Nominal spacing suggests a cell size; occupancy decides it. At 16.9 points per square metre the spacing is 0.24 m, yet at 0.25 m cells 78% of cells had no ground point.

The sweep is the tool:

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%

Choose the finest cell size whose empty fraction is close to the floor, and interpolate the remainder — recording that you did.

Minimum and mean per-cell DTM reductions compared on what each resists and what it is vulnerable to.
The minimum is standard because classification errs upward. It pays for that with noise sensitivity.

Edge cases or notes

  • DSM uses all returns; DTM uses ground only. That is why their hole fractions differ.
  • CHM inherits the DTM's holes, always the larger of the two.
  • A per-cell CHM cannot be negative. An interpolated one can, and that is diagnostic.
  • Do not clip negative CHM to zero without investigating the DTM.
  • The minimum reduction is vulnerable to low noise. Filter class 7 and 18 first.
  • Report the observation count per cell alongside every surface.
  • Water returns nothing. Expect permanent holes and do not interpolate across them.
  • Normalise the point cloud for canopy structure; a CHM has already thrown the distribution away.

FAQ

What is the difference between a DSM and a DTM?

A DSM is the top of everything — canopy, roofs, bare ground — from all returns. A DTM is the bare earth, from ground-classified returns only.

What is a canopy height model?

DSM minus DTM: height above ground per cell. It inherits the DTM's holes, so it is undefined wherever no ground point was observed.

Why does my DTM have more holes than my DSM?

Because it needs a ground-classified point in each cell, not just any return. Measured here: 7.47% of cells had no return at all and 10.93% had no ground point.

Should I use the minimum or the mean for a DTM cell?

Minimum is standard and vulnerable to low noise; mean is more robust and slightly high on slopes. They differed by 8 cm on average here.

Why is my CHM negative in places?

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.

What cell size should I use?

The finest whose empty-cell fraction is close to the irreducible floor. On a 16.9 points per square metre survey that was about 1 m for the ground.

How do I get canopy cover rather than canopy height?

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