Validating Nodata and Band Integrity After Warp

Reprojecting a raster is the one operation in a coverage pipeline that can produce a technically perfect output which is completely wrong. The warp succeeds, the file is valid, every tile renders, the pyramid is complete — and the imagery is in the wrong bands, or three quarters of the scene has quietly become nodata because the target grid clipped it, or the nodata value itself was dropped so that a sentinel like -9999 is now being averaged into overviews as if it were an elevation. None of these raise an error. This guide covers the assertions to run immediately after a warp so that a coverage which changed shape is rejected rather than published, and it belongs to raster and tile pipeline observability under geospatial observability architecture fundamentals.

Three silent warp defects compared against a correct output Four panels show a source scene and three defective warp outputs. The correct output preserves the scene footprint, three bands in the order red green blue, and a declared nodata sentinel. The first defect shows nodata expansion, where the target grid clips the scene so that a large border region becomes nodata. The second shows band reordering, where the band order becomes blue green red and the imagery renders in wrong colours. The third shows a dropped nodata declaration, where the sentinel value survives as data and contaminates overview statistics. Each panel names the assertion that catches it. Each defect produces a valid file — and a different measurable signature valid pixels correct nodata 4% · bands RGB sentinel declared nodata expansion nodata 4% → 79% caught by nodata_ratio delta B G R band reordering count unchanged, order wrong caught by per-band statistics min = −9999 treated as elevation sentinel dropped nodata_ratio reads 0% caught by range assertion A ratio of zero is as suspicious as a ratio of eighty percent — assert both bounds.

Problem framing: why a warp changes shape silently

Reprojection is a resampling operation onto a new grid, and three things can shift without producing a failure.

Footprint against grid. The target grid is defined by an extent and a resolution. If the source scene does not tile neatly onto it — because the extent was inherited from a different scene, or the resolution was rounded — the warp fills the difference with nodata. The output is a valid raster of the requested shape; it simply contains less data than it should. This is by far the most common silent defect, and it is invisible unless you compare nodata ratios before and after.

Band identity. Band order is positional, not named, in most raster formats. A source that changes band order, or a warp invoked with an explicit band list that no longer matches, produces a file with the correct band count and the wrong band content. Rendering succeeds; colours are wrong; nothing counts an error.

Nodata declaration. The nodata value is metadata, not data. If it is not carried through the warp, the sentinel value survives as an ordinary sample. Overviews then average it in, statistics are poisoned, and a terrain layer acquires a plateau at minus nine thousand nine hundred and ninety-nine metres. The nodata ratio reads zero, which looks like the healthiest possible value.

All three are contract violations against the source, which means the right check compares the output to the input, not to an absolute threshold.

Implementation: assert the output against the source

Run the assertions in the warp worker, before the output is promoted. Each is cheap relative to the warp itself.

# warp_assertions.py — reject a coverage whose shape changed during reprojection.
from dataclasses import dataclass
from osgeo import gdal
from opentelemetry import metrics

gdal.UseExceptions()
meter = metrics.get_meter("gis.raster")
nodata_ratio = meter.create_observable_gauge("gis.raster.nodata_ratio", callbacks=[])
warp_rejected = meter.create_counter("gis.raster.warp_rejected_total")


@dataclass
class BandProfile:
    index: int
    nodata: float | None
    valid_ratio: float
    minimum: float
    maximum: float
    mean: float


def profile(path: str) -> list[BandProfile]:
    ds = gdal.Open(path)
    out = []
    for i in range(1, ds.RasterCount + 1):
        band = ds.GetRasterBand(i)
        # approx_ok=False so the stats reflect real pixels, not a decimated read;
        # a decimated read can miss a nodata border entirely on a large scene.
        mn, mx, mean, _std = band.ComputeStatistics(False)
        # GetMaskBand gives valid-pixel coverage whether nodata is declared
        # by value or by an alpha/mask band, which the raw value test misses.
        mask = band.GetMaskBand().ComputeBandStats()[0] / 255.0
        out.append(BandProfile(i, band.GetNoDataValue(), mask, mn, mx, mean))
    return out


def assert_warp_integrity(src: str, dst: str, layer: str,
                          nodata_tolerance: float = 0.05) -> None:
    before, after = profile(src), profile(dst)

    # 1. Band count is a hard contract — never "close enough".
    if len(before) != len(after):
        warp_rejected.add(1, {"layer": layer, "reason": "band_count"})
        raise ValueError(f"{layer}: band count {len(before)} → {len(after)}")

    for b0, b1 in zip(before, after):
        # 2. The nodata DECLARATION must survive, or the sentinel becomes data.
        if (b0.nodata is None) != (b1.nodata is None):
            warp_rejected.add(1, {"layer": layer, "reason": "nodata_declaration"})
            raise ValueError(f"{layer} band {b1.index}: nodata declaration lost")

        # 3. Valid coverage must not collapse — this catches grid clipping.
        lost = b0.valid_ratio - b1.valid_ratio
        if lost > nodata_tolerance:
            warp_rejected.add(1, {"layer": layer, "reason": "nodata_expansion"})
            raise ValueError(
                f"{layer} band {b1.index}: valid coverage {b0.valid_ratio:.3f} "
                f"→ {b1.valid_ratio:.3f}")

        # 4. Band identity by statistical fingerprint — a reorder changes which
        #    band carries which distribution even though the count is intact.
        if b0.mean > 0 and abs(b1.mean - b0.mean) / b0.mean > 0.15:
            warp_rejected.add(1, {"layer": layer, "reason": "band_identity"})
            raise ValueError(
                f"{layer} band {b1.index}: mean {b0.mean:.1f} → {b1.mean:.1f} "
                "— band order or scaling changed")

