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, 18.3 per square metre, 22.9% classified ground — at 1 m cells, 3.67% of the surveyed cells had no ground point (quick answer). One cell in 27 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))


dtm, counts, transform = dtm_from_points(x, y, z, classification)
inside = surveyed(transform, dtm.shape)
print(f"  inside the survey footprint: {np.isnan(dtm[inside]).mean():.2%} empty, "
      f"elevation {np.nanmin(dtm):.1f} .. {np.nanmax(dtm):.1f} m")

  877 x 876 at 1.0 m, 10.93% empty, median 4 ground points per cell
  inside the survey footprint: 3.67% empty, elevation 15.9 .. 194.3 m

The function's own "empty" counts the whole bounding box, corners included; the second line counts only the ground that was surveyed.

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.

Metric CRS, noise removed, cell size from occupancy

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 841 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)
print(f"{x.max() - x.min():.0f} x {y.max() - y.min():.0f} m in UTM, "
      f"{las.x.max() - las.x.min():.0f} x {las.y.max() - las.y.min():.0f} in Web Mercator")

877 x 876 m in UTM, 1613 x 1613 in Web Mercator

The UTM figure is a bounding box of a slightly rotated square; the patch itself is 841 m across.

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])
print(f"{ground.sum():,} ground points; noise classes present: "
      f"{ {c: int((classification == c).sum()) for c in (7, 18)} }")

2,965,430 ground points; noise classes present: {7: 2, 18: 0}

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

Counted over the surveyed cells only (point density explained runs the sweep):

0.25 m:  76.53% of cells have no ground point
0.50 m:  32.99%
1.00 m:   3.67%
2.00 m:   0.13%
5.00 m:   0.00%

Take the finest size whose empty fraction you can live with: 1 m here, or 2 m for a nearly gap-free product. Counted over the bounding box instead, every row gains the 7.7% of corners that were never surveyed, and the sweep appears to hit a floor that is not in the data.

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, surveyed 1 m cells
  holes    0.00% vs 3.67%
  bias    -0.040 m mean, p1 -0.33 m, p99 +0.00 m

(The filter guide measures this.) Last returns fill every hole here, 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.

A DTM with interpolation and an observation mask

The data used on this page

The USGS 3DEP lidar survey of Brooks Camp, Katmai National Park, Alaska (2012), from the USGS Entwine Point Tile archive on AWS, AK_BrooksCamp_2012: every node under octree node 4-3-13-7 — 287 LAZ files, 12,968,770 points — joined into tile.laz, reprojected from Web Mercator to UTM zone 5N, and loaded as arrays. footprint is the convex hull of the points: the survey's square patch is turned 2.4° against the UTM grid, so a grid over its bounding box has empty corners that were never surveyed, and surveyed() marks the cells that were.

import json
import os
import urllib.request

import laspy

EPT = "https://s3-us-west-2.amazonaws.com/usgs-lidar-public/AK_BrooksCamp_2012"
ROOT = (4, 3, 13, 7)                                     # depth, x, y, z of the octree node


