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, 10.93% of cells had no ground point, against a floor of 7.47% with no returns at all.

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):
    holes = ~np.isfinite(dtm)
    print(f"  holes           {holes.mean():7.2%}")
    print(f"  single-point cells {(counts == 1).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

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.

Step-by-step solution #

1. Holes: work out which kind you have #

cells with no return at all :  7.47%   water, or outside the flight line
returns but no ground point :  3.46%   canopy or buildings blocking the ground
ground observed             : 89.07%

The two 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 — water absorbs the near-infrared laser and returns nothing.

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 #

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 is about 7.3%. Chasing a 0.5 m DTM on this survey means filling 38% of it — at which point most 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

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.

Code examples #

Example 1 — separating the hole types #

import numpy as np


def classify_holes(x, y, classification, cell=1.0, ground_class=2):
    """Which empty cells are water, and which are unclassified 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)

    no_data = ~any_return
    no_ground = any_return & ~ground

    print(f"  no return at all        {no_data.mean():7.2%}  "
          "(water or outside the flight line — not recoverable)")
    print(f"  returns but no ground   {no_ground.mean():7.2%}  "
          "(canopy or buildings — may be recoverable)")
    print(f"  ground observed         {ground.mean():7.2%}")
    return {"no_data": no_data, "no_ground": no_ground, "ground": ground}

  no return at all           7.47%  (water or outside the flight line — not recoverable)
  returns but no ground      3.46%  (canopy or buildings — may be recoverable)
  ground observed           89.07%

Filling only no_ground and leaving no_data 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

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.

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 "
              "(no returns at all)")
    return filled, observed

Passing never_fill=no_data 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.

Explanation #

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.

The measured floor of 7.47% of cells with no return at all does not move with cell size — it was 7.29% at 5 m cells, essentially the same. That flatness is the signature of real absence rather than sampling.

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.

Edge cases or notes #

• 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 — 7.47% here) and cells with returns but no ground-classified point (canopy, buildings — a further 3.46%).

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 7.51% empty against an irreducible floor of about 7.3%. Below 1 m the hole fraction rises steeply.

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.