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, counting only the cells the survey covered (Example 1):

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

The DTM has holes where the DSM has none, and the CHM inherits the DTM's. 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)


ground = classification == 2
bounds = (x.min(), y.min(), x.max(), y.max())
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
print(f"CHM {np.nanmin(chm):.1f} .. {np.nanmax(chm):.1f} m")

CHM 0.0 .. 22.7 m

DSM from all returns, DTM from ground, CHM between

1. DSM: the maximum of all returns

Every pulse contributes, so a DSM is limited only by overall density. In the survey measured here, at 18.3 points per square metre, no 1 m cell inside the survey was empty. (Over the grid's bounding box 7.47% were — the corners outside the survey, which a grid over a rotated patch always has.)

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 3.67% empty cells where the DSM has none. They are cells with returns but no ground point — closed canopy, dense cover, 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, 3.67% of surveyed cells 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.

All three surfaces from one pass

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. surveyed() marks the grid cells inside the survey's footprint: the patch is turned 2.4° against the UTM grid, so a grid over its bounding box has corners that were never surveyed, and they are not counted as holes.

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)

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

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}


surfaces = lidar_surfaces(x, y, z, classification)
inside = surveyed(surfaces["transform"], surfaces["dsm"].shape)
for name in ("dsm", "dtm", "chm"):
    print(f"  {name.upper()} inside the footprint: {np.isnan(surfaces[name][inside]).mean():6.2%} empty")
print(f"  negative CHM cells: {int((surfaces['chm'] < 0).sum())}")

  877 x 876 at 1.0 m
  DSM    16.0..  203.1 m,  7.47% empty
  DTM    15.9..  194.3 m, 10.93% empty
  CHM     0.0..   22.7 m, 10.93% undefined
  DSM inside the footprint:  0.00% empty
  DTM inside the footprint:  3.67% empty
  CHM inside the footprint:  3.67% empty
  negative CHM cells: 0

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


normalised, supported = normalise_heights(x, y, z, surfaces["dtm"], surfaces["transform"], 1.0)

  96.9% of points sit over an observed ground cell
  height above ground -2.02 .. 22.72 m
  1,710 points more than 0.5 m below ground — low noise, or a DTM cell dragged down by a bad point

The "below ground" count is a free diagnostic. With a per-cell minimum DTM no ground point can sit below its own cell's ground, so the 1,710 points more than half a metre under it are other returns — unclassified points lower than every ground return in their cell. They are low noise the vendor did not flag, or ground the classifier missed, and either way they are worth a look before trusting heights near zero.

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)


p95, cover = canopy_metrics(x, y, normalised, bounds=bounds)

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

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.

Why the DTM is the hard one

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: inside the survey no 1 m cell lacks a return, and 3.67% lack a ground return. That 3.67% 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 18.3 points per square metre the spacing is 0.23 m, yet at 0.25 m cells 77% of surveyed cells had no ground point.

The sweep is the tool (point density explained runs it):

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

Choose the finest cell size whose empty fraction you can accept, and interpolate the remainder — recording that you did.

Return selection, holes and negative heights

• 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 often returns nothing. Expect holes there and do not interpolate across them — but count only cells inside the survey before blaming water.
• Normalise the point cloud for canopy structure; a CHM has already thrown the distribution away.

Internal links

• LiDAR point clouds explained: returns, classes and intensity — the returns these surfaces reduce
• Point density explained — choosing the cell size
• How to create a DTM from LiDAR ground points in Python — the implementation
• How to create a canopy height model from LiDAR — the CHM in practice
• My LiDAR DTM has holes, spikes or terraces — filling and its consequences
• How to rasterise a point cloud to a grid in Python — the general reduction
• Digital elevation models explained: DEM, DSM and DTM — the raster products in general
• How to measure building heights from LiDAR in Python — a normalised-height application
• 3D spatial data models explained: 2.5D, meshes and solids — where a surface model sits among the 3D models

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 at 1 m: every surveyed cell had a return, and 3.67% 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 you can accept, counted inside the survey. On an 18.3 points per square metre survey that was about 1 m for the ground: 3.67% empty, and 0.13% at 2 m.

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.