My LiDAR DTM Has Holes, Spikes or Terraces

Problem statement

Three defects appear in almost every lidar terrain model, and they have different causes and different fixes.

Holes — cells with no value. Measured on a real survey at 1 m cells, 3.67% of the cells it covered had no ground point, and none had no return at all. (A bounding-box count says 10.93% and 7.47%, because a grid over a rotated survey has empty corners.)

Spikes — single cells far below or above their neighbours. The per-cell minimum is maximally sensitive to one low point.

Terraces — flat steps across a slope, from a cell size too coarse for the terrain, or from a fill that propagated one value across a gap.

Quick answer

Diagnose before fixing:

import numpy as np


def diagnose_dtm(dtm, counts, cell=1.0, surveyed=None):
    s = np.ones(dtm.shape, bool) if surveyed is None else surveyed   # cells that count
    holes = ~np.isfinite(dtm)
    print(f"  holes           {holes[s].mean():7.2%}")
    print(f"  single-point cells {(counts == 1)[s].mean():7.2%} "
          "(one bad point sets the whole cell)")

    padded = np.pad(dtm, 1, constant_values=np.nan)
    neighbours = np.stack([padded[a:a + dtm.shape[0], b:b + dtm.shape[1]]
                           for a in range(3) for b in range(3)
                           if (a, b) != (1, 1)])
    residual = dtm - np.nanmedian(neighbours, axis=0)

    print(f"  residual vs 8-neighbour median: "
          f"p1 {np.nanpercentile(residual, 1):+.2f} m, "
          f"p99 {np.nanpercentile(residual, 99):+.2f} m")
    print(f"  cells more than 2 m below their neighbours: "
          f"{int(np.nansum(residual < -2)):,}")
    return residual


residual = diagnose_dtm(dtm, counts, surveyed=inside)

  holes             3.67%
  single-point cells   6.94% (one bad point sets the whole cell)
  residual vs 8-neighbour median: p1 -0.29 m, p99 +0.32 m
  cells more than 2 m below their neighbours: 1

A cell built from one ground point, sitting two metres below its neighbours, is a spike. A cell sitting two metres below with forty supporting points is terrain.

Which holes, cell size first, spikes and terraces

1. Holes: work out which kind you have

Measured at 1 m (Example 1), with a third kind of empty cell that a bounding-box grid adds:

outside the survey footprint :  7.67% of the grid
surveyed, no return at all   :  0.00% of surveyed cells
returns but no ground point  :  3.67%
ground observed              : 96.33%

The kinds need different treatment. Cells with returns but no ground point may be recoverable with a better ground filter, or by accepting last returns. Cells with no returns at all cannot be improved by any processing — still water often returns nothing to the laser — and this survey had none. Cells outside the footprint were never surveyed and should not be in the count.

Interpolating a lake into a hillside is the classic result of not making this distinction.

2. Holes: choose the cell size before choosing a fill

Counted over surveyed cells (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%

Chasing a 0.5 m DTM on this survey means filling a third of it — at which point much of the product is a model, not a measurement.

3. Spikes: filter noise classes, then use a low percentile

The per-cell minimum takes the single lowest point. One multipath return below the surface sets that cell.

usable = ~np.isin(classification, [7, 18])       # declared noise
print(f"{(~usable).sum()} declared noise points")

2 declared noise points

This survey declares only two noise points, so the declared classes are not the problem. The fix that helps is replacing the strict minimum with a low percentile:

# 10th percentile of ground points in the cell, where there are enough
elevation = np.percentile(cell_points, 10) if len(cell_points) >= 5 \
    else cell_points.min()

With a median of four ground points per cell at 1 m, a percentile is only meaningful in the better-populated cells — which is exactly where the strict minimum is riskiest.

4. Spikes: check against the neighbourhood, not against a global threshold

A ground point three metres below its neighbours is a spike on a plateau and normal in a ravine. Compare each cell with the median of its eight neighbours and flag large negative residuals.

5. Terraces: check the cell size and the fill method

Terraces come from two places. A cell size much larger than the terrain's roughness quantises a smooth slope into steps. And a nearest-neighbour fill propagates one value across a gap, producing a plateau with a hard edge.

Nearest-neighbour fill is honest — the facets are visible and read as artefacts. A smooth interpolation over the same gap looks like terrain and is not. Choose deliberately, and ship the mask either way.

Separating hole types and filtering spikes

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. The setup also builds the 1 m lowest-ground-point DTM with its per-cell ground counts, and inside marks the cells within the survey's footprint: the patch is turned 2.4° against the UTM grid, so the grid's bounding box has corners that were never surveyed.

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)

from rasterio.transform import from_origin

# the 1 m lowest-ground-point DTM and its ground-point counts, over the bounding box
CELL = 1.0
left, top = x.min(), y.max()
width, height = int(np.ceil((x.max() - left) / CELL)), int(np.ceil((top - y.min()) / CELL))
transform = from_origin(left, top, CELL, CELL)
ground = classification == 2
flat = (np.clip(((top - y[ground]) / CELL).astype(int), 0, height - 1) * width
        + np.clip(((x[ground] - left) / CELL).astype(int), 0, width - 1))
