DEM/DSM Generation & Raster Export Automation

The final stage of a drone photogrammetry pipeline is where geometry becomes a deliverable: a dense point cloud or reconstructed mesh has to be resampled into survey-grade raster surfaces that a GIS, a hydrologist, or a stockpile-volume report can consume directly. For UAV operators, surveying technicians, and Python GIS developers the hard engineering problem here is not calling one tool once — it is producing a Digital Surface Model (DSM), a Digital Terrain Model (DTM, the bare-earth surface), and their derivatives deterministically: same point cloud, same configuration, same bytes out, with an explicit horizontal and vertical coordinate reference system (CRS) baked into every file. A DEM that is internally smooth but sits 30 m too high because someone rasterized ellipsoidal heights against an orthometric expectation is worse than no DEM at all, because it looks correct.

Python is the right automation layer for this because the entire raster toolchain — PDAL, rasterio, GDAL, and NumPy — exposes first-class bindings that let a single orchestrator classify ground points, rasterize surfaces, fill voids, compress and tile Cloud-Optimized GeoTIFFs (COGs), and derive contours and hillshades without shelling out to disconnected GUI applications. Each stage validates its inputs, pins its CRS and NoData contract, and writes an inspectable artifact. The dense point clouds consumed here are produced upstream by the automated image alignment and feature matching workflows, and the horizontal and vertical datum they carry is established during ground control point optimization and coordinate sync; this page assumes both are already correct and focuses on turning classified points into validated raster products.

Figure 1 — The DSM (first-return surface) and DTM (bare-earth surface) are rasterized from the same classified cloud but with different point selections, then share the void-filling, COG-export, and derivative stages that turn them into validated deliverables.

The four stages below map to the deep-dive guides in this section: point classification and rasterization with PDAL, COG export with rasterio, contour and hillshade derivation with GDAL, and the validation discipline that catches broken exports before they ship. Read this page for the architecture and the contracts between stages; follow each inline link when you need the full runnable implementation.

# Ground and Surface Classification

Every raster surface a drone survey produces is a selection of points followed by an interpolation. A DSM keeps the highest return in each neighbourhood — rooftops, canopy, parked vehicles — because it models the visible surface. A DTM keeps only points classified as bare ground and discards everything that sits on top of it. Getting these two surfaces right therefore begins with classification, not rasterization: if vegetation and building points leak into the ground class, the DTM will show phantom mounds where a hedge used to be, and no amount of downstream smoothing recovers the true terrain.

PDAL’s Simple Morphological Filter (filters.smrf) and Progressive Morphological Filter (filters.pmf) are the workhorses here. SMRF slides a morphological opening across a provisional surface grid, flagging any point that rises more than threshold metres above the opened surface as non-ground; the window, slope, and scalar parameters control how aggressively it climbs over real terrain relief without swallowing low vegetation. Because the classification result is the contract that the rasterization stage depends on, it must be tuned per scene — steep vineyards, flat quarries, and dense urban blocks each want different windows — and the full parameter treatment lives in generating DSM and DTM from point clouds with PDAL.

import json
import pdal

def classify_and_rasterize_dtm(las_path: str, out_tif: str,
                               resolution: float = 0.25,
                               target_crs: str = "EPSG:32633+5703") -> int:
    """Classify ground with SMRF, keep class 2, and rasterize a bare-earth DTM.

    target_crs is a compound code: horizontal (UTM 33N) + vertical (EGM2008),
    so the DTM carries an explicit orthometric datum, not raw ellipsoidal height.
    """
    pipeline = {
        "pipeline": [
            {"type": "readers.las", "filename": las_path},
            {"type": "filters.reprojection", "out_srs": target_crs},
            # SMRF ground classification -> writes Classification == 2 for ground
            {"type": "filters.smrf", "scalar": 1.2, "slope": 0.15,
             "threshold": 0.45, "window": 18.0},
            # Keep only ground returns for the bare-earth surface
            {"type": "filters.range", "limits": "Classification[2:2]"},
            {"type": "writers.gdal", "filename": out_tif,
             "resolution": resolution, "radius": resolution * 1.5,
             "output_type": "idw", "gdaldriver": "GTiff",
             "nodata": -9999, "data_type": "float32"},
        ]
    }
    p = pdal.Pipeline(json.dumps(pipeline))
    n_points = p.execute()   # returns the number of points that passed the pipeline
    return n_points

