Generating DSM and DTM from Point Clouds with PDAL

Once dense matching finishes, a drone photogrammetry job leaves you with a multi-million-point LAS or LAZ cloud and a deliverable obligation: two rasters that survey clients actually consume. The first is a Digital Surface Model (DSM), the elevation of everything the sensor saw — canopy, rooflines, parked vehicles. The second is a Digital Terrain Model (DTM, also called a DTM/bare-earth or DEM), the ground surface after every above-ground return has been stripped away. Producing both from the same cloud, reproducibly, on a headless machine, is a pipeline problem, and PDAL is the correct tool because a PDAL pipeline is a declarative JSON document: the same JSON yields the same raster on a field laptop and a cloud worker, with no hidden GUI state. This page builds that pipeline end to end and gates the output against CRS, resolution, and NoData-coverage assertions. It sits inside the broader DEM/DSM generation and raster export automation stage, and the GeoTIFFs it emits feed directly into exporting cloud-optimized GeoTIFFs with Rasterio.

Audience and prerequisites. This guide targets Python 3.10+ on a 64-bit OS with at least 16 GB RAM for typical survey-scale clouds. You should understand the difference between a return and a classification, know that LAS classification code 2 means ground, and be comfortable reading a raster’s affine transform and CRS. Ground classification, rasterization interpolation, and void filling are treated as three separable stages so each can be tuned and unit-tested on its own, exactly as the parent raster export automation stage expects.

# Prerequisites

Install the stack in an isolated environment and pin the versions below. The single most common cross-machine discrepancy is a python-pdal binding built against a different PDAL core than the one on PATH; keep the binding and the CLI in the same conda environment so pdal --version and pdal.__version__ agree.

Conda-forge is strongly preferred over pip for PDAL and GDAL because the C++ dependencies (GEOS, PROJ, libLAS) must match exactly; a pip-installed GDAL over a conda PDAL is the usual source of segfaults in writers.gdal. Lock the environment and validate it once against a small reference cloud with known extent before any production run.

# Conceptual architecture

The pipeline forks the same input cloud into two rasterization paths and merges them at export. The DSM path rasterizes all returns with an aggressive aggregate (max or idw) so the highest surface wins. The DTM path first runs a morphological ground filter — filters.smrf (Simple Morphological Filter) or the older filters.pmf (Progressive Morphological Filter) — to stamp classification code 2 onto ground points, keeps only those points with filters.range, then rasterizes them. Both rasters are void-filled to close small NoData gaps, and both are written as GeoTIFFs. Because the ground filter is expensive and the rasterization is cheap, the classification runs once and both the ground-only and all-returns rasters are derived from that single classified stream.

Each stage has a narrow contract — points in, points out for the filters; points in, a raster out for writers.gdal; a raster in, a filled raster out for the void step — so every part can be validated on its own and recombined behind one command-line entry point. The sections below implement them in order.

# Step 1: Define the pipeline JSON

A PDAL pipeline is a JSON object with a single "pipeline" key whose value is an array of stages. The array is read top to bottom: a reader, zero or more filters, and one or more writers. Stages are objects keyed by "type" (readers.las, filters.smrf, writers.gdal), and PDAL infers the connection order from array position, so you rarely need explicit "tag" wiring for a linear flow. Build the JSON as a native Python dictionary rather than a hand-typed string — that way you interpolate filenames and resolution safely and never ship a trailing-comma syntax error to the field. The DSM pipeline below rasterizes every return; keep it deliberately minimal so the interpolation behaviour is easy to reason about.

def build_dsm_pipeline(infile: str, outfile: str, resolution: float = 0.25) -> dict:
    """Return a PDAL pipeline dict that rasterizes ALL returns into a DSM.

    output_type='max' takes the highest point in each cell so rooftops and
    canopy survive; 'idw' would smooth the surface. resolution is the cell
    size in the point cloud's own horizontal units (metres for a UTM cloud).
    """
    return {
        "pipeline": [
            {"type": "readers.las", "filename": infile},
            {
                "type": "writers.gdal",
                "filename": outfile,
                "resolution": resolution,
                "output_type": "max",     # DSM: keep the top surface
                "radius": resolution * 1.5,  # search radius for cell membership
                "gdaldriver": "GTiff",
                "nodata": -9999.0,
            },
        ]
    }

The radius deserves attention now because it decides how large a NoData hole appears where point coverage thins — a value near 1.5 * resolution is a safe default that we tune in the parameter deep-dive. If holes still appear, the dedicated fixing holes and voids in drone-derived DSMs workflow covers both radius tuning and post-hoc fill.

# Step 2: Run the pipeline via pdal.Pipeline in Python