dtm = np.full(width * height, np.inf)
np.minimum.at(dtm, flat, z[ground])
dtm[np.isinf(dtm)] = np.nan
dtm = dtm.reshape(height, width)
counts = np.bincount(flat, minlength=width * height).reshape(height, width)
inside = surveyed(transform, dtm.shape)

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

Example 1 — separating the hole types

import numpy as np


def classify_holes(x, y, classification, surveyed_cells, cell=1.0, ground_class=2):
    """Which empty cells were never surveyed, which had no return, and which had no ground?"""
    left, top = x.min(), y.max()
    width = int(np.ceil((x.max() - left) / cell))
    height = int(np.ceil((top - y.min()) / cell))

    def occupancy(mask):
        col = np.clip(((x[mask] - left) / cell).astype(int), 0, width - 1)
        row = np.clip(((top - y[mask]) / cell).astype(int), 0, height - 1)
        grid = np.zeros(width * height, bool)
        grid[row * width + col] = True
        return grid.reshape(height, width)

    any_return = occupancy(np.ones(len(x), bool))
    ground = occupancy(classification == ground_class)

    outside = ~surveyed_cells
    no_data = surveyed_cells & ~any_return
    no_ground = any_return & ~ground

    print(f"  outside the footprint   {outside.mean():7.2%}  of the grid (never surveyed)")
    s = surveyed_cells
    print(f"  no return at all        {no_data[s].mean():7.2%}  "
          "(water or a gap — not recoverable)")
    print(f"  returns but no ground   {no_ground[s].mean():7.2%}  "
          "(canopy or buildings — may be recoverable)")
    print(f"  ground observed         {ground[s].mean():7.2%}")
    return {"outside": outside, "no_data": no_data, "no_ground": no_ground, "ground": ground}


holes = classify_holes(x, y, classification, inside)

  outside the footprint     7.67%  of the grid (never surveyed)
  no return at all          0.00%  (water or a gap — not recoverable)
  returns but no ground     3.67%  (canopy or buildings — may be recoverable)
  ground observed          96.33%

Filling only no_ground and leaving no_data and outside empty is the honest default. It also produces a DTM whose holes are lakes, which looks right on a hillshade.

Example 2 — a spike filter that uses the neighbourhood

import numpy as np
from scipy.ndimage import generic_filter


