Python LiDAR

Building a Normalized DSM

TL;DR: nDSM = DSM − DTM on identical grids. Either write both rasters with the same origin_x, origin_y, width, height and resolution in writers.gdal and subtract them with rasterio, or skip the subtraction by computing HeightAboveGround with filters.hag_nn and rasterizing its per-cell maximum. Clamp small negative values to zero and keep NoData where either input is missing.

# Context and Motivation

This guide is part of DSM Generation. A DSM describes the top of everything — roofs, canopy, bridges — in absolute elevation, which makes it awkward for asking “how tall is that?”. A 12 m building on a hill and a 12 m building in a valley have very different DSM values. Subtracting the terrain gives heights above ground everywhere, and that normalized surface is what building-height maps, canopy analysis, urban morphology and change detection actually need. Over vegetation, the nDSM is essentially a canopy height model; over built areas, it is a building-height map.

Getting it right is mostly about alignment: two rasters that differ by half a cell produce ghost edges along every wall.

DSM − DTM = heights Three stacked profiles. The DSM traces roofs and tree tops on a sloping hillside. The DTM traces the sloping ground alone. Their difference, the nDSM, shows the building and tree as heights on a flat zero baseline, so a 12 metre building reads as 12 metres whatever the terrain beneath it. DSM DTM nDSM 12 m 9 m zero baseline everywhere

# Prerequisites and Assumptions

  • A classified point cloud, or an existing DSM and DTM of the same area.
  • PDAL, rasterio and NumPy.
  • The same CRS for both rasters, and both built at the same resolution.

# Step-by-Step Implementation

# Step 1 — Define one grid

Compute an origin and size from the tile bounds, rounded to the resolution, and use them for every raster.

# Step 2 — Write DSM and DTM on that grid

DSM: maximum Z of first returns. DTM: IDW of ground returns. Both with the shared origin_x, origin_y, width, height.

# Step 3 — Subtract

Read both with rasterio, subtract where both are valid, and set NoData elsewhere.

# Step 4 — Clean small negatives

Values slightly below zero come from interpolation differences between the two surfaces; clamp values between −0.5 m and 0 to zero, and treat larger negatives as errors to inspect.

# Step 5 — Alternative: rasterize HAG directly

filters.hag_nn followed by writers.gdal with dimension: "HeightAboveGround" and output_type: "max" yields an nDSM in one pipeline, without a separate DTM.

# Complete Working Example

python
"""nDSM from aligned DSM and DTM rasters, plus the one-pipeline HAG alternative."""
from __future__ import annotations

import json
import math

import numpy as np
import pdal
import rasterio

SRC, RES = "tiles/t_0431.laz", 0.5


def grid(src: str, res: float) -> dict:
    b = pdal.Pipeline(json.dumps({"pipeline": [src]})).quickinfo["readers.las"]["bounds"]
    ox, oy = math.floor(b["minx"] / res) * res, math.floor(b["miny"] / res) * res
    return {"origin_x": ox, "origin_y": oy,
            "width": math.ceil((b["maxx"] - ox) / res), "height": math.ceil((b["maxy"] - oy) / res)}


def rasters(src: str, res: float) -> None:
    g = grid(src, res)
    common = {"resolution": res, "data_type": "float32", "nodata": -9999, **g}
    pdal.Pipeline(json.dumps({"pipeline": [
        src,
        {"type": "filters.range", "limits": "Classification![7:7],Classification![18:18]", "tag": "clean"},
        {"type": "filters.range", "inputs": ["clean"], "limits": "ReturnNumber[1:1]", "tag": "first"},
        {"type": "writers.gdal", "inputs": ["first"], "filename": "out/dsm.tif",
         "output_type": "max", "radius": res * 0.71, **common},
        {"type": "filters.range", "inputs": ["clean"], "limits": "Classification[2:2]", "tag": "ground"},
        {"type": "writers.gdal", "inputs": ["ground"], "filename": "out/dtm.tif",
         "output_type": "idw", "window_size": 6, **common},
    ]})).execute()


def ndsm(dsm: str, dtm: str, dst: str) -> None:
    with rasterio.open(dsm) as a, rasterio.open(dtm) as b:
        assert a.transform == b.transform and a.shape == b.shape, "grids differ"
        s, t = a.read(1, masked=True), b.read(1, masked=True)
        profile = a.profile
    h = s - t
    h = np.ma.where((h < 0) & (h > -0.5), 0.0, h)
    with rasterio.open(dst, "w", **profile) as out:
        out.write(h.filled(-9999).astype("float32"), 1)
    print(f"nDSM: max {h.max():.1f} m, cells < −0.5 m: {(h < -0.5).sum()}")


