Generating a High-Resolution DTM from ALS Data: A Production PDAL Workflow

Airborne Laser Scanning (ALS) delivers the vertical structure data that underpins forest ecology, hydrological routing, and slope analysis — but turning a raw point cloud into a sub-meter bare-earth surface is a narrow, repeatable engineering task with its own failure modes. This page walks the exact PDAL-driven pipeline for producing a high-resolution Digital Terrain Model (DTM) from ALS returns: format-safe ingestion, statistical noise suppression, algorithmic ground filtering, and min-return rasterization. It is one concrete recipe within Digital Terrain Model Generation, which sits inside the wider Canopy Height Modeling & Terrain Extraction workflow. The hard parts — leaf-on canopy attenuation, steep gradients, and preserving micro-topographic depressions that drive flow and understory light — are all parameter-calibration problems, and this guide names the parameters that matter.

When to use this approach

ALS-to-DTM is not a one-algorithm problem; the ground filter you pick is dictated by terrain and vegetation structure. Use the table below to choose a starting configuration before you tune.

If you only need height-normalized points rather than a gridded surface, do not rasterize at all — feed classified ground returns into LiDAR Point Cloud Preprocessing and compute height-above-ground there instead. Rasterization is only the right move when a downstream raster operation (hillshade, flow accumulation, CHM subtraction) needs a continuous grid.

Minimal reproducible example

The end-to-end pipeline below ingests a compressed tile, strips noise, classifies ground with CSF, and writes a 0.5 m bare-earth GeoTIFF. It is a single PDAL pipeline executed from Python so the whole sequence stays reproducible and version-controllable.

Step 1 — Bind the lazrs backend and inspect the tile

Reading compressed LAZ with laspy frequently raises ValueError: invalid point format or silently truncates extended classification fields. Bind the lazrs backend explicitly for parallel, specification-compliant decompression, then check what classification you actually received:

import laspy

# Bind lazrs for parallel decompression and extended-attribute support
with laspy.open("forest_tile.laz", laz_backend=laspy.LazBackend.LazrsParallel) as fh:
    las = fh.read()

print(f"Points:  {len(las.points):,}")
print(f"Format:  {las.header.point_format.id}")

# ALS data rarely arrives fully classified; unclassified (1) and
# vegetation (3-5) usually dominate before ground filtering.
import numpy as np
codes, counts = np.unique(las.classification, return_counts=True)
for code, count in zip(codes, counts):
    print(f"  class {code:>2}: {count:,}")

Verify custom classification schemes against the ASPRS LAS Specification before trusting any code values.

Step 2 — Run the full PDAL pipeline

This declarative pipeline chains statistical outlier removal, noise reclassification, CSF ground filtering, and min-return rasterization in one pass:

import json
import pdal

pipeline_def = {
    "pipeline": [
        {"type": "readers.las", "filename": "forest_tile.laz"},
        {
            "type": "filters.outlier",
            "method": "statistical",
            "mean_k": 8,
            "multiplier": 2.5,
        },
        {
            "type": "filters.assign",
            "assignment": "Classification[Classification==7]=18",
        },
        {
            "type": "filters.csf",
            "ignore": "Classification[18:18]",
            "class_threshold": 0.3,
            "cloth_resolution": 1.0,
            "rigidness": 3,
            "iterations": 500,
        },
        {
            "type": "writers.gdal",
            "filename": "dtm_0.5m.tif",
            "resolution": 0.5,
            "output_type": "min",
            "data_type": "float32",
            "where": "Classification == 2",
            "gdalopts": "COMPRESS=DEFLATE,PREDICTOR=2",
        },
    ]
}

pipeline = pdal.Pipeline(json.dumps(pipeline_def))
count = pipeline.execute()
print(f"{count:,} points processed -> dtm_0.5m.tif")

The filters.outlier stage flags statistical outliers as ASPRS code 7, then filters.assign reclassifies them to code 18 (high noise) so the CSF stage can ignore them. The where: "Classification == 2" clause on the writer guarantees only ground-classified returns reach the grid. For karst or terraced sites, swap the filters.csf stage for filters.smrf with slope=0.15, window=18.0, threshold=0.5, and scalar=1.25 — its morphological opening/closing handles abrupt elevation steps that make the cloth model produce false terraces.

