⚙️ How-To LiDAR & point clouds Terrain & elevation LiDAR & Point Clouds · step 10 of 14 Updated

How to Create a DTM from LiDAR Ground Points in Python

Problem statement

A digital terrain model from lidar is a per-cell reduction of the ground-classified returns. The code is short. The decisions are not:

  • What cell size? Fine enough to be useful, coarse enough to be occupied.
  • Which reduction? Minimum, mean or something more careful.
  • What about the empty cells? At any usable resolution there will be some.
  • Is the CRS metric? If not, every distance in the process is wrong.

Measured on a real 3DEP survey — 12,968,770 points, 16.9 per square metre, 22.9% classified ground — at 1 m cells:

DTM  15.9 .. 194.3 m,  10.93% of cells with no ground point

One cell in nine has to be filled or left empty, and that decision changes what the product is.

Quick answer 

import numpy as np
import rasterio
from rasterio.transform import from_origin


def dtm_from_points(x, y, z, classification, cell=1.0, crs="EPSG:32605"):
    """Per-cell minimum of ground-classified returns."""
    ground = classification == 2
    if not ground.any():
        raise ValueError("no ground-classified points — run a ground filter first")

    left, top = x.min(), y.max()
    width = int(np.ceil((x.max() - left) / cell))
    height = int(np.ceil((top - y.min()) / cell))

    col = np.clip(((x[ground] - left) / cell).astype(int), 0, width - 1)
    row = np.clip(((top - y[ground]) / cell).astype(int), 0, height - 1)

    dtm = np.full(width * height, np.inf)
    np.minimum.at(dtm, row * width + col, z[ground])
    dtm[np.isinf(dtm)] = np.nan

    counts = np.zeros(width * height)
    np.add.at(counts, row * width + col, 1)

    print(f"  {width} x {height} at {cell} m, "
          f"{np.isnan(dtm).mean():.2%} empty, "
          f"median {np.median(counts[counts > 0]):.0f} ground points per cell")
    return (dtm.reshape(height, width), counts.reshape(height, width),
            from_origin(left, top, cell, cell))

np.minimum.at does the whole reduction in one vectorised pass. A Python loop over 3 million ground points takes minutes; this takes under a second.

Ground-classified returns falling into grid cells, with each cell taking the minimum elevation and empty cells left as no data.
One reduction per cell, and a hole wherever no ground point landed.

Step-by-step solution

1. Reproject to a metric CRS first

The survey used here is stored in EPSG:3857, where its extent reads as 1,613 m and is really 877 m. A DTM built in those units has cells 1.92 times larger than requested, and any slope derived from it is wrong by the same factor.

from pyproj import Transformer
transformer = Transformer.from_crs(3857, 32605, always_xy=True)
x, y = transformer.transform(las.x, las.y)

2. Filter noise before reducing

The per-cell minimum is maximally sensitive to a single low point. Drop the noise classes explicitly:

ground = (classification == 2) & ~np.isin(classification, [7, 18])

This survey contained exactly two points classified as low noise — a reminder that classified noise is the easy case. Unclassified low points are the ones that hurt.

3. Choose the cell size from occupancy 

0.25 m:  78.32% of cells have no ground point
0.50 m:  38.11%
1.00 m:  10.93%
2.00 m:   7.51%
5.00 m:   7.33%

The floor at about 7.3% is water, which returns nothing — no cell size removes it. Take the finest size whose empty fraction is close to that floor: 1 m here, or 2 m for a gap-free product.

4. Decide what to do with the empty cells

Three defensible options, and one bad one:

  • Leave them as NoData. Honest, and awkward for anything that needs a continuous surface.
  • Interpolate, and record it. Fill from surrounding ground cells and ship a mask saying which cells are observed.
  • Coarsen until they mostly vanish. Loses resolution everywhere to fix a problem in one cell in nine.
  • Fill with zero. Never. A zero elevation in a valley reads as a hole in every hillshade and every flow model.

5. Consider using last returns to reduce holes — carefully 

DTM from last returns vs from classified ground
  holes    7.48% vs 10.93%
  bias    -0.040 m mean, p1 -0.33 m, p99 +0.00 m

Last returns fill a third of the holes, and the resulting surface is never higher than the classified-ground one — because taking a minimum over a superset can only go down. The 4 cm mean bias is downward, driven by low returns the classifier had excluded.

That is the trade: fewer holes, a slightly lower and noisier surface. For hydrological work the extra continuity often wins; for volume calculations the bias may not.

