Raster and Elevation Data Quality Checks

A digital elevation model (DEM) that fails quietly is more dangerous than a vector layer that fails loudly. Nothing raises an exception when a tile is shifted half a pixel, when a void sits in the middle of a floodplain, or when a product mixes ellipsoidal and orthometric heights — the raster opens, renders, and produces plausible contours that are wrong by a metre. This topic covers the checks that catch those defects automatically: grid and transform profiling, void accounting inside the footprint, horizontal and vertical reference validation, tile alignment, and the overview and compression checks that decide whether a raster is usable at scale. It sits alongside the vector work in Core Spatial QC Fundamentals & Standards and feeds the same severity model used everywhere else in the pipeline.

Prerequisites

• GDAL 3.4+ and rasterio 1.3+ in the validation environment. Confirm with gdalinfo --version and python -c "import rasterio; print(rasterio.__version__)". Rasterio 1.3 is the first release with stable windowed reads over /vsicurl/ sources, which matters when the rasters live in object storage.
• numpy 1.24+ for the array statistics, and optionally rio-cogeo 5.x if you validate cloud-optimized GeoTIFF (COG) structure.
• A declared product specification. Every check below compares the file against an expectation: the target pixel size, the grid origin, the horizontal coordinate reference system (CRS), the vertical datum, the nodata value and the acceptable void percentage. Without a written specification these checks have nothing to assert against.
• The data footprint, either supplied as a polygon or derivable from the valid-data mask. Void accounting is meaningless without it, because nodata outside the footprint is correct behaviour.
• A canonical CRS decision for the programme, already made in Coordinate Reference System Precision Standards. Raster reprojection is lossy in a way vector reprojection is not — resampling changes pixel values — so the target CRS should be chosen once and enforced at ingestion.

Conceptual Foundation

A raster is three things at once, and each carries its own class of defect. It is an array of values with a data type and a nodata sentinel; it is a grid positioned in the world by an affine transform; and it is a product with a specification stating what the values mean. Most raster QC failures are a disagreement between those three layers rather than corruption of the pixels themselves.

The affine transform is where the subtle errors live. Six numbers place the array in the world: pixel width, row rotation, upper-left x, column rotation, pixel height (negative for a north-up raster), and upper-left y. A tile whose origin is not an exact multiple of its pixel size is not aligned to the product grid, and every downstream mosaic, zonal statistic or tile join will resample it — introducing a half-pixel blur that nobody attributes to the transform three months later. Checking alignment costs two modulo operations against the declared grid origin.

Nodata is the second recurring source of confusion. The value is a sentinel, not a measurement, and it means “no observation here” — which is correct outside a survey footprint and a defect inside it. Distinguishing the two requires the footprint, so a void check is really a mask intersection: count cells that are nodata and inside the footprint, divide by the footprint cell count, and compare that ratio against the product tolerance. A raster reported as “12% nodata” tells you nothing; the same raster reported as “0.4% void inside footprint” is an acceptance decision.

Elevation adds a third axis and a third reference system. A height is meaningless without knowing whether it is measured from an ellipsoid or from a geoid model — the two differ by tens of metres in many parts of the world, and by enough to change a flood extent almost everywhere. The vertical datum is declared in a compound CRS when the producer bothers, and absent otherwise. Any product used for drainage, flood modelling or engineering design should treat a missing vertical datum the same way a cadastral pipeline treats a missing horizontal one: as a blocker, routed to the same dead-letter path described in Categorizing and Prioritizing Spatial Errors.

Step-by-Step Implementation

Step 1 — Profile the grid from the header

Every contract question that can be answered from the header should be answered there, before a single pixel is read. This is the cheapest stage in raster QC and it catches the majority of publishing mistakes.

# raster_qc/profile.py — header-only grid profile
import rasterio
from rasterio.crs import CRS


def profile_raster(path: str) -> dict:
    """Read the contract-relevant properties of a raster without touching pixel data."""
    with rasterio.open(path) as src:
        t = src.transform
        return {
            "path": path,
            "driver": src.driver,
            "width": src.width,
            "height": src.height,
            "band_count": src.count,
            "dtype": src.dtypes[0],
            "nodata": src.nodata,
            "crs": src.crs.to_string() if src.crs else None,
            "epsg": src.crs.to_epsg() if src.crs else None,
            "pixel_size": (abs(t.a), abs(t.e)),
            "origin": (t.c, t.f),
            "rotation": (t.b, t.d),
            "overview_levels": src.overviews(1),
            "block_shape": src.block_shapes[0],
        }


if __name__ == "__main__":
    import json
    import sys
    print(json.dumps(profile_raster(sys.argv[1]), indent=2, default=str))

