Create a DEM from a Point Cloud in PyQGIS

Turning a cloud into a raster is the step where LiDAR stops being LiDAR and becomes a layer the rest of your workflow can use. It is also the step where most of the mistakes get baked in permanently, because a DEM at the wrong resolution, or filled by the wrong method, looks entirely credible. Nobody opens a hillshade and thinks "those terraces are interpolation artefacts"; they think the hillside has terraces.

This recipe belongs to Point Cloud & LiDAR Workflows in PyQGIS. It covers picking a resolution from the data rather than from habit, the two gridding algorithms and when each is right, filling holes without inventing landscape, and mosaicking a tiled survey into one seamless raster.

Resolution is decided by point density, not by tasteAt two metres every cell contains many ground returns and the surface is smooth but coarse. At one metre nearly every cell contains at least one return and the detail is real. At a quarter of a metre most cells contain no return at all, so the raster is mostly holes and whatever fills them is invention.4 ground returns per square metre, gridded three ways2.0 m cells~16 returns per cellsmooth, but detail is gone1.0 m cells~4 returns per cellthe honest choice here0.25 m cellsmost cells emptythe fill becomes the maprule of thumb: cell size ≈ 1 / √(ground returns per m²), rounded upand ground return density under canopy is a fraction of the survey headline figure

Prerequisites

  • QGIS 3.34 LTR or newer with the PDAL provider.
  • A cloud with ground returns classified, or a decision to build a surface model instead — see filtering and classifying a point cloud.
  • Enough disk for the output: a 1 m float32 raster over 1 km² is about 4 MB before compression.

Measure the density before choosing a resolution

The survey specification quotes an average pulse density over the whole flight. What you need is the density of the returns you are actually gridding, over the area you are actually gridding, and those two numbers can differ by a factor of ten under woodland.

from qgis.core import QgsPointCloudLayer

cloud = QgsPointCloudLayer("/data/clean/tile_0345_ground.copc.laz", "ground", "copc")
extent = cloud.extent()
area = extent.width() * extent.height()

density = cloud.pointCount() / area
print(f"{density:.2f} ground returns per m² over {area / 1e6:.2f} km²")
print("suggested cell size:", round((1 / density) ** 0.5, 2), "m")

Breakdown: Dividing by the bounding-box area is a slight underestimate of density wherever the tile is not fully covered, which errs in the safe direction — it suggests a coarser cell than strictly necessary. The square root converts a per-area density into a spacing: at 4 points/m² the mean spacing is 0.5 m, so a 0.5 m cell holds about one point on average and a 1 m cell holds about four. One point per cell is not enough; aim for a cell size where the typical count is a handful.

Grid it: binning or triangulation

There are two export algorithms and the difference between them is what happens in a cell with no points.

pdal:exportraster bins: each cell takes a statistic of the points that fall inside it, and a cell with no points is nodata. pdal:exportrastertin triangulates the points into a TIN and samples the resulting surface, so every cell inside the data's hull gets a value — interpolated across the facet it sits on.

import processing

dtm = processing.run("pdal:exportraster", {
    "INPUT": "/data/clean/tile_0345_ground.copc.laz",
    "RESOLUTION": 1.0,
    "TILE_SIZE": 1000,
    "OUTPUT": "/data/output/dtm_0345.tif",
})["OUTPUT"]

Breakdown: With the cloud already filtered to ground, no FILTER_EXPRESSION is needed — one less thing to get wrong. TILE_SIZE is the algorithm's internal chunk in map units, not the output tiling; raising it uses more memory and slightly fewer edge computations. The output is a single-band float32 GeoTIFF in the cloud's CRS, which means it inherits any CRS problem the cloud had, so fix that before this step rather than after.

Use the TIN variant when the binned output is holey:

dtm_tin = processing.run("pdal:exportrastertin", {
    "INPUT": "/data/clean/tile_0345_ground.copc.laz",
    "RESOLUTION": 1.0,
    "OUTPUT": "/data/output/dtm_0345_tin.tif",
})["OUTPUT"]

Breakdown: Triangulation is considerably slower and it is honest about what it is doing — it interpolates linearly between measured points, which is the most defensible thing to do in a gap. Where it goes wrong is across genuinely unmeasured areas such as a lake or a building footprint removed by the ground filter: it will span them with a flat triangle, producing a smooth ramp where there is no data at all. Mask those areas afterwards rather than trusting the fill.

Deciding what to do with the holes

A gap is data tooLeaving a gap as nodata is truthful and awkward downstream. Interpolating across it with a triangulation produces a continuous raster and quietly invents a surface under buildings and water. Filling only gaps below a size threshold and masking the rest keeps the raster usable while refusing to invent large areas.Three answers, one of which is usually rightleave it as nodatatruthful; breaks slopeand hillshade at the edgesinterpolate everythingcontinuous and convenientinvents ground under lakesfill small, mask largea threshold you can defendand a mask you can publishwhichever you choose, record it — a DEM with no gap policy stated is not reproduciblethe layer metadata is the natural place for it

The middle route — fill gaps up to a size you are willing to defend, and leave the rest — is a single GDAL algorithm.