Step 3 — Tile large surveys to stay inside memory

A single regional survey will exhaust RAM if loaded whole. Wrap the pipeline in a per-tile function and process tiles independently; PDAL pools memory per execution so workers stay isolated:

import pathlib
from concurrent.futures import ProcessPoolExecutor


def build_dtm(tile: pathlib.Path) -> str:
    out = pathlib.Path("dtm") / f"{tile.stem}.tif"
    defn = json.loads(json.dumps(pipeline_def))           # deep copy
    defn["pipeline"][0]["filename"] = str(tile)
    defn["pipeline"][-1]["filename"] = str(out)
    pdal.Pipeline(json.dumps(defn)).execute()
    return out.name


tiles = sorted(pathlib.Path("raw_tiles").glob("*.laz"))
with ProcessPoolExecutor(max_workers=4) as pool:
    for name in pool.map(build_dtm, tiles):
        print(f"wrote {name}")

Mosaic the per-tile rasters afterward with gdal_merge.py or rasterio.merge, keeping the cell size identical so the tiles align on the same grid.

Parameter reference

output_type selects PDAL’s per-cell aggregation (min, max, mean, idw, count, stdev), while data_type sets the on-disk pixel dtype — conflating the two is a common copy-paste bug that only surfaces when GDAL refuses to open the result.

Expected output and verification

A correct run produces a single-band float32 GeoTIFF at the requested cell size, populated everywhere ground returns exist and NaN/nodata only in genuine voids (water bodies, building footprints, deep canopy with no ground penetration). Validate it before trusting it downstream:

import numpy as np
import rasterio

with rasterio.open("dtm_0.5m.tif") as src:
    dtm = src.read(1, masked=True)
    assert src.res == (0.5, 0.5), f"unexpected cell size {src.res}"
    assert src.crs is not None, "DTM is missing a CRS — check input header"

    # A clean bare-earth tile should be mostly filled; large void
    # fractions usually mean the ground filter was too aggressive.
    void_fraction = dtm.mask.mean()
    print(f"void fraction: {void_fraction:.1%}")
    assert void_fraction < 0.15, "too many voids — relax class_threshold"

    print(f"elevation range: {dtm.min():.1f}-{dtm.max():.1f} m")

For a visual check, derive a hillshade and overlay the classified ground points: systematic voids or a terraced look in the hillshade signal misclassification. Cross-reference the surface against USGS 3D Elevation Program benchmarks where coverage exists, and confirm hydrological plausibility by running flow accumulation with richdem or whitebox — micro-topographic depressions like vernal pools must survive into the final grid.

Common pitfalls

• Min-return picks up low noise. output_type="min" is correct for bare earth but will latch onto a single sub-ground noise return if outlier removal was skipped — always run filters.outlier first and reclassify noise to code 18 so it is excluded.
• CSF sags into leaf-on understory. Default class_threshold=0.5 classifies dense shrub returns as ground in leaf-on broadleaf stands, inflating the DTM by tens of centimetres. Drop to 0.25–0.35 and coarsen cloth_resolution to 1.0.
• Cell size finer than pulse density. Requesting a 0.25 m grid from a 4 pt/m² survey produces empty cells that interpolation must guess — match resolution to true density to avoid fabricated micro-topography.
• output_type vs data_type swap. Writing output_type="float32" (invalid) or data_type="min" (invalid) yields a file GDAL silently mishandles or refuses to open — keep aggregation and pixel dtype distinct.

Frequently Asked Questions

Related

• Digital Terrain Model Generation — the parent workflow this recipe belongs to
• Canopy Height Modeling & Terrain Extraction — the full point-cloud-to-canopy pipeline
• LiDAR Point Cloud Preprocessing — noise filtering and classification upstream of this step
• Normalizing LiDAR Point Clouds with PDAL — compute height-above-ground from the DTM
• Canopy Height Model Creation — subtract this DTM from a DSM to derive the CHM