Verification: run the profiler against a known-good tile. rotation must be (0.0, 0.0) for a north-up product, nodata must not be None, and epsg must equal your canonical code. A None nodata on a product that clearly has empty margins means the sentinel is undeclared — every consumer will read those cells as real values, often 0, and a zero-elevation cliff will appear at the tile edge.

Step 2 — Assert the grid against the product specification

The profile is data; the specification turns it into a pass or fail. Keep the specification in version control next to the rules, as described in Defining Spatial Data Quality Policies.

# raster_qc/rules.py — assert a profile against the declared product spec
from dataclasses import dataclass


@dataclass(frozen=True)
class RasterSpec:
    epsg: int
    pixel_size: float          # metres, square pixels
    grid_origin: float         # product grid must be a multiple of this, in CRS units
    dtype: str
    require_nodata: bool = True
    require_overviews: bool = True


def check_profile(profile: dict, spec: RasterSpec) -> list[dict]:
    findings: list[dict] = []

    def fail(rule: str, severity: str, message: str) -> None:
        findings.append({"rule": rule, "severity": severity, "message": message})

    if profile["epsg"] != spec.epsg:
        fail("RAS_CRS_001", "blocker",
             f"declared EPSG:{profile['epsg']} != canonical EPSG:{spec.epsg}")

    px, py = profile["pixel_size"]
    if abs(px - spec.pixel_size) > 1e-6 or abs(py - spec.pixel_size) > 1e-6:
        fail("RAS_RES_001", "blocker",
             f"pixel size {px}x{py} != specified {spec.pixel_size}")

    ox, oy = profile["origin"]
    off_x = ox % spec.grid_origin
    off_y = oy % spec.grid_origin
    if min(off_x, spec.grid_origin - off_x) > 1e-6 or min(off_y, spec.grid_origin - off_y) > 1e-6:
        fail("RAS_ALIGN_001", "blocker",
             f"origin ({ox}, {oy}) is not on the {spec.grid_origin} m product grid")

    if profile["rotation"] != (0.0, 0.0):
        fail("RAS_ROT_001", "blocker", "raster is rotated; product requires a north-up grid")

    if spec.require_nodata and profile["nodata"] is None:
        fail("RAS_NODATA_001", "blocker", "no nodata value declared")

    if profile["dtype"] != spec.dtype:
        fail("RAS_TYPE_001", "warning",
             f"dtype {profile['dtype']} != specified {spec.dtype}")

    if spec.require_overviews and not profile["overview_levels"]:
        fail("RAS_OVR_001", "warning", "no overviews present; clients will read full resolution")

    return findings

Verification: feed the checker a deliberately shifted tile — add half a pixel to the origin — and confirm RAS_ALIGN_001 fires. Alignment is the rule most often missing from raster pipelines, and the one whose absence is hardest to notice after the fact.

Step 3 — Measure voids inside the footprint

This is the one check that needs pixel data, and it should be windowed so memory stays flat regardless of raster size.

# raster_qc/voids.py — void accounting inside a declared footprint
import numpy as np
import rasterio
from rasterio.features import geometry_mask


def void_stats(path: str, footprint_geom: dict) -> dict:
    """Fraction of the footprint that is nodata, computed window by window."""
    footprint_cells = 0
    void_cells = 0

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

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

            missing = np.isnan(data) if np.issubdtype(data.dtype, np.floating) and nodata is None \
                else (data == nodata)

            footprint_cells += int(inside.sum())
            void_cells += int((missing & inside).sum())

    pct = (100.0 * void_cells / footprint_cells) if footprint_cells else 0.0
    return {
        "footprint_cells": footprint_cells,
        "void_cells": void_cells,
        "void_pct": round(pct, 4),
    }

Verification: run it twice — once with the true footprint and once with a footprint clipped to the data extent. The two void percentages should be close; a large divergence means the supplied footprint is wrong, not that the raster is bad. Void tolerance is a product decision, typically 0.1% for photogrammetric surfaces and up to 2% for lidar-derived models in dense vegetation.

Step 4 — Check statistics for physical plausibility

Statistics answer a different question from structure: not “is the file correct” but “are these values possible”. They are cheap when read from overviews.

# raster_qc/stats.py — plausibility bounds on an elevation surface
import numpy as np
import rasterio