def remove_spikes(dtm, counts, drop_below=1.5, min_support=3):
    """Flag cells far below their neighbours, especially thinly supported ones."""
    def neighbour_median(window):
        centre = window[len(window) // 2]
        others = np.delete(window, len(window) // 2)
        others = others[np.isfinite(others)]
        return np.nan if others.size < 3 else np.nanmedian(others)

    medians = generic_filter(dtm, neighbour_median, size=3, mode="constant",
                             cval=np.nan)
    residual = dtm - medians

    spike = (residual < -drop_below) & np.isfinite(residual)
    thin = spike & (counts < min_support)

    print(f"  {int(spike.sum()):,} cells more than {drop_below} m below "
          f"their neighbours ({spike.mean():.3%})")
    print(f"  of those, {int(thin.sum()):,} have fewer than {min_support} "
          "supporting ground points")

    cleaned = dtm.copy()
    cleaned[thin] = medians[thin]
    print(f"  replaced {int(thin.sum()):,} thinly supported spikes")
    return cleaned, spike


cleaned, spike = remove_spikes(dtm, counts)

  2 cells more than 1.5 m below their neighbours (0.000%)
  of those, 2 have fewer than 3 supporting ground points
  replaced 2 thinly supported spikes

Replacing only the thinly supported spikes is the important restriction. A cell three metres below its neighbours with forty ground points in it is a ditch, and smoothing it away destroys real terrain. This survey was clean: two cells fell more than 1.5 m below their neighbours, and both rested on fewer than three ground points.

Example 3 — a fill with a distance limit and a mask

import numpy as np
from scipy.ndimage import distance_transform_edt


def fill_dtm(dtm, cell=1.0, max_distance_m=5.0, method="nearest",
             never_fill=None):
    """Fill small gaps, leave large ones, and record which is which."""
    observed = np.isfinite(dtm)
    if never_fill is not None:
        fillable = ~observed & ~never_fill
    else:
        fillable = ~observed

    distance = distance_transform_edt(~observed) * cell
    indices = distance_transform_edt(~observed, return_distances=False,
                                     return_indices=True)

    filled = dtm.copy()
    target = fillable & (distance <= max_distance_m)
    filled[target] = dtm[tuple(indices)][target]

    print(f"  {int(observed.sum()):,} observed, "
          f"{int(target.sum()):,} filled (within {max_distance_m} m), "
          f"{int((~observed & ~target).sum()):,} left empty")
    if never_fill is not None:
        print(f"  {int(never_fill.sum()):,} cells excluded from filling "
              "(never surveyed, or no returns at all)")
    return filled, observed


filled, observed = fill_dtm(cleaned, never_fill=holes["outside"] | holes["no_data"])
print(f"  surveyed cells still empty: {np.isnan(filled[inside]).mean():.2%}")

  684,275 observed, 26,019 filled (within 5.0 m), 57,958 left empty
  58,898 cells excluded from filling (never surveyed, or no returns at all)
  surveyed cells still empty: 0.00%

Passing never_fill is what keeps the lakes as lakes. Without it, the nearest-neighbour fill extends the shoreline elevation across the water and produces a flat plateau where the lake was. Here it held the grid's unsurveyed corners — the only cells with no returns — and every hole inside the survey was within 5 m of a ground point and got filled.

Why the minimum reduction is fragile

Why the minimum reduction is fragile

The per-cell minimum is chosen because ground classification errs upward: low vegetation and building edges are occasionally accepted as ground, and genuine ground points are rarely invented. Taking the lowest is the best defence against that.

The cost is total sensitivity in the other direction. A single point below the true surface — multipath, a sensor artefact, a misregistered flight line — sets the cell entirely, regardless of how many good points are also there.

With a median of four ground points per cell at 1 m, most cells have very little redundancy. A 10th percentile is barely different from the minimum at four points, which is why the honest fix is often a coarser cell rather than a cleverer statistic.

Why some holes cannot be fixed

Water absorbs near-infrared strongly. A pulse hitting still water returns nothing detectable, so the cell is empty and no reprocessing changes that.

A floor that does not move with cell size is the signature of real absence rather than sampling — but check what it is before calling it water. Counted over the bounding box, this survey showed a floor of 7.47% at 1 m and 7.29% at 5 m; counted over the surveyed cells it vanished, because every one of those cells was outside the survey.

Cells with returns but no ground point are a different matter. There, the laser reached something; the ground classifier either could not see the ground through the canopy or declined to accept a candidate. Some of those are recoverable with a different filter.

Why nearest-neighbour fill is often the right choice

A smooth interpolation across a hole produces something that looks like terrain: continuous, plausibly sloped, and indistinguishable from measurement in the output raster.

Nearest-neighbour fill produces visible facets. That is uglier and more honest — a reader looking at the hillshade can see where the data stops.

Either way the mask is what matters. A second band recording observed versus filled costs one band and lets every downstream user decide.

Why terracing is usually a cell-size problem

If the cell is much larger than the horizontal scale of the terrain's variation, each cell reduces a range of elevations to one number, and the result is a staircase.

The tell is that the step height matches the elevation change across one cell. If a 5 m DTM on a 20% slope has 1 m steps, that is exactly the cell size times the gradient, and the fix is a finer cell — if the point density supports one.

Where it does not, the honest response is to accept the coarser product rather than to interpolate a smooth surface that implies detail the survey never measured.

Water, fill distance and no-ground cells

• Separate "no returns" from "no ground points". Only the second is recoverable.
• Never fill across water. Pass an exclusion mask.
• Limit the fill distance. Filling a 2 m gap is reasonable; a 200 m one is fabrication.
• Ship an observation mask as a second band.
• Filter noise classes before any minimum-based reduction.
• Only smooth thinly supported spikes. A well-supported low cell is terrain.
• Terraces usually mean the cell is too coarse for the slope.
• Never fill with zero. It reads as a hole in every hillshade and every flow model.

Internal links

• How to create a DTM from LiDAR ground points in Python — where these defects originate
• Point density explained — the occupancy sweep that sets the cell size
• DSM, DTM and CHM from LiDAR: how each surface is derived — how holes propagate into a CHM
• How to filter a point cloud by class and return number — removing noise classes
• How to create a canopy height model from LiDAR — where negative values appear
• How to fill sinks and derive flow direction from a DEM — a different kind of hole
• Slope values are wrong or absurdly steep — what a spiky DTM does downstream
• Contours come out jagged or full of tiny rings — the same noise seen through contours

FAQ

Why does my lidar DTM have holes?

Two reasons with different fixes: cells where no pulse returned at all (water, flight-line edges — none here) and cells with returns but no ground-classified point (canopy, buildings — 3.67% of the surveyed 1 m cells here). Count only cells inside the survey, or its bounding box's corners join the first group.

Should I fill the holes?

Fill the small ones, limit the fill distance, exclude water, and ship a mask saying which cells were observed.

Why does my DTM have single-cell pits?

The per-cell minimum takes the lowest point, so one low noise return sets the whole cell. Filter noise classes, and replace only thinly supported outliers with the neighbourhood median.

Why is my DTM terraced?

Usually the cell size is too coarse for the slope, so each cell reduces a range of elevations to one value. Check whether the step height equals the cell size times the gradient.

Can I just smooth the DTM?

Smoothing treats holes, spikes and terraces identically and fixes none properly. It also removes real terrain detail — a genuine ditch looks exactly like a spike to a smoothing filter.

What cell size avoids holes?

On this survey, 2 m left 0.13% of the surveyed cells empty and 1 m left 3.67%. Below 1 m the hole fraction rises steeply: a third at 0.5 m.

Why should I never fill with zero?

Zero is a valid elevation. A filled zero in a valley reads as a deep pit in every hillshade, contour set and flow-direction model.