The contract this stage emits is precise: a classified cloud in a known compound CRS where Classification == 2 marks bare ground and every point carries a ReturnNumber the DSM stage can filter on. Validate that contract at the boundary — assert that the ground class is non-empty and that the reprojected envelope falls inside the project’s expected bounds — because a malformed filters.smrf result becomes a wrong DTM two stages later, far from its cause. If the pipeline JSON itself is rejected, the schema-level fixes are catalogued in fixing PDAL pipeline JSON schema errors, and residual vegetation that survives classification is addressed in removing vegetation from DTM ground surfaces.

# DSM/DTM Rasterization and Void Filling

Once points are classified, rasterization maps them onto a regular grid. PDAL’s writers.gdal bins points into cells and computes one or more statistics per cell — min, max, mean, idw, or count. A DSM typically uses max (or first-return filtering plus max) so tall features are preserved; a DTM uses idw (inverse-distance weighting) or mean over ground points to produce a smooth continuous surface. The two parameters that govern quality are resolution (the cell size in CRS units) and radius (how far writers.gdal searches for points to populate a cell). Set resolution finer than your point spacing and cells will be starved of data; set radius too large and you smear detail across gaps. A defensible starting point ties the cell size to point density $\rho$ (points per m²):

$$\text{resolution} \gtrsim \frac{1}{\sqrt{\rho}}, \qquad \text{radius} \approx 1.5 \times \text{resolution}$$

Even a well-tuned grid leaves NoData cells wherever occlusion, water absorption, or sparse returns starved the interpolator. These voids must be filled deterministically before export, because most consumers (volume calculations, contouring, flood models) cannot tolerate holes. rasterio.fill.fillnodata runs an inverse-distance interpolation across the NoData mask with a bounded search distance, which keeps the fill local and reproducible rather than hallucinating terrain across a wide river. The void-filling parameters — and when to prefer a coarser fallback grid over aggressive filling — are covered in fixing holes and voids in drone-derived DSMs.

import numpy as np
import rasterio
from rasterio.fill import fillnodata

def fill_dem_voids(src_path: str, dst_path: str,
                   max_search: int = 100, smoothing: int = 0) -> float:
    """Fill NoData voids with bounded inverse-distance interpolation.

    Returns the void fraction BEFORE filling so the caller can gate on it:
    a surface that was 40% holes should be re-gridded, not patched.
    """
    with rasterio.open(src_path) as src:
        band = src.read(1)
        profile = src.profile.copy()
        nodata = src.nodata if src.nodata is not None else -9999.0

    valid_mask = (band != nodata).astype("uint8")   # 1 = keep, 0 = interpolate
    void_fraction = 1.0 - valid_mask.mean()

    filled = fillnodata(band, mask=valid_mask,
                        max_search_distance=max_search,
                        smoothing_iterations=smoothing)

    profile.update(dtype="float32", nodata=nodata)
    with rasterio.open(dst_path, "w", **profile) as dst:
        dst.write(filled.astype("float32"), 1)
    return void_fraction

Returning the pre-fill void fraction is deliberate: it turns a silent quality problem into a value the orchestrator can threshold on. A 3% void fraction is routine; a 40% void fraction means the grid resolution is finer than the data supports, and patching it produces a smooth-looking lie. Gate on the fraction, do not just fill and forget.

# Cloud-Optimized GeoTIFF Export

A raster that will live in object storage and be read by web maps and remote clients should be a Cloud-Optimized GeoTIFF: internally tiled, with pre-computed overview pyramids and a header laid out so an HTTP range request can fetch one tile without downloading the whole file. Producing a valid COG deterministically is where most export pipelines quietly fail — a file can be a perfectly good GeoTIFF and still fail COG validation because its overviews are missing or its blocks are striped rather than tiled. rasterio drives GDAL’s dedicated COG driver, which enforces the layout rules for you, and the choice between rasterio and raw GDAL for this step is weighed in GDAL vs Rasterio for DEM and orthomosaic export.

import rasterio
from rasterio.shutil import copy as rio_copy

def write_cog(src_path: str, cog_path: str,
              blocksize: int = 512,
              compress: str = "DEFLATE",
              predictor: int = 3) -> None:
    """Rewrite a plain GeoTIFF as a validated Cloud-Optimized GeoTIFF.

    predictor=3 is the floating-point predictor, correct for continuous DEMs;
    use predictor=2 (horizontal) for integer orthomosaic bands instead.
    """
    cog_options = {
        "driver": "COG",
        "COMPRESS": compress,          # DEFLATE / ZSTD (lossless) or LZW
        "PREDICTOR": predictor,        # 3 = float predictor for elevation rasters
        "BLOCKSIZE": blocksize,        # internal tile size in pixels
        "OVERVIEWS": "AUTO",           # build overview pyramid automatically
        "OVERVIEW_RESAMPLING": "AVERAGE",
        "BIGTIFF": "IF_SAFER",         # switch to BigTIFF past the 4 GB limit
    }
    with rasterio.open(src_path) as src:
        rio_copy(src, cog_path, **cog_options)

