Asserting Raster Band Statistics in pytest

9 min read

A band statistic — a mean, a minimum, a valid-cell count — is the cheapest summary a raster offers and the easiest to assert wrongly. This guide sits beneath raster and grid assertions and shows how to write statistics assertions with rasterio 1.3+ and NumPy that actually discriminate: masking the fill before computing, asserting the population alongside the value, choosing the bound form from the band’s semantics rather than by habit, and casting before any arithmetic so a difference cannot wrap into a positive number.

The failure this prevents is not a wrong number that looks wrong. It is a plausible number computed over the wrong population, which passes review, passes CI, and is discovered months later when a total does not reconcile.

Root cause: the statistic and the population are two claims

band.mean() answers a question, and the question includes an unstated clause: over which cells? If nodata is unmasked, the answer covers fill as well as measurements. If a window was read, it covers that window. If a mask was applied upstream, it covers whatever survived. Every one of those produces a number, and none of them raises.

The consequence is that asserting a value alone is only half an assertion. The other half is the count of cells that contributed — the population — and it is the half that discriminates. A mean of 11.9 over 8,412 valid cells and a mean of 11.9 over 10,000 cells including 1,588 zeros are different results from a suite’s point of view, and only the count tells them apart.

One array, three populations, three means A single band is read three ways. Reading it raw includes every cell, so the declared fill value contributes to the mean and the contributing count equals the full cell total. Reading it with masking enabled excludes the fill cells, so the mean covers measurements only and the contributing count is smaller. Reading a window covers only a sub-region, producing a third answer that is equally valid and equally different. All three calls return a number and none of them raises an error, so the only way a test can tell which population it measured is to assert the contributing cell count alongside the value. one band 10 000 cells, 1 588 fill read raw — fill included mean 11.9 · contributing cells 10 000 read masked — fill excluded mean 14.2 · contributing cells 8 412 read a window mean 13.4 · contributing cells 2 048 all three return a number none of them raises Asserting only the value cannot distinguish these three. Asserting the contributing count alongside it distinguishes all of them, and costs one line.

Parameter reference

Element Choice that discriminates Common mistake
Read call src.read(1, masked=True) src.read(1) — fill included
Population Assert valid_count explicitly Assert only the statistic
Bound form Relative epsilon for magnitudes Absolute for a quantity spanning orders
Bound form Absolute for bounded quantities Relative near zero
Arithmetic dtype Cast to float64 before differences uint8 subtraction wraps
Range check Assert min and max against physical limits Trust the mean alone
Failure message Statistic, value, expected, population “statistics differ”

Step-by-step implementation

The gate targets rasterio 1.3+, NumPy and pytest 7+.

Step 1 — Read masked, always

Masking is a read-time decision, and getting it wrong contaminates everything downstream. Make it the only way the fixture is opened.

import rasterio
import numpy as np

def read_band(path: str, band: int = 1) -> np.ma.MaskedArray:
    """Read one band with the declared nodata masked out."""
    with rasterio.open(path) as src:
        if src.nodata is None:
            raise AssertionError(f"{path}: no nodata declared — masking is not possible")
        return src.read(band, masked=True)

Raising when nodata is undeclared is deliberate. A raster with no declared fill either genuinely has none — in which case the contract should say so — or has fill that nothing will mask, which is the defect this whole guide exists to prevent.

Step 2 — Compute the statistic and its population together

Returning both from one function makes it impossible to assert one without the other.

from dataclasses import dataclass

@dataclass(frozen=True)
class BandStats:
    mean: float
    minimum: float
    maximum: float
    valid_cells: int
    total_cells: int

def band_stats(arr: np.ma.MaskedArray) -> BandStats:
    valid = int(arr.count())            # unmasked cells only
    return BandStats(
        mean=float(arr.mean()),
        minimum=float(arr.min()),
        maximum=float(arr.max()),
        valid_cells=valid,
        total_cells=int(arr.size),
    )

Step 3 — Assert the population first

The order matters for the same reason it matters elsewhere: a wrong population makes every value assertion meaningless, so failing on it first names the cause rather than a consequence.

import pytest

EXPECTED_VALID = 8412
EXPECTED_MEAN = 14.2
REL_EPS = 0.01                     # 1% — the band is a magnitude, not a bounded ratio

def test_band_statistics(dem_path):
    stats = band_stats(read_band(dem_path))

    assert stats.valid_cells == EXPECTED_VALID, (
        f"population changed: {stats.valid_cells} valid of {stats.total_cells} "
        f"(expected {EXPECTED_VALID}) — the mean below is not comparable"
    )
    assert abs(stats.mean - EXPECTED_MEAN) / EXPECTED_MEAN <= REL_EPS, (
        f"mean {stats.mean:.3f} outside {REL_EPS:.1%} of {EXPECTED_MEAN} "
        f"over {stats.valid_cells} cells"
    )

Step 4 — Bound the range as well as the centre

A mean can be correct while individual cells are impossible. For any physically bounded quantity, the range assertion is the one that catches an overshoot from interpolation or a wrapped dtype.

def test_band_values_are_physically_possible(reflectance_path):
    stats = band_stats(read_band(reflectance_path))
    assert 0.0 <= stats.minimum, f"negative reflectance: {stats.minimum}"
    assert stats.maximum <= 1.0, f"reflectance above unity: {stats.maximum}"