python-pdal executes a pipeline by constructing pdal.Pipeline from the JSON string and calling .execute(), which returns the number of points that passed through and populates .metadata with a JSON tree describing every stage’s output — bounds, spatial reference, point count. Wrap execution so a malformed stage raises a clear, catchable error instead of aborting the batch, and always serialise your dict with json.dumps rather than trusting str(). A pipeline that will not even construct is almost always a JSON-shape or stage-name problem; the fixing PDAL pipeline JSON schema errors page walks through validating it in isolation.

import json
import logging
import pdal

logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")
logger = logging.getLogger(__name__)


def run_pipeline(spec: dict) -> pdal.Pipeline:
    """Execute a PDAL pipeline dict and return the executed Pipeline object."""
    pipeline = pdal.Pipeline(json.dumps(spec))
    try:
        count = pipeline.execute()  # raises RuntimeError on a bad stage
    except RuntimeError as exc:
        logger.error("PDAL pipeline failed: %s", exc)
        raise
    logger.info("Pipeline processed %d points", count)
    return pipeline


dsm = build_dsm_pipeline("cloud.laz", "dsm.tif", resolution=0.25)
run_pipeline(dsm)

.execute() streams points through the whole array in one pass, so peak memory is driven by the widest filter, not by the file size on disk. For clouds that exceed RAM, tile the input first — the streaming LAS tiles to avoid out-of-memory kills approach lets writers.gdal accumulate into one grid across tiles without ever holding the full cloud.

# Step 3: Classify ground with filters.smrf

The DTM path diverges here. filters.smrf implements the Simple Morphological Filter: it estimates a ground surface by progressively opening the point cloud with a growing structuring element and marking points that sit within a slope-and-elevation tolerance of that surface as classification 2. After SMRF stamps the classes, filters.range with "limits": "Classification[2:2]" keeps only ground, and the surviving points feed the same writers.gdal — this time with output_type="idw" so the bare-earth surface is smoothly interpolated rather than spiky. Vegetation that survives the filter is the single most common DTM defect; tuning SMRF against it is covered in removing vegetation from DTM ground surfaces.

def build_dtm_pipeline(infile: str, outfile: str, resolution: float = 0.25) -> dict:
    """Return a PDAL pipeline dict that classifies ground and rasterizes a DTM."""
    return {
        "pipeline": [
            {"type": "readers.las", "filename": infile},
            {
                "type": "filters.smrf",   # Simple Morphological Filter
                "slope": 0.15,            # max terrain slope (rise/run)
                "window": 18.0,           # max structuring-element size, metres
                "threshold": 0.45,        # elevation tolerance, metres
                "scalar": 1.25,           # elevation-scaling factor
            },
            {"type": "filters.range", "limits": "Classification[2:2]"},  # ground only
            {
                "type": "writers.gdal",
                "filename": outfile,
                "resolution": resolution,
                "output_type": "idw",     # DTM: smooth bare-earth interpolation
                "radius": resolution * 2.0,
                "gdaldriver": "GTiff",
                "nodata": -9999.0,
            },
        ]
    }


run_pipeline(build_dtm_pipeline("cloud.laz", "dtm.tif", resolution=0.25))

Prefer filters.smrf over filters.pmf for drone photogrammetry clouds: SMRF handles the sharp breaklines and steep cut slopes typical of construction sites better and needs fewer parameters. Reach for filters.pmf only when you must reproduce a legacy LiDAR workflow that specified it.

# Step 4: Rasterize with writers.gdal

writers.gdal is where points become pixels, and its three decisive knobs are resolution, radius, and output_type. resolution sets the ground sample distance of the output grid in the cloud’s horizontal units. radius is the distance from a cell centre within which points are gathered to compute that cell’s value; too small and cells go NoData, too large and detail smears across the grid. output_type selects the per-cell aggregate. For a DSM you want max (the true top surface) or idw if you prefer a smoother canopy; for a DTM you want idw over the sparse ground points. The inverse-distance weight applied by idw for a cell with contributing points at distances $d_i$ is

$$w_i = \frac{d_i^{-p}}{\sum_j d_j^{-p}}, \qquad z_{\text{cell}} = \sum_i w_i\, z_i$$

with power $p = 2$ by default, so nearer points dominate. Choosing the cell size is a coverage calculation: for a cloud spanning $X$ metres with $N$ ground points, a resolution finer than roughly $\sqrt{X^2 / N}$ guarantees empty cells and a hole-riddled raster, so match resolution to the actual point density rather than to the client’s wish list. You can request multiple aggregates at once by passing output_type as an array (["min", "max", "mean", "idw"]), which writes a multi-band raster — useful when you want to inspect max and idw side by side before committing to one.