Two details carry most of the value. Using the mask band rather than testing pixels against the nodata value covers coverages that express validity through an alpha or mask band instead of a sentinel — a distinction that silently defeats value-based checks on exactly the datasets most likely to use masks. And computing exact statistics rather than approximate ones matters because a decimated read samples the interior of a scene and can miss a nodata border entirely, which is the defect you are hunting.

The band-identity check via mean comparison is deliberately loose. Resampling shifts means slightly; a fifteen percent tolerance passes ordinary resampling and fails a swap between spectrally distinct bands. For datasets where two bands have genuinely similar distributions, add a correlation check against the source band instead, or carry band descriptions in metadata and assert on those.

Assertion order from cheapest and most decisive to most tolerant Four assertions are arranged in a descending ladder. Band count is checked first as a hard contract with no tolerance. The nodata declaration is checked second, also with no tolerance, because losing it converts a sentinel into data. Valid-coverage loss is checked third with a five percent tolerance, catching grid clipping. Band identity by statistical fingerprint is checked last with a fifteen percent tolerance, because resampling legitimately moves band statistics a little. A note explains that ordering the checks this way makes the rejection reason specific rather than generic. Order the assertions so the rejection reason is specific 1. band count hard contract · tolerance 0 · reason=band_count a changed count is never acceptable 2. nodata declaration present hard contract · tolerance 0 · reason=nodata_declaration losing it turns a sentinel into data 3. valid-coverage loss tolerance 5% · reason=nodata_expansion catches target-grid clipping 4. band identity fingerprint tolerance 15% · reason=band_identity resampling moves means a little; a swap moves them a lot Valid-coverage ratio across a month of warps, with one rejected scene Daily valid-coverage ratios are plotted for a coverage pipeline. Most days sit close to ninety-six percent with small variation caused by cloud. One day drops to twenty-two percent, far outside the variation, and is marked as a rejected scene caused by a target-grid offset. A tolerance band five percentage points below the rolling mean is drawn, and the note states that cloud variation stays inside it while a grid fault does not. Cloud variation stays inside the band; a grid offset does not tolerance band rejected — target grid offset 96% 22% The rejected scene never entered the pyramid, so completeness shows the corresponding tiles as missing.

Verification: manufacture each defect and confirm rejection

The assertions are only trustworthy if you have seen them fire. Each defect is easy to manufacture from a known-good scene.

Produce nodata expansion by warping into a target extent deliberately offset from the source footprint; the valid-coverage assertion must reject it, with reason="nodata_expansion" on the counter. Produce band reordering by warping with an explicit band list in reversed order; the identity assertion must reject it. Produce a dropped declaration by warping with the nodata value unset in the output creation options; the declaration assertion must reject it and, importantly, must do so before the coverage assertion, since a missing declaration makes the coverage figure meaningless.

Then run the assertions against a normal daily warp and confirm they pass with margin. A tolerance that only just passes ordinary operation will reject a legitimately cloudy scene next week.

Gotchas

Comparing against a fixed nodata threshold. “Reject if nodata exceeds twenty percent” fails on a scene that is legitimately mostly ocean and passes on a scene that lost half its data but started at five percent. Always compare against the source’s own ratio.

Approximate statistics. Fast statistics decimate the read and can entirely miss a nodata border, which is the shape most warp clipping produces. Use exact computation for the assertion even if approximate values are fine for dashboards.

Assuming nodata is a value. Coverages using alpha or mask bands have no nodata value at all, and a value-based check reports zero nodata on a scene that is three-quarters transparent. Read validity through the mask band.

Running the assertions after publication. These checks are a gate, not a monitor. Running them on the published coverage tells you about a problem your users already have, and the whole point of asserting against the source is that the source is still available at warp time — the trust-boundary discipline described in defining spatial data trust boundaries.

FAQ

What tolerance should valid-coverage loss use?

Start at five percent and calibrate against a month of normal runs. Resampling at a different resolution genuinely changes edge coverage slightly, so a zero tolerance flaps. If your warps are grid-aligned and same-resolution, tighten to one percent — a tighter bound on a stable pipeline is strictly more useful.

Can I use band descriptions instead of statistics for identity?

Yes, and prefer it where the format and toolchain preserve descriptions reliably. Asserting that band three is still described as nir is far more direct than inferring identity from a mean. The statistical fingerprint is the fallback for formats and pipelines that drop descriptions, which is many of them.

Should a rejected warp halt the whole layer?

Halt that scene, not the layer. A single bad source scene should be quarantined and reported while other scenes continue, since coverage pipelines routinely process hundreds of scenes per run. If the rejection rate across scenes crosses a threshold, that is a layer-level incident and belongs on the pager.

How do these metrics feed the pyramid checks?

A scene rejected at warp never enters the pyramid, so completeness will show the corresponding tiles as missing. That is the correct outcome and the two signals should be read together: a completeness gap coinciding with a spike in gis.raster.warp_rejected_total is a diagnosed incident rather than a mystery, which is why both are described in raster and tile pipeline observability.

Does this apply to vector tiles?

Only partly. Vector tiles have no bands or nodata, but the analogous assertion — that the feature count and attribute schema of the tile match the source layer — serves the same purpose, and the schema half of it is the subject of schema and attribute drift detection.