filled = processing.run("gdal:fillnodata", {
    "INPUT": dtm,
    "BAND": 1,
    "DISTANCE": 6,
    "ITERATIONS": 0,
    "OUTPUT": "/data/output/dtm_0345_filled.tif",
})["OUTPUT"]

Breakdown: DISTANCE is in pixels, so on a 1 m raster this fills anything within six metres of a measured cell and leaves larger voids alone — which is roughly the right behaviour for gaps under isolated trees while refusing to fill a lake. ITERATIONS controls a smoothing pass over the filled values; zero keeps the interpolation crisp. Record the distance you used in the layer's metadata, as covered in reading and writing layer metadata, because six months later nobody will remember.

Mosaicking a tiled survey

Grid each tile, then merge. Gridding the whole survey in one call works and is simple, but it gives you no restart point and no parallelism.

import glob
import os

rasters = []
for tile in sorted(glob.glob("/data/clean/*_ground.copc.laz")):
    out = tile.replace("/clean/", "/output/").replace("_ground.copc.laz", ".tif")
    if not os.path.exists(out):
        processing.run("pdal:exportraster", {
            "INPUT": tile, "RESOLUTION": 1.0, "OUTPUT": out,
        })
    rasters.append(out)

processing.run("gdal:buildvirtualraster", {
    "INPUT": rasters,
    "RESOLUTION": 0,
    "OUTPUT": "/data/output/dtm_survey.vrt",
})

Breakdown: A VRT is the right mosaic here: it is a small XML file referencing the tiles, so building it is instant, it costs no extra disk, and every downstream algorithm treats it as one seamless raster. RESOLUTION: 0 means "average of the inputs", which is correct only when every input has the same cell size — pass 1 (highest) if they do not, and expect resampling. If the tiles were gridded with matching cell sizes and aligned origins the mosaic has no seams; if the origins do not align, a half-cell offset produces a visible grid of hairlines under hillshade.

Check the result before you use it

from qgis.core import QgsRasterLayer

dem = QgsRasterLayer("/data/output/dtm_survey.vrt", "dtm")
stats = dem.dataProvider().bandStatistics(1)
print(f"min {stats.minimumValue:.1f} max {stats.maximumValue:.1f} mean {stats.mean:.1f}")

Breakdown: Three numbers catch most disasters. A minimum near zero over hilly ground means nodata is being read as a value rather than as nodata. A maximum in the thousands over lowland means noise points survived the filter. A mean far from the survey's stated ground elevation means the vertical datum is not what you assumed. All three are five seconds of checking against hours of rework, and the same statistics call underpins calculating raster statistics.

QGIS version compatibility

pdal:exportraster and pdal:exportrastertin both arrived in 3.32 with the parameter names used here and are unchanged through 3.44. gdal:fillnodata and gdal:buildvirtualraster are long-standing GDAL wrappers, though fillnodata gained a MASK_LAYER parameter in 3.28 that is worth using where you already have a water mask.

Troubleshooting

  • The raster is all nodata. No points matched — either the cloud is empty over this extent or a filter expression matched nothing.
  • Terraced, stepped hillsides in the hillshade. The DEM was written as an integer type, or the source Z was rounded. Keep float32 throughout.
  • Hairline seams across the mosaic. Tile origins are not aligned to a common grid. Snap the extents to a multiple of the cell size when gridding.
  • A flat plateau where a lake is. TIN interpolation spanned an unmeasured area. Mask water before gridding, or grid with binning and fill only small gaps.
  • The output is enormous. No compression. Add -co COMPRESS=DEFLATE -co PREDICTOR=3 through the algorithm's creation options; float DEMs compress well with predictor 3.
  • Slope looks wrong by orders of magnitude. The DEM is in a geographic CRS, so horizontal units are degrees while heights are metres — the same trap described in terrain and interpolation analysis.

Conclusion

Compute the density, derive the cell size from it, bin when the returns are dense and triangulate when they are not, fill only the gaps you can defend and record the threshold, then mosaic through a VRT and check three statistics before anything else touches the raster. None of those steps is difficult; skipping them is what produces a DEM that is wrong in a way nobody notices until it has been published.

Frequently Asked Questions

Should I build a DEM or a DSM? A terrain model, built from ground returns, for anything about water, earthworks or accessibility. A surface model, built from the highest return per cell, for anything about what a signal, a view or the sun can reach. They are different products and neither substitutes for the other.

Can I get a canopy height model directly? Not in one call. Grid the surface and the terrain separately at the same resolution and origin, then subtract with the raster calculator. Matching the grids exactly is what makes the subtraction meaningful.

Which statistic does pdal:exportraster use per cell? It writes several bands when asked, including mean, minimum, maximum, count and inverse-distance-weighted values; the default output is the one most builds label as mean. Check parameterDefinitions() on your build rather than assuming, since the available outputs grew across releases.

Is a VRT mosaic slower than a real merged raster? Marginally, on random access across many tiles. For sequential work it is indistinguishable, and the flexibility of not duplicating hundreds of gigabytes usually wins. Materialise it with gdal:translate only when profiling says it matters.