Detecting Voids and Nodata Gaps in DEMs

A void is a hole where the survey claims coverage: cells with no elevation value inside the area the producer says was flown. They come from water bodies that absorbed the lidar pulse, from shadowed slopes in a photogrammetric match, from clouds, from failed returns over glass roofs, and from processing steps that dropped cells nobody re-checked. This guide finds them, measures them the way a contract can be judged against, turns them into reviewable polygons, and decides which ones may be filled — the operational half of the checks described in Raster and Elevation Data Quality Checks.

Prerequisites

• GDAL 3.4+ with the Python bindings, or rasterio 1.3+ plus numpy 1.24+. The examples use rasterio for reading and rasterio.features for mask polygonisation; the gdal_fillnodata utility appears once, in the fill step.
• A declared nodata value on the raster. If the file does not declare one, stop: the void measurement is not defined, and that missing declaration is the finding to report.
• The product tolerance, expressed as two numbers — maximum total void percentage inside the footprint, and maximum area of any single void. Without both, “how many voids are too many” is decided per reviewer.
• Enough disk for a mask, roughly one bit per cell if written as a packed mask, or one byte per cell as an unpacked GeoTIFF. On a 40,000 × 40,000 tile that is 200 MB unpacked, so write masks compressed.

Step-by-Step Procedure

Step 1 — Confirm the nodata declaration before measuring anything

# voids/step1_declaration.py
import numpy as np
import rasterio


def nodata_declaration(path: str) -> dict:
    """Report how the raster represents missing data, and whether it is consistent."""
    with rasterio.open(path) as src:
        declared = src.nodata
        sample = src.read(1, out_shape=(min(src.height, 2048), min(src.width, 2048)))

    floaty = np.issubdtype(sample.dtype, np.floating)
    has_nan = bool(np.isnan(sample).any()) if floaty else False
    sentinel_hits = int((sample == declared).sum()) if declared is not None else 0

    return {
        "declared_nodata": declared,
        "dtype": str(sample.dtype),
        "nan_present": has_nan,
        "sentinel_cells_in_sample": sentinel_hits,
        "mixed_sentinels": bool(has_nan and declared is not None and sentinel_hits > 0),
    }

Verification: a healthy elevation tile reports a declared sentinel and mixed_sentinels: false. If declared_nodata is None but the sample contains -9999 or a large negative plateau, the file has an undeclared sentinel — report it as a blocker and do not attempt a void measurement, because every consumer will read those cells as elevations.

Step 2 — Derive the data footprint from the valid-data mask

The footprint is what makes a void measurable. Deriving it from the data itself avoids inheriting an optimistic extent from the delivery paperwork.

# voids/step2_footprint.py
import rasterio
from rasterio.features import shapes
from shapely.geometry import shape
from shapely.ops import unary_union


def data_footprint(path: str, simplify_m: float = 2.0):
    """Polygonise the valid-data mask and return its outer footprint."""
    with rasterio.open(path) as src:
        mask = src.dataset_mask()          # 255 where valid, 0 where nodata
        transform = src.transform

    polys = [
        shape(geom)
        for geom, value in shapes(mask, mask=(mask == 255), transform=transform)
        if value == 255
    ]
    if not polys:
        return None

    footprint = unary_union(polys)
    # Keep only the outer boundary: interior rings here ARE the voids.
    if footprint.geom_type == "Polygon":
        outer = shape({"type": "Polygon", "coordinates": [list(footprint.exterior.coords)]})
    else:
        outer = unary_union([
            shape({"type": "Polygon", "coordinates": [list(p.exterior.coords)]})
            for p in footprint.geoms
        ])
    return outer.simplify(simplify_m, preserve_topology=True)

Verification: the footprint area should be within a few percent of the tile area for an interior tile, and noticeably smaller for an edge tile. Note the comment on interior rings — dropping them is deliberate, because those holes are exactly what the next step counts.

Step 3 — Measure void coverage inside the footprint

