Filling Voids in DEM Tiles

Fill small voids by interpolation, large ones from a secondary DEM, and record what was filled:

from scipy import ndimage
import numpy as np

void = ~np.isfinite(dem)
labels, n = ndimage.label(void)
sizes = ndimage.sum(void, labels, range(1, n + 1))
small = np.isin(labels, np.flatnonzero(sizes <= 50) + 1)    # up to ~50 cells: interpolate
large = void & ~small                                        # the rest: secondary source

Filling is not optional before deriving terrain products: a single unfilled cell becomes a ring of undefined slope around it. This page belongs to terrain analysis and DEM-derived products in Satellite Processing Workflows & Index Pipelines.

Why Void Size Decides the Method

Fill methods across two void sizes Across a small void of a few cells, nearest-valid and interpolated fills are nearly indistinguishable. Across a large void spanning a valley, nearest-valid fill produces two flat terraces meeting in a step at the centre, which becomes a false cliff in slope; interpolation produces a smooth surface; a secondary DEM restores the real valley shape. small void: any method works large void: method matters void void over a valley nearest-valid interpolated secondary DEM

Small voids are dominated by their surroundings, so any sensible fill converges on the same answer. Large voids hide real terrain — a valley, a ridge — that no interpolation can reconstruct from the edges alone, and the only faithful fill is elevation measured by a different sensor. Choosing the method by void size is therefore the whole strategy.

Environment & Setup

Package Version pin Used for
numpy >=1.23 Masks and arithmetic
scipy >=1.11 Labelling voids, distance transforms, interpolation
rioxarray >=0.15 Aligning a secondary DEM onto the primary grid
rasterio >=1.3.0 Reading and writing, including fillnodata 
pip install "numpy>=1.23" "scipy>=1.11" "rioxarray>=0.15" "rasterio>=1.3.0"

Complete Working Example

import numpy as np
import rioxarray
from rasterio.fill import fillnodata
from scipy import ndimage


def fill_dem(primary_path: str, secondary_path: str | None = None, *,
             small_max_cells: int = 50, blend_px: int = 5):
    dem = rioxarray.open_rasterio(primary_path, masked=True).squeeze("band")
    z = dem.values.astype("float32")
    void = ~np.isfinite(z)
    filled_mask = np.zeros(z.shape, dtype="uint8")          # 0 measured, 1 interpolated, 2 secondary

    labels, n = ndimage.label(void)
    sizes = ndimage.sum(void, labels, range(1, n + 1)) if n else np.array([])
    small = np.isin(labels, np.flatnonzero(sizes <= small_max_cells) + 1)
    large = void & ~small

    # 1. Large voids from a secondary DEM, offset-corrected and blended at the edge
    if secondary_path and large.any():
        sec = rioxarray.open_rasterio(secondary_path, masked=True).squeeze("band")
        s = sec.rio.reproject_match(dem, resampling=1).values.astype("float32")
        ring = ndimage.binary_dilation(large, iterations=blend_px) & ~void
        bias = float(np.nanmedian(z[ring] - s[ring])) if ring.any() else 0.0
        s = s + bias                                        # remove vertical offset between sources
        z = np.where(large & np.isfinite(s), s, z)
        filled_mask[large & np.isfinite(s)] = 2

    # 2. Remaining voids by inverse-distance interpolation from the edges
    remaining = ~np.isfinite(z)
    if remaining.any():
        z = fillnodata(np.nan_to_num(z, nan=0.0), mask=(~remaining).astype("uint8"),
                       max_search_distance=100, smoothing_iterations=0)
        filled_mask[remaining & (filled_mask == 0)] = 1

    return dem.copy(data=z), filled_mask

Correcting the vertical bias between the two sources before patching is what prevents a visible step around every filled void. Different DEMs have different datums, different surface-versus-terrain characteristics and different systematic errors, so their elevations over the same ground commonly differ by a few metres; the median difference on a ring of valid cells around the void measures that offset locally. The same principle — calibrate on the overlap before combining — appears in harmonising Landsat and Sentinel-2 reflectance.

Where Voids Cluster

