DTM Raster Generation with PDAL

A Digital Terrain Model is the deliverable most LiDAR projects are ultimately commissioned to produce: a continuous bare-earth elevation surface with vegetation, buildings, and vehicles stripped away. Once ground returns have been isolated by an algorithm such as SMRF or PMF, the remaining task is rasterization — collapsing an irregular scatter of classified ground points into a regular grid of elevation values written to a GeoTIFF. PDAL performs this final step with writers.gdal, a rasterizing sink that bins points into cells, interpolates a value per cell, and hands the array to GDAL for encoding. This guide is part of Ground Filtering and DTM/DSM Generation with PDAL, and it focuses on the parameters that decide whether the resulting terrain model is survey-usable or riddled with holes and artefacts.

# Prerequisites

Confirm the following before rasterizing a terrain model:

• PDAL 2.5 or later with Python bindings (pip install pdal) and a working GDAL 3.x underneath.
• A ground-classified point cloud — Classification code 2 already assigned by SMRF Ground Classification or an equivalent filter. Rasterizing an unclassified cloud produces a first-surface DSM, not a DTM.
• A projected, metric CRS. DTM cell sizes only make sense in metres or feet; a cloud still in EPSG:4326 degrees must be reprojected first. See Spatial Reprojection.
• rasterio installed for validation (pip install rasterio), plus numpy.
• Knowledge of your ground-point spacing — the average distance between ground returns sets the finest resolution you can rasterize without holes. The Point Density Metrics guide covers how to measure it.

# Core Workflow Architecture

Producing a DTM with PDAL is a five-phase lifecycle. Each phase maps onto one or two pipeline stages, and the whole sequence executes in a single streaming pass.

1. Ingest and confirm CRS — a readers.las stage loads the classified tile. If the horizontal CRS is geographic or missing, prepend a filters.reprojection stage targeting a projected code such as EPSG:6339 (NAD83(2011) / UTM zone 12N).
2. Select ground returns — a filters.range stage with limits: "Classification[2:2]" discards every non-ground point. This is the single most important step: skip it and every rooftop and tree crown contaminates the surface. The mechanics of this stage are covered in Pipeline Filtering Logic.
3. Rasterize — writers.gdal bins the surviving points into a grid defined by resolution, interpolates a value per cell using output_type, and searches within radius for contributing points.
4. Interpolate small gaps — window_size triggers an in-writer moving-window fill that closes isolated NoData cells from their populated neighbours before the array is encoded.
5. Encode and tag — gdaldriver: "GTiff" with gdalopts sets compression and tiling, nodata marks empty cells, and the CRS carried by the points is written into the GeoTIFF header automatically.

Stage order is not negotiable. filters.range must precede writers.gdal, and any reprojection must precede the range filter so that the classification codes are evaluated on the same points that get rasterized. This is the same buffer-passing discipline described in PDAL Stage Chaining.

# Full Implementation

The module below ingests a classified LAZ tile, keeps ground returns, and writes a 1 m inverse-distance-weighted DTM as a compressed, tiled GeoTIFF. It carries typed signatures, logging, and defensive error handling so it can drop into a batch harness unchanged.

import json
import logging
import pdal
from pathlib import Path

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


def build_dtm_pipeline(
    input_path: str,
    output_tif: str,
    resolution: float = 1.0,
    output_type: str = "idw",
    radius: float | None = None,
    window_size: int = 3,
    nodata: float = -9999.0,
    ground_class: int = 2,
) -> pdal.Pipeline:
    """
    Build a PDAL pipeline that rasterizes ground-classified points into a DTM GeoTIFF.

    Parameters
    ----------
    input_path  : path to a classified LAS/LAZ tile (Classification 2 = ground)
    output_tif  : destination GeoTIFF path
    resolution  : output cell size in ground units (metres for a UTM CRS)
    output_type : cell interpolator — "idw", "mean", "min", "max", "count", "stdev", "all"
    radius      : point search radius; None lets PDAL use resolution * sqrt(2)
    window_size : moving-window radius (in cells) used to fill isolated NoData holes
    nodata      : value written to empty cells
    ground_class : ASPRS classification code to retain (2 = ground)
    """
    gdal_stage: dict = {
        "type": "writers.gdal",
        "filename": str(output_tif),
        "resolution": resolution,
        "output_type": output_type,
        "window_size": window_size,
        "nodata": nodata,
        "gdaldriver": "GTiff",
        "gdalopts": "COMPRESS=DEFLATE,TILED=YES,BLOCKXSIZE=256,BLOCKYSIZE=256,ZLEVEL=6",
        "data_type": "float32",
    }
    if radius is not None:
        gdal_stage["radius"] = radius

    pipeline_def = [
        {"type": "readers.las", "filename": str(input_path)},
        {"type": "filters.range", "limits": f"Classification[{ground_class}:{ground_class}]"},
        gdal_stage,
    ]
    return pdal.Pipeline(json.dumps({"pipeline": pipeline_def}))