Bound form follows the band's semantics Three band semantics with the bound form each requires. A magnitude such as elevation or population density spans several orders of magnitude, so a relative epsilon is correct because it scales with the value, whereas an absolute bound is simultaneously too tight at the low end of the range and too loose at the high end. A bounded ratio such as reflectance or a percentage occupies a fixed interval, so an absolute bound applies uniformly while a relative bound becomes meaningless as the value approaches zero. A count is discrete and usually admits an exact assertion. A closing note records that choosing the bound form by habit rather than by the band's semantics produces assertions that fail for correct data at one end of the range. BAND SEMANTICS USE BECAUSE Magnitude elevation, density, count per cell relative epsilon spans orders of magnitude — an absolute bound is too tight low and too loose high Bounded ratio reflectance, percentage, index absolute bound fixed interval — a relative bound becomes meaningless as the value approaches zero Discrete class land cover, category code exact set equality any interpolation is itself the defect — an intermediate value must not exist Choosing by habit rather than by semantics produces bounds that fail for correct data at one end of the range and pass for defects at the other.

Verify the fix

Provoke the population failure once to confirm the assertion order works:

pytest -q tests/test_band_stats.py -v

Temporarily read without masked=True and confirm the count assertion fails first, with a message naming both counts. If instead the mean assertion fails first, the ordering is wrong and a future reader will be sent to investigate a value that was never the problem.

The dtype trap in difference assertions

Comparing two bands is the most common raster assertion after a single statistic, and it has a failure mode that produces no error and a confidently wrong answer. Unsigned integer arithmetic wraps: subtracting 10 from 4 in uint8 yields 250, not −6. A difference raster computed that way has a plausible dtype, a plausible range, and values that are the arithmetic complement of the truth in exactly the cells where the second band exceeds the first.

uint8 subtraction wraps; casting first does not Two computations of the same difference. In the first, two unsigned eight-bit bands are subtracted directly; where the second value of ten exceeds the first value of four the result wraps to two hundred and fifty instead of negative six, and the resulting difference raster is plausible in dtype and range while being wrong in every cell where the second band was larger. In the second, both bands are cast to a signed or floating type before the subtraction, so the same cells yield negative six as intended. A closing note records that a mean computed over the wrapped raster is also plausible, so the defect is caught only by a range assertion or by the cast itself. Direct subtraction on uint8 band A cell = 4 band B cell = 10 A - B # uint8 wraps below zero result = 250 dtype plausible range plausible Cast first band A cell = 4 band B cell = 10 A.astype("f8") - B signed arithmetic result = −6 as intended A mean over the wrapped raster is also plausible. Only a range assertion — or performing the cast in the first place — catches this.

Choosing which statistics are worth asserting

Not every summary a band offers is worth a test. Three are, one is conditional, and two are usually noise.

The valid-cell count is the most valuable single number in a raster suite, because it is sensitive to almost every structural defect — a changed mask, a truncated read, a shifted window, an extra fill value — and insensitive to the ordinary variation that makes value assertions brittle. If a suite asserts one thing about a band, this is the one.

The mean is the workhorse: sensitive enough to catch a scaling error, a unit change or a wrong band, and stable enough not to move when a handful of cells do. Bound it relatively for magnitudes and absolutely for ratios.

The min and max are the range assertion, and they catch what the mean cannot: interpolation overshoot, wrapped arithmetic, and impossible values in a physically bounded quantity. They are cheap and they fail loudly.

The standard deviation is conditional. It catches a smoothing or sharpening that leaves the mean untouched, which is genuinely useful after any resampling — and it is unstable on small fixtures, so it belongs on production-shaped data in a scheduled run rather than in a fast lane.

The median and the mode are usually noise in a continuous band, because both move discontinuously with small changes in the data and produce failures nobody can attribute. The exception is a categorical band, where the mode is the dominant class and asserting it is exactly right.

The general principle is to prefer statistics that move for the reasons you care about and stay still otherwise. A suite full of assertions that fail for legitimate data variation trains a team to widen thresholds, and once that habit exists the sensitivity that would have caught a real defect has already been spent.

Failure modes and edge cases

  1. A fully-masked band. If every cell is nodata, mean() on a masked array returns a masked value rather than a number, and comparisons against it are neither true nor false. Assert valid_cells > 0 before any statistic.
  2. Float nodata compared by equality. A nodata value of -9999.0 stored as float32 may not compare exactly equal to the declared float64 constant, leaving a few fill cells unmasked. Compare with a tolerance, or declare the value in the band’s own dtype.
  3. NaN in the data. NaN is not masked by a nodata declaration and propagates through every statistic, turning a mean into NaN — which compares False against any bound, so an assertion using <= silently passes. Check for finiteness explicitly.
  4. Relative bounds near zero. A relative epsilon on a quantity whose expected value is close to zero produces an absurdly tight bound. Switch to an absolute bound whenever the expected value can approach zero.
  5. Statistics cached in the file. Some drivers store precomputed statistics in metadata, and those can be stale relative to the pixel data. Compute from the array rather than reading the cached values, or assert that the two agree.
  6. Comparing across different windows. Two statistics computed over different windows are not comparable however close their values look. Recording the window alongside the statistic makes that visible in the report rather than in someone’s assumption.

Conclusion

A band statistic assertion is two claims, and most suites make only one of them. Reading masked, asserting the contributing cell count before the value, choosing the bound form from the band’s semantics, and casting before any difference converts a check that passes over plausible-looking errors into one that names the population it measured — the statistics layer of raster and grid assertions.