Fixing Holes and Voids in Drone-Derived DSMs

You rasterized a dense cloud into a DSM, opened it in QGIS, and found dark craters punched through the surface: a rectangular void over the retention pond, ragged holes across the glazed atrium roof, and a scatter of pinholes over freshly poured concrete. Those cells hold your NoData sentinel (-9999), and any downstream volume calculation, contour, or hillshade will either fail or produce garbage where they fall. This page explains why those voids form, gives one focused routine that closes them, and shows where the fix ends and a re-processing decision begins.

# Why NoData voids appear in a DSM

writers.gdal can only interpolate a cell from points that exist near it. A void appears wherever no point landed within the cell’s search radius, and drone photogrammetry produces point-free zones for two distinct reasons. The first is physical: dense stereo matching needs texture, so calm water, clean glass, wet asphalt, and uniform white membrane roofs return almost no matchable features and therefore no points. No amount of rasterization can invent a surface the sensor never reconstructed. The second is sampling: where coverage is merely thin — a flight-line gap, a single missed overlap strip — points exist but are sparser than the cell grid, so cells fall between them and go empty.

The two causes want different fixes. Sampling voids close by widening the writers.gdal radius so each cell reaches a real neighbour; texture voids have no nearby point at any radius and must be filled by interpolating across the hole from its rim. Conflating them is the usual mistake: cranking radius to swallow a genuine texture void smears the entire raster to erase one pond. The correct order is to set radius for the sampling floor, then fill the physical voids as a separate raster post-process. This void handling is the most common defect coming out of generating DSM and DTM from point clouds with PDAL, and it must be resolved before the raster reaches exporting cloud-optimized GeoTIFFs with Rasterio.

Figure 1 — The raster cannot tell you which void it is holding. Filling first and asking later produces a smooth surface that invents elevations over water and hides a mis-tuned ground filter.

# Minimal reproducible solution

The routine below reads the holed DSM, builds a boolean mask of NoData cells, and fills them with rasterio.fill.fillnodata, which interpolates each void inward from its valid rim using an inverse-distance weighting bounded by max_search_distance (in pixels). Bounding the search distance is the safety valve: it lets small pinholes and thin gaps close while leaving a genuinely huge void — a whole lake — untouched, so you never hallucinate a flat plane across an area the sensor never saw.

import numpy as np
import rasterio
from rasterio.fill import fillnodata

NODATA = -9999.0


def fill_dsm_voids(src_path: str, dst_path: str,
                   max_search_distance: float = 12.0,
                   smoothing_iterations: int = 0) -> float:
    """Fill small NoData voids in a DSM and return the fraction filled.

    max_search_distance is in PIXELS: a cell is only filled if a valid
    cell lies within that many pixels, so large physical voids stay NoData.
    """
    with rasterio.open(src_path) as src:
        band = src.read(1)                       # single elevation band
        profile = src.profile
        # Mask is True where data is VALID (fillnodata's convention).
        valid = band != NODATA
        before = (~valid).mean()                 # NoData fraction before fill

        filled = fillnodata(
            band,
            mask=valid.astype(np.uint8),         # 1 = keep, 0 = interpolate
            max_search_distance=max_search_distance,
            smoothing_iterations=smoothing_iterations,
        )
        after = (filled == NODATA).mean()        # residual NoData after fill

    profile.update(dtype="float32", nodata=NODATA, count=1)
    with rasterio.open(dst_path, "w", **profile) as dst:
        dst.write(filled.astype(np.float32), 1)  # preserves CRS + transform

    return before - after                        # fraction of cells recovered

The mask argument follows rasterio’s convention where 1/True marks cells to keep and 0 marks cells to interpolate — inverting it fills the wrong pixels and destroys your real surface, so keep the valid = band != NODATA line exactly as written. Leave smoothing_iterations at 0 for survey DSMs where you must not alter measured elevations; raise it to 1–2 only for visualization products where a smoother filled patch matters more than metric fidelity. If you would rather stay inside the GDAL toolchain, the command-line equivalent is gdal_fillnodata.py -md 12 dsm.tif dsm_filled.tif, which uses the same algorithm with -md as the pixel search distance.

Figure 2 — Fill methods differ in how they are wrong, not in whether they are. Choosing one is choosing which artefact you would rather explain to whoever measures volumes from the result.

# Edge-case matrix

The lake case is the one that bites: a large physical void filled with a bounded search stays honestly empty, but an unbounded fill produces a plausible-looking ramp that no one notices until a volume report is wrong. When water covers a meaningful area, digitise a mask and hold those cells at NoData deliberately rather than trusting the distance bound alone.

Figure 3 — The fill is defensible only if it is labelled. Writing the mask in the same pass costs almost nothing and is the difference between a surface with known provenance and a surface with a story.

# Verification snippet

Confirm the fill closed the sampling voids without touching the elevation of cells that already had data — a fill that shifts valid pixels means the mask convention was inverted.

import numpy as np
import rasterio

NODATA = -9999.0
with rasterio.open("dsm.tif") as a, rasterio.open("dsm_filled.tif") as b:
    orig, fill = a.read(1), b.read(1)

valid_before = orig != NODATA
# 1. Valid cells must be untouched by the fill.
assert np.allclose(orig[valid_before], fill[valid_before]), "fill altered real data"
# 2. NoData fraction must drop.
before = (~valid_before).mean()
after = (fill == NODATA).mean()
assert after < before, "fill recovered no cells — check the mask convention"
print(f"NoData {before:.2%} -> {after:.2%}, recovered {before - after:.2%}")

If assertion 1 fails, mask was inverted; if assertion 2 fails with no change, either there were no voids within max_search_distance or the sentinel did not match the stored NoData value.

Finally, report the filled fraction with the deliverable. A DSM that is 2% interpolated and one that is 25% interpolated are different products, and only the metadata distinguishes them once the surface is continuous.

# When to escalate

Filling is cosmetic repair of a sampling problem; it cannot recover a surface the reconstruction never built. Return to generating DSM and DTM from point clouds with PDAL and reprocess when:

• Voids cover a large, structured region (a whole flight strip), which points to a coverage or overlap failure upstream, not a rasterization gap — re-fly or re-match rather than fill.
• The holes trace real texture-poor surfaces you must measure (water level, glass roof height); no interpolation is metrically valid there, so capture control on those surfaces or accept NoData.
• Widening radius alone would close the voids, in which case fix it at the source in writers.gdal and re-rasterize, because a correctly sampled grid is always better than a filled one.

# Related

• Generating DSM and DTM from Point Clouds with PDAL
• Removing Vegetation from DTM Ground Surfaces
• Exporting Cloud-Optimized GeoTIFFs with Rasterio

← Generating DSM and DTM from Point Clouds with PDAL