def run_dtm(pipeline: pdal.Pipeline) -> int:
    """Validate, execute, and return the number of ground points rasterized."""
    try:
        pipeline.validate()
    except RuntimeError as exc:
        logging.error("DTM pipeline failed validation: %s", exc)
        raise

    count = pipeline.execute()
    if count == 0:
        raise RuntimeError(
            "Zero ground points reached writers.gdal. Check that the input is "
            "classified and that Classification 2 exists in the tile."
        )
    logging.info("Rasterized %d ground points into the DTM.", count)
    return count


if __name__ == "__main__":
    INPUT = Path("ground_classified_utm12n.laz")
    OUTPUT = Path("dtm_1m_idw.tif")

    if not INPUT.exists():
        raise FileNotFoundError(f"Input tile not found: {INPUT}")

    pipe = build_dtm_pipeline(str(INPUT), str(OUTPUT), resolution=1.0, output_type="idw")
    run_dtm(pipe)
    logging.info("DTM written to %s", OUTPUT)

The equivalent declarative JSON, useful for the pdal pipeline CLI or for storing pipelines as version-controlled assets, is:

{
  "pipeline": [
    {
      "type": "readers.las",
      "filename": "ground_classified_utm12n.laz"
    },
    {
      "type": "filters.range",
      "limits": "Classification[2:2]"
    },
    {
      "type": "writers.gdal",
      "filename": "dtm_1m_idw.tif",
      "resolution": 1.0,
      "output_type": "idw",
      "window_size": 3,
      "radius": 1.5,
      "nodata": -9999.0,
      "gdaldriver": "GTiff",
      "gdalopts": "COMPRESS=DEFLATE,TILED=YES,BLOCKXSIZE=256,BLOCKYSIZE=256",
      "data_type": "float32"
    }
  ]
}

# Code Breakdown

# Keeping ground with filters.range

filters.range with limits: "Classification[2:2]" passes only points whose classification equals 2. The bracket syntax is an inclusive range; [2:2] is a single value. If your workflow labels water (class 9) as part of the terrain, widen the expression to "Classification[2:2],Classification[9:9]" — comma-separated ranges are ORed. Because this filter runs before the rasterizer, writers.gdal never sees a rooftop point, so the interpolator cannot pull a cell upward toward a building edge.

# resolution versus radius

resolution fixes the grid: a value of 1.0 in a UTM CRS produces 1 m cells, while 0.5 produces 0.5 m cells and quadruples the pixel count. radius is the search neighbourhood used to gather points for each cell. Left unset, PDAL derives it as resolution * sqrt(2) — just enough to reach the corner of a cell. Explicitly setting radius to roughly 1.5 times the resolution gives the interpolator more points to work with and reduces speckle, at the cost of mild smoothing. Setting it much larger blurs breaklines and drags valley elevations up toward ridge points.

# output_type interpolators

For a production terrain model, idw is the standard. The trade-offs between idw and mean for gap-prone tiles are examined in depth in IDW vs Mean Interpolation for DTM Gaps.

# window_size gap filling

window_size is an integer radius, expressed in cells, for a moving-window interpolation that runs after binning but before encoding. A value of 3 lets PDAL fill an empty cell if it can see populated cells within three cells in any direction. This closes the small, scattered holes that appear at the finest resolutions without a separate post-processing tool. Larger values invent terrain across genuine voids, so keep window_size modest and reach for a dedicated fill only when needed — see Filling NoData Voids in DTM Rasters.

# GDAL driver, options, and NoData

gdaldriver: "GTiff" selects the GeoTIFF encoder. gdalopts is a comma-separated creation-option string passed straight to GDAL: COMPRESS=DEFLATE shrinks the file losslessly, TILED=YES with 256-pixel blocks makes windowed reads fast, and ZLEVEL=6 balances size against write time. nodata: -9999.0 is the sentinel written to every cell no point or interpolation could fill; downstream tools must be told to honour it. data_type: "float32" keeps sub-centimetre vertical precision at half the size of float64.

# Parameter Reference Table

# Validation and Integrity Checks

Never ship a DTM you have not opened. Load the GeoTIFF with rasterio and confirm its shape, CRS, NoData handling, and elevation distribution:

import numpy as np
import rasterio


def validate_dtm(path: str, max_void_fraction: float = 0.05) -> dict:
    """Open a DTM GeoTIFF and assert basic integrity properties."""
    with rasterio.open(path) as ds:
        band = ds.read(1, masked=True)
        total = band.size
        void = int(band.mask.sum())
        void_fraction = void / total

        report = {
            "shape": (ds.height, ds.width),
            "crs": ds.crs.to_string() if ds.crs else None,
            "resolution": ds.res,
            "nodata": ds.nodata,
            "void_fraction": round(void_fraction, 4),
            "z_min": float(band.min()),
            "z_max": float(band.max()),
            "z_mean": float(band.mean()),
        }

    assert report["crs"] is not None, "DTM has no CRS — reproject before rasterizing"
    assert report["nodata"] is not None, "DTM has no NoData value set"
    assert void_fraction <= max_void_fraction, (
        f"{void_fraction:.1%} of cells are NoData (limit {max_void_fraction:.0%}). "
        "Coarsen resolution, widen radius, or run a void-fill pass."
    )
    return report