The two compression knobs matter for DEMs specifically. COMPRESS=DEFLATE (or ZSTD) is lossless, which is mandatory for elevation data — never ship a lossy-JPEG DEM. The PREDICTOR=3 floating-point predictor exploits the fact that neighbouring elevation values differ by small amounts, shrinking file size substantially with no accuracy cost; using the wrong predictor for the data type either bloats the file or corrupts it. The full matrix of compression, predictor, blocksize, and overview-resampling choices for both DEMs and RGB orthomosaics is worked through in exporting Cloud-Optimized GeoTIFFs with rasterio and, for the orthomosaic-specific tuning, in setting compression and overviews for orthomosaic COGs.

Export is also where the CRS and vertical-datum discipline established during ground control point optimization has to be re-asserted. GDAL will happily write a file whose horizontal CRS is UTM but whose vertical datum tag is absent, leaving the consumer to guess whether elevations are ellipsoidal or orthometric. Write a compound CRS (for example EPSG:32633+5703) into the output profile so the vertical reference travels with the pixels, not in a README nobody reads.

# Contour and Hillshade Derivation

The DEM is rarely the final human-facing product; contours and hillshades are what land in a plan set or a briefing map. GDAL’s DEMProcessing computes hillshade, slope, aspect, roughness, and colour relief from a DEM in a single call, and ContourGenerate walks the grid to emit vector isolines at a fixed interval. Deriving these from the DEM rather than re-processing the point cloud keeps them perfectly registered to the surface they describe — a contour that does not lie on its own DEM is a data-integrity bug waiting to be noticed on site.

Hillshade is an analytic function of the local surface gradient and a fixed illumination direction. For a solar zenith $z_s$, solar azimuth $a_s$, cell slope $S$, and cell aspect $A$, the shaded value is

$$H = 255\,\bigl(\cos z_s \cos S + \sin z_s \sin S \cos(a_s - A)\bigr),$$

which is why the azimuth and altitude parameters, not the DEM, decide where shadows fall — a point worth remembering when two hillshades of the same site look inconsistent.

from osgeo import gdal, ogr

gdal.UseExceptions()

def derive_hillshade(dem_path: str, hs_path: str,
                     z_factor: float = 1.0,
                     azimuth: float = 315.0, altitude: float = 45.0) -> None:
    """Standard north-west illuminated hillshade with edge handling."""
    gdal.DEMProcessing(
        hs_path, dem_path, "hillshade",
        options=gdal.DEMProcessingOptions(
            zFactor=z_factor, azimuth=azimuth, altitude=altitude,
            computeEdges=True, format="GTiff",
        ),
    )

def derive_contours(dem_path: str, gpkg_path: str, interval: float = 1.0) -> None:
    """Generate elevation contours at a fixed interval into a GeoPackage."""
    src = gdal.Open(dem_path)
    band = src.GetRasterBand(1)
    drv = gdal.GetDriverByName("GPKG")
    dst = drv.Create(gpkg_path, 0, 0, 0, gdal.GDT_Unknown)
    layer = dst.CreateLayer("contours", srs=None, geom_type=ogr.wkbLineString)
    layer.CreateField(ogr.FieldDefn("elev", ogr.OFTReal))
    # band, interval, base, fixedLevels, useNoData, noDataValue, layer, idField, elevField
    gdal.ContourGenerate(band, interval, 0.0, [], 0, 0.0, layer, -1, 0)
    dst = None   # flush and close the GeoPackage

Two choices dominate derivative quality. z_factor must be 1.0 when horizontal and vertical units match (both metres); it exists precisely to correct the pathological case of a DEM in geographic degrees with elevation in metres, where an uncorrected hillshade renders as noise. The contour interval should match the vertical accuracy of the survey — emitting 0.1 m contours from a DEM whose vertical RMSE is 0.15 m manufactures precision that is not in the data. Tooling trade-offs between GDAL’s command-line utilities and the rasterio/NumPy path for these derivatives are compared in GDAL vs Rasterio for DEM and orthomosaic export.

# Stage Parameter Reference

These are the thresholds, driver options, and datum settings that govern behaviour across the four stages. Treat the defaults as starting points to tune per platform and per scene, not universal constants.

# Failure Modes and Diagnostics

Raster-export failures are usually deterministic and detectable in Python well before a human notices a warped surface or a black tile. The patterns below recur across drone DEM jobs; each has a programmatic symptom and a concrete remediation, and the full runbook is troubleshooting DEM and raster export failures.

