Article

Asserting Raster Outputs Without Golden Files

Compare the properties a user depends on — CRS, transform, shape, band statistics and nodata — rather than the file’s bytes. A four-assertion helper covers everything a golden GeoTIFF would, survives GDAL upgrades that change metadata tags or compression defaults, and reports which property differed instead of files are not identical. It builds on the Testing Geospatial CLI Tools guide, part of the broader CLI Architecture & Design Patterns reference.

Prerequisites

  • Python 3.10 or later
  • pip install pytest rasterio numpy pyproj
  • An existing test that writes a raster to tmp_path; the helper below replaces whatever assertion follows it

Why a Golden File Fails for Reasons You Do Not Care About

A GeoTIFF carries far more than pixels. It carries the library version that wrote it, the compression predictor, block layout, the order metadata tags were written in, and — for some drivers — a timestamp. None of those are things your tool decided, and all of them can change when a dependency is upgraded.

What in a GeoTIFF your tool actually decided Your tool decides the CRS, the geotransform, the shape, the band values and the nodata value. GDAL decides the writer version tag, the compression predictor, the block layout when unspecified, the tag ordering and any embedded overview structure. A byte comparison tests both halves. your tool decided these the CRS the geotransform and shape the band values the nodata value and mask GDAL decided these the writer version tag compression predictor defaults block layout when unspecified tag ordering in the header

A checksum over the file tests both columns with equal weight. When the suite goes red after a dependency bump, the failure carries no information about which column moved.

Complete Working Implementation

# tests/raster_assertions.py
"""Property assertions for raster outputs, in place of golden-file comparison."""
from __future__ import annotations

from dataclasses import dataclass
from pathlib import Path

import numpy as np
import rasterio
from pyproj import CRS


def _coordinate_tolerance(crs) -> float:
    """Absolute tolerance in the CRS's own units."""
    unit = (CRS(crs).axis_info[0].unit_name or "").lower()
    if unit.startswith("degree"):
        return 1e-8
    if "foot" in unit:
        return 1e-2
    return 1e-3


@dataclass(frozen=True)
class RasterExpectation:
    epsg: int
    shape: tuple[int, int] | None = None          # (height, width)
    bounds: tuple[float, float, float, float] | None = None
    band_count: int = 1
    dtype: str | None = None
    nodata: float | None = None
    value_range: tuple[float, float] | None = None   # (min, max) of valid data
    min_valid_fraction: float = 0.0                  # guards against an all-nodata result


def assert_raster(path: Path, expected: RasterExpectation) -> None:
    with rasterio.open(path) as src:
        # 1. georeferencing
        assert src.crs is not None, "output has no CRS"
        actual_epsg = CRS(src.crs).to_epsg()
        assert actual_epsg == expected.epsg, (
            f"CRS: got EPSG:{actual_epsg}, expected EPSG:{expected.epsg}"
        )

        # 2. geometry
        assert src.width > 0 and src.height > 0, "output is degenerate"
        if expected.shape is not None:
            assert (src.height, src.width) == expected.shape, (
                f"shape: got {(src.height, src.width)}, expected {expected.shape}"
            )
        if expected.bounds is not None:
            tol = _coordinate_tolerance(src.crs)
            for got, want, name in zip(src.bounds, expected.bounds,
                                       ("left", "bottom", "right", "top")):
                assert abs(got - want) < tol, (
                    f"bounds.{name}: got {got!r}, expected {want!r} (tol {tol})"
                )

        # 3. structure
        assert src.count == expected.band_count, (
            f"band count: got {src.count}, expected {expected.band_count}"
        )
        if expected.dtype is not None:
            assert src.dtypes[0] == expected.dtype, (
                f"dtype: got {src.dtypes[0]}, expected {expected.dtype}"
            )
        if expected.nodata is not None:
            assert src.nodata == expected.nodata, (
                f"nodata: got {src.nodata!r}, expected {expected.nodata!r}"
            )

        # 4. data
        band = src.read(1, masked=True)
        valid = band.compressed()
        fraction = valid.size / band.size if band.size else 0.0
        assert fraction >= expected.min_valid_fraction, (
            f"only {fraction:.1%} of pixels are valid; "
            f"expected at least {expected.min_valid_fraction:.1%}"
        )
        assert np.isfinite(valid).all(), "output contains NaN or infinity"
        if expected.value_range is not None and valid.size:
            lo, hi = expected.value_range
            assert valid.min() >= lo and valid.max() <= hi, (
                f"values span [{valid.min()}, {valid.max()}], "
                f"expected within [{lo}, {hi}]"
            )

Using it reads as a description of the contract rather than a comparison:

# tests/test_warp.py
from tests.raster_assertions import RasterExpectation, assert_raster


def test_warp_to_web_mercator(raster_path, tmp_path, runner, app):
    out = tmp_path / "out.tif"
    result = runner.invoke(app, ["raster", "warp", str(raster_path), str(out),
                                 "--dst-crs", "EPSG:3857"])
    assert result.exception is None
    assert result.exit_code == 0

    assert_raster(out, RasterExpectation(
        epsg=3857,
        band_count=1,
        dtype="float32",
        min_valid_fraction=0.90,     # warping leaves some edge nodata
        value_range=(0.0, 65535.0),
    ))

Each property answers a different question about the output, and they fail independently. Laid out against the bugs they catch, it is clear why no single one of them is sufficient.