if __name__ == "__main__":
    print(validate_dtm("dtm_1m_idw.tif"))

Three checks matter most. First, CRS must be present — a null CRS means the source points were untagged, and the raster is spatially meaningless. Second, the NoData fraction should sit below a project threshold; a sudden jump usually means the resolution outran the ground-point spacing. Third, the elevation range should match known site relief — a z_max in the thousands over a floodplain betrays unfiltered noise or a stray high point that ground classification missed.

# Performance Tuning

Rasterization cost is dominated by the output grid size, which scales with the inverse square of resolution. Halving the cell size quadruples the cells and the memory GDAL holds while building the raster. The table below shows representative figures for a 40 million ground-point tile covering roughly 1 km² on a 16-core workstation:

Practical guidance:

• Match resolution to ground-point spacing. If ground returns average 0.7 m apart, a 0.25 m DTM is mostly interpolation and holes; 1 m is honest. Measure spacing with Point Density Metrics before choosing.
• Prefer float32 over float64 unless you genuinely need sub-millimetre vertical precision — it halves both RAM and file size with no visible quality loss for airborne LiDAR.
• Keep TILED=YES so downstream slope, aspect, and hillshade derivations read windows efficiently rather than scanning whole scanlines.
• Tile large areas rather than rasterizing a whole survey at once; a fixed bounds per tile keeps peak memory bounded and lets a batch system parallelise across cores, as covered in Parallel Execution.

# Common Errors and Troubleshooting

writers.gdal: could not create raster / empty output Root cause: no points survived filters.range, so there is nothing to rasterize. Fix: confirm the tile is actually classified — run pdal info --metadata and check for a non-zero count of Classification 2 points before rasterizing.

DTM cell size is one degree, not one metre Root cause: the point cloud is still in a geographic CRS (EPSG:4326), so resolution: 1.0 means one degree. Fix: insert filters.reprojection targeting a projected CRS such as EPSG:6339 before the range filter, as detailed in Spatial Reprojection.

Sawtooth or terraced terrain Root cause: output_type: "min" combined with sparse points quantises the surface, or the vertical values were stored as integers. Fix: switch to idw, set data_type: "float32", and add a small radius so each cell draws from several points.

Speckled NoData holes at fine resolution Root cause: resolution is finer than the ground spacing and window_size is 0. Fix: raise window_size to 3–4 to interpolate across small gaps, widen radius, or coarsen the grid. Persistent large voids need the techniques in Filling NoData Voids in DTM Rasters.

Downstream tools treat -9999 as real elevation Root cause: some consumers ignore the GeoTIFF NoData tag. Fix: confirm nodata is set in the writer and, if a specific tool still mishandles it, mask explicitly with rasterio’s masked reads or re-tag with gdal_edit.py -a_nodata.

# Frequently Asked Questions

What does output_type control in writers.gdal?

output_type selects the statistic each raster cell reports from the ground points that fall inside its search radius. idw computes an inverse-distance-weighted mean and is the standard choice for a smooth bare-earth DTM; mean is an unweighted average; min returns the lowest Z, which suppresses residual low vegetation but is sensitive to negative noise; max, count, and stdev are diagnostic bands. Passing output_type: "all" writes every statistic as a separate band.

How are resolution and radius related in a PDAL DTM?

resolution sets the cell size in ground units and radius sets how far around each cell centre PDAL looks for contributing points. If radius is left unset PDAL defaults it to resolution multiplied by the square root of two, which just covers the cell’s diagonal. Raising radius above that fills more cells at the cost of smoothing and cross-feature bleed; lowering it leaves more NoData holes in sparse areas.

Does writers.gdal need the point cloud to already be in a projected CRS?

It should be. writers.gdal rasterizes in whatever horizontal units the incoming points carry, so a geographic CRS in degrees produces a cell size of, for example, one degree rather than one metre. Reproject to a metric projected CRS with filters.reprojection before rasterizing, and the GeoTIFF inherits that CRS in its header automatically.

Why is my DTM full of NoData cells?

Empty cells appear where no ground point fell within the search radius. Common causes are a resolution finer than the ground-return spacing, an overly tight radius, or aggressive ground filtering that removed too many points. Coarsen the resolution, widen radius, raise window_size to interpolate across small voids, or post-process with a dedicated void-fill pass.

# Related

• Ground Filtering and DTM/DSM Generation with PDAL — parent overview of ground filtering and terrain modelling
• Generating a DTM GeoTIFF with writers.gdal — the complete step-by-step recipe
• IDW vs Mean Interpolation for DTM Gaps — how the two interpolators treat sparse cells
• Filling NoData Voids in DTM Rasters — closing holes with window_size and post-processing
• DSM Generation — the first-return surface counterpart to a DTM
• SMRF Ground Classification — how ground points are labelled before rasterization
• Point Density Metrics — measuring ground spacing to choose a resolution