Empty DTM cells left as no data, interpolated with a mask, or filled from last returns, with the trade-offs of each.
Every fill is a model. The mask is what keeps the product honest.

Code examples

Example 1 — DTM with interpolation and an observation mask 

import numpy as np
import rasterio
from rasterio.transform import from_origin
from scipy.ndimage import distance_transform_edt


def dtm_with_fill(x, y, z, classification, cell=1.0, crs="EPSG:32605",
                  max_fill_distance=None, out_path=None):
    """DTM plus a mask band recording which cells were observed."""
    dtm, counts, transform = dtm_from_points(x, y, z, classification, cell, crs)
    observed = np.isfinite(dtm)

    # nearest-observed fill, limited by distance so large voids stay empty
    indices = distance_transform_edt(~observed, return_distances=False,
                                     return_indices=True)
    distance = distance_transform_edt(~observed)
    filled = dtm[tuple(indices)]

    if max_fill_distance is not None:
        too_far = distance > (max_fill_distance / cell)
        filled[too_far] = np.nan
        print(f"  left {too_far.sum():,} cells empty "
              f"(further than {max_fill_distance} m from any ground point)")

    print(f"  filled {int((~observed & np.isfinite(filled)).sum()):,} cells "
          f"({(~observed & np.isfinite(filled)).mean():.2%})")

    if out_path:
        profile = dict(driver="GTiff", height=dtm.shape[0], width=dtm.shape[1],
                       count=2, dtype="float32", crs=crs, transform=transform,
                       nodata=np.nan, compress="deflate", tiled=True)
        with rasterio.open(out_path, "w", **profile) as dst:
            dst.write(filled.astype("float32"), 1)
            dst.write(observed.astype("float32"), 2)
            dst.set_band_description(1, "elevation (m)")
            dst.set_band_description(2, "1 = observed, 0 = interpolated")
            dst.update_tags(cell_size=cell, reduction="min",
                            source_class="2 (ground)")
    return filled, observed, transform

Nearest-neighbour fill is deliberately crude: it produces visible facets, which is honest. A smooth interpolation over a large void looks like terrain and is not.

max_fill_distance is the important parameter. Filling a 2 m gap between ground points is reasonable; filling a 200 m lake is fabrication.

Example 2 — comparing reductions on your own data 

import numpy as np


def compare_reductions(x, y, z, classification, cell=1.0):
    """Minimum, mean and a low-percentile reduction, side by side."""
    ground = classification == 2
    left, top = x.min(), y.max()
    width = int(np.ceil((x.max() - left) / cell))
    height = int(np.ceil((top - y.min()) / cell))
    col = np.clip(((x[ground] - left) / cell).astype(int), 0, width - 1)
    row = np.clip(((top - y[ground]) / cell).astype(int), 0, height - 1)
    flat = row * width + col
    zg = z[ground]

    minimum = np.full(width * height, np.inf)
    np.minimum.at(minimum, flat, zg); minimum[np.isinf(minimum)] = np.nan

    total = np.zeros(width * height); counts = np.zeros(width * height)
    np.add.at(total, flat, zg); np.add.at(counts, flat, 1)
    mean = np.where(counts > 0, total / np.maximum(counts, 1), np.nan)

    difference = np.abs(minimum - mean)
    print(f"  min vs mean: mean |diff| {np.nanmean(difference):.3f} m, "
          f"p99 {np.nanpercentile(difference, 99):.3f} m, "
          f"max {np.nanmax(difference):.2f} m")
    return minimum.reshape(height, width), mean.reshape(height, width)
  min vs mean: mean |diff| 0.078 m, p99 0.326 m, max 0.94 m

An 8 cm typical difference means the choice barely matters on this terrain. Run it on your own data before assuming that transfers — on steep slopes with large cells the gap grows quickly.

Example 3 — a fallback ground filter when the survey has none 

import numpy as np