def elevation_stats(path: str, min_valid: float, max_valid: float) -> dict:
    """Summary statistics from the coarsest overview, with plausibility flags."""
    with rasterio.open(path) as src:
        levels = src.overviews(1)
        factor = levels[-1] if levels else 1
        data = src.read(
            1, out_shape=(src.height // factor or 1, src.width // factor or 1),
        ).astype("float64")
        if src.nodata is not None:
            data[data == src.nodata] = np.nan

    valid = data[~np.isnan(data)]
    if valid.size == 0:
        return {"empty": True}

    stats = {
        "min": float(valid.min()),
        "max": float(valid.max()),
        "mean": float(valid.mean()),
        "p01": float(np.percentile(valid, 1)),
        "p99": float(np.percentile(valid, 99)),
        "sampled_from_overview": factor,
    }
    stats["below_min_valid"] = stats["min"] < min_valid
    stats["above_max_valid"] = stats["max"] > max_valid
    # A spike of exactly 0 in a terrain model usually means undeclared nodata.
    stats["suspect_zero_plateau"] = bool(((valid == 0).mean() > 0.02) and min_valid > 0)
    return stats

Verification: on a national terrain model the 1st and 99th percentiles should sit well inside the country’s real elevation range. A minimum of -9999 proves an undeclared nodata sentinel; a maximum of 32767 proves an integer overflow during a format conversion.

Step 5 — Compare adjacent tiles at their shared edge

Tiles are validated individually and consumed as a mosaic, so the seam is a first-class check.

# raster_qc/seams.py — compare the shared edge of two adjacent tiles
import numpy as np
import rasterio


def edge_mismatch(left_path: str, right_path: str, tolerance: float = 0.05) -> dict:
    """Mean absolute difference along the shared column of two horizontally adjacent tiles."""
    with rasterio.open(left_path) as left, rasterio.open(right_path) as right:
        if left.transform.a != right.transform.a:
            return {"comparable": False, "reason": "different pixel sizes"}
        if abs((left.bounds.right) - (right.bounds.left)) > 1e-6:
            return {"comparable": False, "reason": "tiles are not adjacent"}

        rows = min(left.height, right.height)
        left_edge = left.read(1, window=((0, rows), (left.width - 1, left.width)))
        right_edge = right.read(1, window=((0, rows), (0, 1)))

    a = left_edge.astype("float64").ravel()
    b = right_edge.astype("float64").ravel()
    diff = np.abs(a - b)
    return {
        "comparable": True,
        "mean_abs_diff": float(np.nanmean(diff)),
        "max_abs_diff": float(np.nanmax(diff)),
        "exceeds_tolerance": bool(np.nanmean(diff) > tolerance),
    }

Verification: neighbouring tiles from a single production run should differ by well under the vertical accuracy of the product. A systematic offset of a constant value across the whole seam is a datum or calibration difference between production batches, not noise.

Common Failure Modes & Fixes

Performance & Scale Considerations

Header reads are constant time and effectively free — a full grid profile over ten thousand tiles takes minutes, and it answers the CRS, resolution, alignment, nodata and overview questions for the whole collection. Design the pipeline so this stage runs on every tile, every time, and only failures escalate to a pixel read.

Pixel-level checks should always be windowed. Rasterio’s block_windows iterator yields the raster’s native tiling, which means each read touches exactly one compressed block and memory stays bounded by the block size rather than the raster size. A 40 GB elevation mosaic validates in constant memory this way; the same check written as src.read(1) fails on any machine you would want to pay for.

Statistical checks belong on the overviews. A level-8 overview is 1/64th the pixel count and answers plausibility questions with more than enough fidelity — provided the overview level is recorded with the numbers, because comparing statistics sampled at different levels produces drift that looks like data change.

Parallelism is per tile and needs no coordination, which puts raster QC firmly in the well-behaved half of the scaling model described in Batch Processing Large Spatial Datasets. The one exception is the seam check, which needs both neighbours — partition by tile-pair rather than by tile, or run it as a second pass over an index of adjacencies.

Integration with the Validation Pipeline

Raster checks slot into the same DAG as vector rules, at the ingestion and rule-evaluation stages defined in Validation Pipeline Architecture. The header profile runs at ingestion, where a wrong CRS or an unaligned grid should stop the tile before anything expensive happens. Void and statistical checks run in the rule stage, emitting the same result contract as every other rule — feature identifier (here, the tile identifier), rule identifier, severity, message and, where useful, the geometry of the defect as a footprint polygon of the void.

Because the results share a contract, they share the report, the dashboard and the alerting thresholds. A raster void rate belongs on the same quality scorecard as a vector overlap count, and both should feed the observability signals described in Observability and Lineage for Validation.

Frequently Asked Questions

Related

• Detecting Voids and Nodata Gaps in DEMs — mapping and measuring the holes inside a survey footprint with GDAL
• Validating Raster Alignment and Resolution with Rasterio — grid origin, pixel size and tile-seam checks in code
• Checking Cloud-Optimized GeoTIFF Structure — tiling, overviews and header layout for remote reads
• Coordinate Reference System Precision Standards — the canonical CRS decision these checks assert against
• Geometry Validity Checks for Vector Data — the vector half of the same fundamentals

Back to Core Spatial QC Fundamentals & Standards