# Step 5: Verify CRS, resolution, and NoData coverage with rasterio

A raster that wrote without error is not necessarily a correct deliverable. Open both outputs with rasterio and assert the things that silently break downstream: that the CRS survived from the cloud, that the pixel size matches the requested resolution, and that NoData does not swallow more of the grid than you can tolerate. The NoData fraction is the single most useful health metric — a DTM that is 60% NoData means SMRF rejected most of the scene as non-ground, and a DSM with large NoData patches signals a radius that is too small for the local point density.

import numpy as np
import rasterio

EXPECTED_EPSG = 32633          # UTM 33N — set to your project CRS
EXPECTED_RES = 0.25
MAX_NODATA_FRACTION = 0.05     # fail if >5% of cells are empty


def verify_raster(path: str) -> None:
    with rasterio.open(path) as src:
        assert src.crs is not None, f"{path}: no CRS written"
        assert src.crs.to_epsg() == EXPECTED_EPSG, f"{path}: wrong CRS {src.crs}"
        # Affine.a is x pixel size; b/-e is y pixel size (negative for north-up).
        assert abs(src.res[0] - EXPECTED_RES) < 1e-6, f"{path}: res {src.res}"

        band = src.read(1, masked=True)
        nodata_fraction = band.mask.mean() if np.ma.is_masked(band) else 0.0
        assert nodata_fraction <= MAX_NODATA_FRACTION, (
            f"{path}: {nodata_fraction:.1%} NoData exceeds "
            f"{MAX_NODATA_FRACTION:.0%} budget"
        )
    print(f"{path}: CRS EPSG:{EXPECTED_EPSG}, res {EXPECTED_RES} m, "
          f"NoData {nodata_fraction:.2%}")


verify_raster("dsm.tif")
verify_raster("dtm.tif")

The CRS assertion catches the most damaging silent failure — a raster written with no spatial reference because the source LAS carried no WKT and no override_srs was supplied — which would make the product unplaceable on any map. The NoData budget catches a too-fine resolution or a too-tight ground filter before the raster reaches a client. When the NoData fraction is close to but over budget because of scattered pinholes rather than large voids, fill the raster and re-verify; the export-side compression and overview settings for the final GeoTIFF are covered in setting compression and overviews for orthomosaic COGs.

# Parameter deep-dive

Every parameter trades output fidelity against coverage or runtime. Tune them against a reference cloud with known ground truth, never by guesswork.

# Troubleshooting

pdal.Pipeline raises RuntimeError: Unable to convert... or the process aborts with no raster written. The pipeline constructed but a stage failed at runtime, most often writers.gdal receiving zero points after a too-strict filters.range. Check the point count from .execute(); if it is 0, your Classification[2:2] filter removed everything because SMRF classified no ground. Loosen the SMRF threshold and confirm the reader actually found points.

The DSM has NoData craters over water, glass roofs, and low-texture asphalt. Dense matching produced no points there, so writers.gdal had nothing to interpolate. This is expected physics, not a bug. Widen radius modestly, then fill the residual holes; the full remedy is in fixing holes and voids in drone-derived DSMs.

The DTM is still bumpy with tree and shrub artefacts after ground filtering. SMRF’s window is too small to remove the vegetation footprint, or threshold/scalar are too loose and admit low canopy as ground. Increase window toward the largest non-ground object and tighten threshold; the parameter-by-parameter tuning is in removing vegetation from DTM ground surfaces.

The output raster opens but has no CRS, so it will not overlay in QGIS. The source LAS carried no spatial reference and none was supplied. Add "override_srs": "EPSG:32633" to readers.las, or set "a_srs" on writers.gdal, using the project CRS you validated during coordinate sync.

filters.smrf runs but takes many minutes on a large cloud. SMRF cost scales with the maximum window; a 50 m window on a dense urban cloud is expensive. Tile the cloud, classify per tile, then rasterize into a shared grid, and cap window at the smallest value that still removes your largest building.

Every cell in the DTM equals the NoData sentinel. filters.range kept classification 2, but SMRF never wrote class 2 because it ran after a filter that had already dropped the needed dimensions, or the classes were reset. Ensure filters.smrf precedes filters.range in the array and that no earlier stage strips the Classification dimension.

# Related

• Fixing Holes and Voids in Drone-Derived DSMs
• Removing Vegetation from DTM Ground Surfaces
• Fixing PDAL Pipeline JSON Schema Errors
• Exporting Cloud-Optimized GeoTIFFs with Rasterio

← DEM/DSM Generation & Raster Export Automation