Which assertion catches which bug The CRS check catches a wrong target CRS. The valid-fraction check catches a transform error that leaves the output mostly nodata. The value-range check catches a scaling error. The dtype check catches a narrowing that clips values. No single assertion catches more than one. wrong CRS bad transform scaling bug dtype clip crs / to_epsg bounds valid fraction value range dtype filled = catches it

The empty column is the useful thing to look for when adding an assertion. If a bug you have actually shipped does not light up any row, that is the assertion you are missing.

Step Annotations

  1. result.exception is None comes before the exit-code assertion. A runner that caught an exception reports exit code 1, and asserting on the code first hides the traceback that would explain it.
  2. to_epsg() rather than string comparison. "EPSG:3857", "epsg:3857" and a full WKT string all describe the same CRS. Normalising through pyproj compares meaning rather than spelling.
  3. min_valid_fraction is the assertion that catches silent disaster. A warp with the wrong transform frequently produces a correctly-shaped raster that is almost entirely nodata. Shape and CRS both pass; this does not.
  4. value_range is deliberately loose. Its job is to catch a scaling bug — values that came out in the thousands when they should be fractions — not to pin exact numbers. A tight range here reintroduces the fragility golden files had.
  5. The tolerance is derived from the CRS, not passed in. A test written for a projected CRS and later re-parametrised over a geographic one keeps working, because the tolerance follows the units.

One Named Gotcha: src.nodata Is None, Not NaN, When Unset

Asserting src.nodata == expected.nodata looks symmetric until the expectation is float("nan"), because nan == nan is False and the assertion fails against an output that is correct. Worse, a driver that does not carry the tag returns None, and None == nan is False too — so the same assertion fails for two unrelated reasons with one message.

Handle the two cases explicitly:

if expected.nodata is not None:
    if np.isnan(expected.nodata):
        assert src.nodata is not None and np.isnan(src.nodata), (
            f"nodata: got {src.nodata!r}, expected NaN"
        )
    else:
        assert src.nodata == expected.nodata

If your tool uses NaN as its nodata sentinel, prefer a distinct value such as -9999.0 in formats that support it. NaN propagates through arithmetic in ways that make downstream statistics quietly wrong, and it costs you this assertion.

A byte comparison and a property comparison also fail differently, and the difference is most of the maintenance cost over a suite’s life.

What each style of failure tells you A golden-file comparison reports that two files are not identical, with no indication of which property changed. A property assertion names the property, the value produced and the value expected, which is enough to act on without opening either file. golden file AssertionError: files differ next step: open both in QGIS property assertion AssertionError: CRS: got EPSG:27700, expected 3857 next step: read one line of code The right-hand message is also the one that stays useful in a CI log nobody can attach a debugger to.

That is the practical argument for the helper carrying explicit messages on every assertion rather than relying on pytest’s expression rewriting: the rewritten output is excellent for scalars and unreadable for a BoundingBox of four floats.

Verification

Confirm the assertions have teeth by breaking one thing at a time and checking that exactly one message appears:

# Deliberately pass the wrong target CRS in the test and confirm the failure
# names the CRS rather than reporting a generic mismatch.
pytest tests/test_warp.py -q 2>&1 | grep -E "^E .*(CRS|bounds|valid|dtype)"

The expected output names the property:

E       AssertionError: CRS: got EPSG:27700, expected EPSG:3857

If a deliberate break produces a failure that does not name the property, the helper is missing a message — which matters more than it sounds, because that message is what a future maintainer reads at 5pm on a Friday.

Extending the Same Idea to Vector Output

The vector case is structurally identical and differs only in which properties matter. A written GeoPackage or FlatGeobuf carries a CRS, a feature count, a geometry type, a field schema and an extent, and every one of those is something your tool decided. What it also carries — the SQLite page size, the order features happen to be stored in, the metadata tables the driver maintains — is GDAL’s business, and comparing it produces the same false failures a raster golden file does.

def assert_vector(path, epsg, feature_count, geometry_type, fields):
    import pyogrio

    info = pyogrio.read_info(path)
    assert info["crs"].endswith(str(epsg)), f"CRS: {info['crs']}"
    assert info["features"] == feature_count, (
        f"feature count: got {info['features']}, expected {feature_count}"
    )
    assert info["geometry_type"] == geometry_type
    assert set(info["fields"]) == set(fields), (
        f"schema drift: {set(info['fields']) ^ set(fields)}"
    )

pyogrio.read_info is worth knowing about here specifically because it reads the header without materialising any features, so the assertion is fast even against a large output. That matters when the same helper is used in a test that writes a hundred thousand features to check a chunking path — the write is the expensive part, and the assertion should not double it.

The one property with no raster equivalent is field-name stability. Attribute names are a public interface for anything joining against the output, and they are easy to change accidentally through a rename in an intermediate dataframe. Asserting on the exact set — with a symmetric difference in the message, so a failure names both the missing and the unexpected fields — catches the whole class in one line.

Geometry itself is best left out of the assertion unless the test is specifically about geometry. Comparing shapes requires a tolerance argument, invites arguments about vertex order, and rarely catches a bug that the feature count and extent do not. Where a test genuinely is about geometry — a clip, a simplification, a buffer — assert on a derived scalar instead: total area, total length, or the count of parts. Those are stable under vertex reordering and read clearly in a failure.