# voids/step3_measure.py
import numpy as np
import rasterio
from rasterio.features import geometry_mask


def measure_voids(path: str, footprint) -> dict:
    """Void cells as a share of footprint cells, computed block by block."""
    inside_total = 0
    void_total = 0

    with rasterio.open(path) as src:
        nodata = src.nodata
        for _ij, window in src.block_windows(1):
            block = src.read(1, window=window)
            tf = src.window_transform(window)

            inside = ~geometry_mask([footprint], out_shape=block.shape,
                                    transform=tf, invert=False)
            if not inside.any():
                continue

            missing = np.isnan(block) if np.issubdtype(block.dtype, np.floating) \
                else np.zeros(block.shape, dtype=bool)
            if nodata is not None:
                missing |= (block == nodata)

            inside_total += int(inside.sum())
            void_total += int((missing & inside).sum())

    with rasterio.open(path) as src:
        cell_area = abs(src.transform.a * src.transform.e)

    return {
        "footprint_cells": inside_total,
        "void_cells": void_total,
        "void_pct": round(100.0 * void_total / inside_total, 4) if inside_total else 0.0,
        "void_area_m2": round(void_total * cell_area, 2),
    }

Verification: run against a tile you have inspected visually. The void area in square metres should match a rough manual estimate from the rendered image; if it is orders of magnitude out, the footprint is wrong or the sentinel test is missing one of the two representations from Step 1.

Step 4 — Vectorise voids so each hole is reviewable

A percentage tells you whether the tile passes. Polygons tell the producer what to re-fly.

# voids/step4_vectorise.py
import geopandas as gpd
import numpy as np
import rasterio
from rasterio.features import shapes
from shapely.geometry import shape


def void_polygons(path: str, footprint, min_area_m2: float = 1.0) -> gpd.GeoDataFrame:
    """One polygon per void inside the footprint, with its area."""
    with rasterio.open(path) as src:
        data = src.read(1)
        nodata = src.nodata
        transform = src.transform
        crs = src.crs

    missing = np.isnan(data) if np.issubdtype(data.dtype, np.floating) \
        else np.zeros(data.shape, dtype=bool)
    if nodata is not None:
        missing |= (data == nodata)

    geoms = [
        shape(geom)
        for geom, value in shapes(missing.astype("uint8"), mask=missing, transform=transform)
        if value == 1
    ]
    gdf = gpd.GeoDataFrame(geometry=geoms, crs=crs)
    if gdf.empty:
        return gdf

    gdf = gdf[gdf.intersects(footprint)].copy()
    gdf["area_m2"] = gdf.area.round(2)
    gdf = gdf[gdf["area_m2"] >= min_area_m2]
    return gdf.sort_values("area_m2", ascending=False).reset_index(drop=True)

Verification: the summed area_m2 of the returned polygons should match void_area_m2 from Step 3 to within the min_area_m2 filter. A large discrepancy means single-cell voids dominate — a scattered speckle pattern, which is a different defect from a few large holes and calls for a different response.

Step 5 — Fill only what the specification allows

# voids/step5_fill.py — fill small voids, record what was filled
import subprocess


def fill_small_voids(src_path: str, out_path: str, mask_path: str, max_distance_cells: int = 8):
    """Interpolate across narrow voids and write a mask of the cells that were filled."""
    # 1. Record where the voids are BEFORE filling — this mask ships with the product.
    subprocess.run(
        ["gdal_calc.py", "-A", src_path, "--outfile", mask_path,
         "--calc", "A==A", "--NoDataValue", "0", "--type", "Byte",
         "--co", "COMPRESS=DEFLATE"],
        check=True,
    )
    # 2. Interpolate only across voids narrower than max_distance_cells.
    subprocess.run(
        ["gdal_fillnodata.py", "-md", str(max_distance_cells), "-si", "1",
         src_path, out_path],
        check=True,
    )
    return {"filled": out_path, "fill_mask": mask_path,
            "max_distance_cells": max_distance_cells}