def simple_ground_filter(x, y, z, coarse_cell=10.0, tolerance=0.5,
                         fine_cell=1.0):
    """Two-pass minimum-surface filter, for clouds with no ground class."""
    def cell_min(px, py, pz, cell):
        left, top = x.min(), y.max()
        width = int(np.ceil((x.max() - left) / cell))
        height = int(np.ceil((top - y.min()) / cell))
        col = np.clip(((px - left) / cell).astype(int), 0, width - 1)
        row = np.clip(((top - py) / cell).astype(int), 0, height - 1)
        grid = np.full(width * height, np.inf)
        np.minimum.at(grid, row * width + col, pz)
        grid[np.isinf(grid)] = np.nan
        return grid.reshape(height, width), row, col, (left, top, width, height)

    coarse, row_c, col_c, meta = cell_min(x, y, z, coarse_cell)
    seed = coarse[row_c, col_c]
    candidate = z <= seed + tolerance

    fine, row_f, col_f, _ = cell_min(x[candidate], y[candidate],
                                     z[candidate], fine_cell)
    print(f"  {candidate.mean():.1%} of points kept as ground candidates")
    print(f"  {np.isnan(fine).mean():.2%} of fine cells empty")
    return fine, candidate

This is a crude approximation of what a real ground filter does — progressive TIN densification or cloth simulation — and it will clip ridges and accept flat roofs. Use it to get a first look, and use a proper filter for anything that matters.

Explanation

Why the minimum is the standard reduction

Ground classification is not perfect, and its errors are asymmetric: low vegetation and building edges are occasionally accepted as ground, while genuine ground points are rarely invented.

So the ground-classified set within a cell tends to contain the true ground plus a few things above it. The minimum is the most likely true ground, and the mean is biased upward by the contaminants.

The cost is sensitivity in the other direction. One multipath return below the surface sets the whole cell, so filtering noise classes first is not optional.

Why filling is a model and should be labelled

A filled cell is not a measurement. Nothing in a single-band GeoTIFF distinguishes an observed 42.3 m from an interpolated one, and every downstream user will treat them identically.

Shipping an observation mask as a second band costs one band and makes the product honest. It also makes it possible for a careful user to mask the derived slope, hillshade or watershed where the terrain was invented.

Why last returns are a tempting shortcut

Every ground point is a last return — measured at 100.0% here. So filtering to last returns keeps all the ground and adds whatever else reached the ground level, and taking a per-cell minimum then approximates the ground surface.

The measurement shows both sides. Holes fall from 10.93% to 7.48%, and the surface is 4 cm lower on average with a 1st percentile of −0.33 m. It cannot be higher, because a minimum over a superset only goes down.

The right choice depends on which error you can tolerate. Continuity matters more for flow routing; unbiased elevation matters more for volumes.

Why the CRS check comes first

Everything in this workflow is metric: cell size, fill distance, the eventual slope. In Web Mercator at 58.6° north, all of those are inflated by 1.92, and none of it raises an error.

The tell is that the extent looks too big for the area you thought you had. Check it once, at load time, and reproject before anything else.

A DTM from last returns having 7.48 percent holes against 10.93 percent from classified ground.
Last returns fill a third of the holes and bias the surface 4 cm downward. Pick the error you can tolerate.

Edge cases or notes

  • Reproject to a projected CRS first. Web Mercator inflated distances by 1.92 here.
  • Drop noise classes before taking a minimum.
  • Choose the cell size from occupancy, not from nominal spacing.
  • Never fill NoData with zero.
  • Limit the fill distance so large voids stay empty.
  • Ship an observation mask as a second band.
  • np.minimum.at is the fast reduction; a Python loop is minutes rather than seconds.
  • Water returns nothing. Do not interpolate a lake into a hillside.

FAQ

How do I make a DTM from lidar in Python?

Reproject to a metric CRS, select ground-classified returns, and take the per-cell minimum with np.minimum.at. Then decide explicitly what to do with the empty cells.

What cell size should I use?

The finest whose empty-cell fraction is close to the irreducible floor. On a 16.9 points per square metre survey that was 1 m, leaving 10.93% empty against a 7.3% floor.

Should I use the minimum or the mean per cell?

Minimum. It resists the upward bias from misclassified low vegetation, at the cost of sensitivity to low noise — so filter noise classes first. They differed by 8 cm on average here.

How do I fill the holes?

Nearest-neighbour or IDW from surrounding ground cells, with a maximum fill distance so large voids stay empty, and an observation mask shipped alongside.

Can I use last returns instead of ground classification?

It reduces holes from 10.93% to 7.48% and lowers the surface by 4 cm on average, because a minimum over a superset can only fall. Useful for continuity, biased for volumes.

Why is my DTM full of spikes?

Low noise points setting the per-cell minimum. Filter noise classes, and consider a low percentile instead of the strict minimum.

What if my cloud has no ground classification?

Run a ground filter. A crude two-pass minimum-surface filter gets a first look; progressive TIN densification or cloth simulation is what production work uses.