def ndsm_from_hag(src: str, dst: str, res: float) -> None:
    pdal.Pipeline(json.dumps({"pipeline": [
        src,
        {"type": "filters.range", "limits": "Classification![7:7],Classification![18:18]"},
        {"type": "filters.hag_nn", "count": 2},
        {"type": "writers.gdal", "filename": dst, "dimension": "HeightAboveGround",
         "output_type": "max", "resolution": res, "radius": res * 0.71,
         "data_type": "float32", "nodata": -9999, **grid(src, res)},
    ]})).execute()


if __name__ == "__main__":
    rasters(SRC, RES)
    ndsm("out/dsm.tif", "out/dtm.tif", "out/ndsm.tif")
    ndsm_from_hag(SRC, "out/ndsm_hag.tif", RES)
Half a cell of misalignment Left: DSM and DTM grids offset by half a cell. Subtracting them produces thin false spikes and trenches along a building's walls, visible as bright and dark lines in the nDSM. Right: with a shared origin and size, the subtraction is clean and the building has sharp, correct edges. grids offset by ½ cell shared origin and size ghost spikes and trenches at walls clean edges

# Key Parameter Table

Setting DSM DTM Why
input first returns class 2 Top surface versus terrain
output_type max idw Highest return; smooth terrain
radius 0.71 × res — Covers each cell
window_size — 6 Fills small ground gaps
grid shared shared Identical origin, width, height
nodata −9999 −9999 Masked in the subtraction

# Verification

  • Grids identical. The transform and shape assertion in ndsm must pass.
  • Ground near zero. On open ground and roads, nDSM values should be within ±0.2 m of zero.
  • Known heights. A few buildings of known height should read correctly to within a few decimetres.
  • Two methods agree. The subtraction and the HAG route should agree closely except at object edges.

# Gotchas and Edge Cases

Water and bridges. Bridges appear as tall objects over rivers; water surfaces may show small positive or negative values from sparse returns. Mask both if the nDSM feeds building or vegetation statistics.

DTM errors become heights. Where ground classification failed — a building accepted as ground — the nDSM shows zero height there. nDSM quality is bounded by DTM quality.

Pits in the DSM. Laser pulses that penetrate gaps in canopy leave low DSM cells, producing pits in the nDSM. Clean the DSM first; see removing pits and spikes from a DSM.

What an nDSM histogram shows A histogram of nDSM cell values on a log count axis. A tall spike at zero represents open ground and roads. A broad hump from 2 to 25 metres represents vegetation. A smaller cluster between 3 and 12 metres represents buildings. A thin tail below zero represents interpolation mismatch and should stay small. 0 m: ground buildings vegetation negatives: keep small

Resolution mismatch. Resampling a 1 m DTM to 0.5 m before subtraction is fine; resampling a 0.5 m DSM down to 1 m loses roof edges. Build both at the target resolution where possible.

# Frequently Asked Questions

What is a normalized DSM?

A raster of heights above the terrain, computed as the digital surface model minus the digital terrain model. Buildings and trees appear with their actual heights regardless of the ground elevation beneath them.

How do I make sure the DSM and DTM align?

Write both with writers.gdal using the same origin_x, origin_y, width, height and resolution, and assert that the transforms and shapes match before subtracting.

Why does my nDSM have negative values?

Small negatives come from differences in how the two surfaces are interpolated and are usually clamped to zero. Large negatives indicate a problem, such as a DSM pit or a ground point misclassified above the true terrain.

Is an nDSM the same as a canopy height model?

Over vegetation, effectively yes. A canopy height model is typically built from vegetation returns only, while an nDSM includes buildings and other objects too.

Related in DSM Generation from LiDAR with PDAL

↑ DSM Generation from LiDAR with PDAL overview --- title: "DSM Generation from LiDAR with PDAL" description: "Building a Digital Surface Model from LiDAR first returns with PDAL... Removing Pits and Spikes from a DSM --- title: "Removing Pits and Spikes from a DSM" description: "Clean a LiDAR digital surface model of data pits (laser penetration into... Building a DSM from First Returns --- title: "Building a DSM from First Returns" description: "A PDAL recipe to build a Digital Surface Model from first-return LiDAR points... DTM vs DSM: Which Surface Model to Generate --- title: "DTM vs DSM: Which Surface Model to Generate" description: "When to build a bare-earth DTM versus a first-return DSM from LiDAR —...