Verification: re-run Step 3 on the filled raster. The remaining void percentage should be the share attributable to voids wider than the fill distance — those are the ones a human must decide about. Ship the fill mask: a downstream hydrological model must be able to exclude interpolated cells, and it cannot do that if the fill is invisible.

Interpreting Results

The two numbers that decide acceptance are void percentage inside footprint and largest single void area, and they fail for different reasons. A high percentage made of thousands of single-cell voids points at sensor noise or an over-aggressive filter in production; a low percentage containing one 4,000 square metre hole points at a specific failure — a cloud, a lake, an aborted flight line — that interpolation must not paper over.

The polygon output tells you which. Sort the voids by area and look at the distribution: a long tail of tiny voids with nothing above a few square metres is a speckle problem, fixable in production settings. A handful of large voids with clean boundaries is a coverage problem, and no amount of interpolation makes the missing survey exist.

Location matters as much as size. A 300 square metre void in the middle of a floodplain changes a modelled water surface; the same void on a hilltop changes nothing anyone will measure. Intersect the void polygons with whatever the product is used for — hydrological network, building footprints, design corridors — and let the intersection drive escalation rather than the raw area alone.

Gotchas & Edge Cases

Water bodies are voids in lidar and features in the specification. A lake surface returns no usable pulse, so the raw model has a hole exactly where the lake is. Many products then set a constant water level. Both are legitimate; what is not legitimate is a pipeline that treats the two identically. Mask known hydrography before measuring, or the same tile passes or fails depending on which production step it came from.

dataset_mask() is not the same as testing the nodata value. Rasterio’s mask honours an alpha band or an internal mask if one exists, and falls back to the nodata value otherwise. On a file with both an internal mask and a stale nodata value, the two disagree. Decide which is authoritative for your programme and use it consistently — mixing them across tiles produces a void statistic that drifts for no data reason.

Filling changes statistics, so measure first. Elevation minimum, maximum and percentiles all move after interpolation. Compute the plausibility statistics before the fill, or you validate the interpolation rather than the survey.

Single-cell voids can be an artefact of the sentinel value. In an integer terrain model, a genuine elevation of exactly -9999 is impossible, but a genuine value of 0 is not. Products that use 0 as nodata are unreliable to measure; report the sentinel choice as a defect rather than working around it.

Compression hides nothing but slows scans. A DEFLATE-compressed tile still needs decompressing to count nodata, so the void scan cost scales with the compressed block layout. Striped rasters make windowed reads pathological — convert to a tiled layout first, as covered in Checking Cloud-Optimized GeoTIFF Structure.

When to Escalate

Send the tile back to the producer, rather than repairing it locally, when:

• A single void exceeds the specification’s maximum area. Interpolation across a large hole invents terrain, and the invention is indistinguishable from measurement once it is in the file.
• Voids align with flight lines or tile edges. Systematic geometry means a production failure, and re-processing the source is cheaper than patching every affected tile.
• The nodata sentinel is undeclared or mixed. This is a delivery defect. Fixing it downstream means every consumer must know your local convention, which defeats the point of a declared sentinel.
• Void locations coincide with the product’s purpose. Holes in the drainage network of a flood model are a coverage failure regardless of what the total percentage says.

For everything else — small, scattered, interpolable voids inside tolerance — fill under the rules in Step 5, ship the fill mask, and record the fill in the run report so the next audit can see what was interpolated and why. The severity assignment follows the model in Classifying Topology Errors by Severity.

Frequently Asked Questions

Related

• Raster and Elevation Data Quality Checks — the full raster check set this void measurement belongs to
• Validating Raster Alignment and Resolution with Rasterio — the grid checks that run before any pixel is read
• Checking Cloud-Optimized GeoTIFF Structure — why a tiled layout makes this scan affordable

Back to Raster and Elevation Data Quality Checks