• Excessive NoData voids. Symptom: the exported DEM is riddled with holes and void_fraction exceeds your gate. Cause: the grid resolution is finer than the point spacing, so cells receive no returns. Detect it by measuring the void fraction directly and re-grid at a coarser cell size rather than mass-filling.

import numpy as np
import rasterio

def void_fraction(path: str, fallback_nodata: float = -9999.0) -> float:
    """Fraction of NoData cells in a single-band raster."""
    with rasterio.open(path) as src:
        arr = src.read(1)
        nd = src.nodata if src.nodata is not None else fallback_nodata
    return float(np.count_nonzero(arr == nd)) / arr.size

• Silent vertical-datum offset. Symptom: the DEM is internally smooth but its mean elevation is off by an amount close to the local geoid separation (tens of metres in many regions). Cause: ellipsoidal heights were rasterized where the deliverable expected orthometric (geoid-referenced) height. Guard by asserting the mean elevation falls inside the site’s known orthometric range and by writing a compound CRS on export.

def assert_orthometric(mean_elev: float, expected: tuple[float, float]) -> None:
    """Fail fast if the DEM mean elevation is outside the site's orthometric range."""
    lo, hi = expected
    if not (lo <= mean_elev <= hi):
        raise ValueError(
            f"Mean elevation {mean_elev:.1f} m outside expected {expected} m; "
            "check for an ellipsoidal vs geoid (vertical datum) mismatch."
        )

• Invalid COG. Symptom: rio_cogeo.cog_validate returns False or a downstream tile server cannot range-read the file. Cause: overviews were never built, or the file is striped rather than internally tiled. Remediate by re-exporting through the COG driver and validating in CI. The specific validator messages are decoded in fixing COG validation errors in GDAL.

from rio_cogeo.cog_validate import cog_validate

is_valid, errors, warnings = cog_validate(cog_path)
if not is_valid:
    raise RuntimeError(f"COG validation failed: {errors}")

• Terrain spikes in the DTM. Symptom: a supposedly bare-earth surface shows high local slope variance where vegetation stood. Cause: filters.smrf left canopy or building points in the ground class. Remediate by tightening smrf.threshold/smrf.window or adding a height-above-ground filter, as detailed in removing vegetation from DTM ground surfaces.

• Cannot allocate memory on export. Symptom: GDAL or rasterio aborts while writing a large raster or building overviews. Cause: a whole-array read on a survey-scale DEM, or BIGTIFF not enabled past 4 GB. Remediate with windowed reads and BIGTIFF="IF_SAFER", per decoding GDAL “Cannot Allocate Memory” on export.

• Empty or black tiles. Symptom: the orthomosaic or DEM COG opens but shows all-NoData or all-black tiles at some zoom levels. Cause: a NoData sentinel mismatch between rasterization and export, or overviews resampled with the wrong method. Diagnose per diagnosing empty or black orthomosaic tiles.

# Integration Checklist

Wire the stages together for a production run by confirming each contract below. These render as interactive toggles.

• Reproject the input cloud to an explicit compound CRS (horizontal + vertical) before classification, never after export
• Run filters.smrf/filters.pmf and assert the ground class (Classification == 2) is non-empty
• Rasterize the DSM from first returns / max and the DTM from ground-only points with idw/mean
• Set resolution no finer than $1/\sqrt{\rho}$ and radius to roughly 1.5× resolution
• Measure the pre-fill void fraction and gate on it; re-grid coarser instead of patching voids above the limit
• Fill remaining voids with a bounded max_search_distance so wide gaps stay NoData
• Export through the COG driver with a lossless codec (DEFLATE/ZSTD) and PREDICTOR=3 for elevation
• Build overview pyramids with AVERAGE resampling for continuous DEMs
• Write the compound CRS and a consistent nodata sentinel into every output profile
• Validate every COG with cog_validate in CI and fail the build on errors
• Derive contours and hillshades from the final DEM, matching contour interval to the vertical RMSE
• Assert the DEM mean elevation against the site’s orthometric range and log the CRS, resolution, and point count with each run

By enforcing these contracts, UAV operators and GIS developers can scale DEM and raster export from a single quarry survey to a fleet of recurring monitoring flights without shipping surfaces that are smooth, tiled, and quietly wrong.

# Related

• Generating DSM and DTM from Point Clouds with PDAL
• Exporting Cloud-Optimized GeoTIFFs with Rasterio
• GDAL vs Rasterio for DEM and Orthomosaic Export
• Troubleshooting DEM and Raster Export Failures
• Automated Image Alignment & Feature Matching Workflows

← Drone Photogrammetry Pipelines