Radar voids fall on the steepest back slopes A side-looking radar illuminates slopes facing it well and receives no return from steep slopes facing away, which fall in radar shadow. The resulting DEM voids concentrate on exactly the steepest mountain terrain, where slope, aspect and illumination correction matter most. Voids are not random radar red: radar shadow, no return, void in the DEM Fill quality matters most exactly where the terrain is most extreme.

Because voids concentrate on steep, shadowed slopes, a poor fill damages exactly the terrain where slope and illumination correction are most consequential. That is the practical argument for a secondary-DEM fill in mountainous areas even when the voids are modest in number: interpolating across a steep back-slope smooths away the very gradient that makes it steep.

Water bodies are the other common void class. For them the right fill is usually a flat surface at the water level, taken from the shoreline elevation, rather than any interpolation — which would otherwise produce a bowl or a dome where the surface is in fact level.

Verification

import numpy as np

filled, mask = fill_dem("dem_30m.tif", "dem_secondary.tif")
assert np.isfinite(filled.values).all(), "voids remain"
print(f"measured {np.mean(mask == 0):.1%}, interpolated {np.mean(mask == 1):.2%}, "
      f"secondary {np.mean(mask == 2):.2%}")

from scipy import ndimage
edge = ndimage.binary_dilation(mask > 0, iterations=1) & (mask == 0)
step = np.abs(ndimage.uniform_filter(filled.values, 3) - filled.values)[edge]
assert np.percentile(step, 99) < 10, "large steps at fill edges — check the vertical bias"
Ship the mask with the fill The filled DEM is stored with a companion mask coding each cell as measured, interpolated or taken from a secondary source. Downstream products can then exclude filled cells, weight them lower, or at least report what share of an area rests on estimates rather than measurements. Filled DEM plus provenance mask filled elevation 1 2 0 measured · 1 interpolated · 2 secondary Without the mask, nobody downstream can tell an estimate from a measurement.

The edge-step check is the one that catches a missed bias correction: filled regions sitting a few metres above or below their surroundings show as a ring of large local differences at the boundary. Rendering a hillshade over the filled areas makes the same defect visible as a crisp outline around each patch.

Common Errors

Filled areas show as flat terraces

Nearest-valid filling was used on large voids. Use a secondary DEM or interpolation for anything beyond a few cells.

A step appears around every patched void

The vertical offset between the two DEMs was not corrected. Measure it on a ring around each void and subtract.

Lakes become bowls

Water voids were interpolated. Fill them flat at the shoreline elevation instead, taken as the median elevation of the cells bordering the water body.

Slope layer still has holes

Filling ran after slope was computed. Fill first, then derive, and mask the cells adjacent to any remaining void so the gradient never touches missing data.

Frequently Asked Questions

Q: Why do DEMs have voids at all? Radar-derived DEMs lose signal in steep terrain where slopes face away from the sensor, and over water and sand where the return is weak. Optical DEMs lose cells under persistent cloud and over featureless surfaces. The voids cluster exactly where terrain is most interesting.

Q: Is nearest-valid filling good enough? For voids a few cells across, yes. For larger voids it produces flat terraces with abrupt edges, which become false cliffs in slope and hillshade. Use interpolation or a secondary DEM for anything bigger.

Q: Should filled cells be flagged? Always. Filled elevation is an estimate of lower quality than measured elevation, and any product derived from it inherits that. A mask of filled cells lets downstream users exclude or down-weight them.

Q: Which secondary DEM should fill a primary one? The best available product with independent errors — ideally from a different sensor type, so its voids fall elsewhere. A radar DEM filled from an optical one, or the reverse, covers most of each other’s gaps.

Q: Can voids be filled from the imagery itself? Only indirectly. Stereo imagery can produce elevation, but that is building a new DEM rather than filling one. For void filling, an existing independent DEM is almost always the practical source.

Q: How should the filled-cell mask be distributed? As a second band in the same file or as a sidecar raster on the same grid, with its codes documented in the tags. A second band keeps the two inseparable, which is usually what you want; a sidecar keeps the elevation file compatible with tools that expect a single band.