def under(key):
    d, x, y, z = map(int, key.split("-"))
    step = 2 ** (d - ROOT[0])
    return d >= ROOT[0] and (x // step, y // step, z // step) == ROOT[1:]


hierarchy = json.load(urllib.request.urlopen(f"{EPT}/ept-hierarchy/0-0-0-0.json"))
nodes = sorted(k for k in hierarchy if under(k))
os.makedirs("nodes", exist_ok=True)
for key in nodes:
    path = f"nodes/{key}.laz"
    if not os.path.exists(path):
        urllib.request.urlretrieve(f"{EPT}/ept-data/{key}.laz", path + ".part")
        os.replace(path + ".part", path)

if not os.path.exists("tile.laz"):                       # the whole subtree as one file
    with laspy.open(f"nodes/{nodes[0]}.laz") as first:
        header = first.header
    with laspy.open("tile.laz", mode="w", header=header) as writer:
        for key in nodes:
            writer.write_points(laspy.read(f"nodes/{key}.laz").points)
with laspy.open("tile.laz") as reader:
    print(len(nodes), "nodes,", f"{reader.header.point_count:,} points,",
          f"{sum(os.path.getsize(f'nodes/{k}.laz') for k in nodes) / 1e6:.1f} MB of nodes")

import numpy as np
from pyproj import Transformer

las = laspy.read("tile.laz")
x, y = Transformer.from_crs(3857, 32605, always_xy=True).transform(las.x, las.y)   # UTM 5N metres
z = np.asarray(las.z)
classification = np.asarray(las.classification)
return_number = np.asarray(las.return_number)
number_of_returns = np.asarray(las.number_of_returns)

import shapely
from scipy.spatial import ConvexHull

hull = ConvexHull(np.c_[x, y])
footprint = shapely.Polygon(np.c_[x, y][hull.vertices])     # the ground the survey covers


def surveyed(transform, shape):
    """True for grid cells whose centre lies inside the survey footprint."""
    rows, cols = np.indices(shape)
    cx, cy = transform * (cols + 0.5, rows + 0.5)
    return shapely.contains_xy(footprint, cx, cy)

287 nodes, 12,968,770 points, 82.8 MB of nodes

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


filled, observed, transform = dtm_with_fill(x, y, z, classification, max_fill_distance=5.0,
                                            out_path="dtm.tif")
inside = surveyed(transform, filled.shape)
print(f"  inside the footprint: {observed[inside].mean():.2%} observed, "
      f"{np.isnan(filled[inside]).mean():.2%} still empty")

  877 x 876 at 1.0 m, 10.93% empty, median 4 ground points per cell
  left 42,556 cells empty (further than 5.0 m from any ground point)
  filled 41,421 cells (5.39%)
  inside the footprint: 96.33% observed, 0.00% still empty

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. Here every hole inside the survey was within 5 m of a ground point and got filled; the 42,556 cells left empty are the unsurveyed corners of the bounding box, which is exactly what the limit is for.

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)


minimum, mean = compare_reductions(x, y, z, classification)

  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


fine, candidate = simple_ground_filter(x, y, z)
ground = classification == 2
print(f"  of the classified ground, {candidate[ground].mean():.1%} kept; "
      f"of the candidates, {ground[candidate].mean():.1%} are classified ground")

  12.3% of points kept as ground candidates
  69.60% of fine cells empty
  of the classified ground, 28.0% kept; of the candidates, 52.1% are classified ground

This is a crude approximation of what a real ground filter does — progressive TIN densification or cloth simulation — and on this survey it shows. With a 10 m seed grid and a 0.5 m tolerance it kept only 28% of the points the vendor classified as ground, and only half of what it kept was ground: the terrain here rises 178 m across 841 m, so a 10 m cell's minimum sits well below most of the cell, and where three quarters of points come from multi-return pulses the lowest return in a cell is often not ground. Use it for a first look where the terrain is flat and open, and a proper filter for anything that matters.

Why the minimum is the standard reduction

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 inside the survey fall from 3.67% to none, 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.

Noise classes and fill limits

• 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.
• Count holes inside the survey footprint. A grid over a rotated patch's bounding box adds corners that were never surveyed — 7.7% here.

Internal links

• DSM, DTM and CHM from LiDAR: how each surface is derived — the three surfaces compared
• Point density explained — choosing the cell size
• My LiDAR DTM has holes, spikes or terraces — diagnosing a bad DTM
• How to filter a point cloud by class and return number — selecting ground points
• laspy returns the wrong coordinates: scale, offset and CRS — the reprojection step
• How to rasterise a point cloud to a grid in Python — the general reduction
• How to calculate slope and aspect from a DEM in Python — what the DTM feeds
• How to fill sinks and derive flow direction from a DEM — where holes matter most

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 you can accept, counted over the cells the survey actually covered. On an 18.3 points per square metre survey that was 1 m, leaving 3.67% empty; 2 m left 0.13%.

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 removed every hole here — 3.67% of surveyed cells — and lowered 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.