Comparing Resampled Rasters with Tolerance

9 min read

Once a raster has been resampled, a cell-by-cell comparison against its source is the wrong instrument. Resampling redistributes values deliberately — that is its entire purpose — so an exact comparison reports differences in almost every cell and none of them are defects. This guide sits beneath raster and grid assertions and shows what to assert instead: the quantities each resampling method actually preserves, why a percentile bound discriminates where a maximum does not, and the two failure classes — cubic overshoot and accidental interpolation of categorical data — that statistical assertions miss entirely.

Root cause: different methods conserve different things

A resampling method is a choice about what to preserve when cells do not line up, and each choice preserves something different. Asserting the wrong invariant produces a test that fails for entirely correct behaviour, which is how resampling assertions acquire their reputation for noise.

Nearest neighbour copies the value of the closest source cell. It preserves the set of values exactly — no value appears in the output that was not in the input — and changes only how many cells hold each one. It is the only method valid for categorical data, and asserting the value set is both a correctness check and an enforcement of that requirement.

Bilinear and cubic compute a weighted combination of neighbouring cells, so intermediate values appear by construction. Neither preserves the value set, both preserve the mean approximately, and cubic can produce values outside the input range entirely.

Average assigns each output cell the mean of the source cells it covers. It preserves the mean well and the total only when the cell areas work out, which is the distinction that decides whether an extensive or an intensive assertion is correct.

What each resampling method preserves Four rows pairing a resampling method with the property it preserves and the assertion that property implies. Nearest neighbour preserves the exact set of distinct values and alters only their counts, so the correct assertion is set equality, and this is the only method valid for categorical data. Bilinear preserves the mean approximately while introducing intermediate values by construction, so the correct assertion is a relative bound on the mean. Cubic preserves the mean approximately and can overshoot beyond the range of its inputs at sharp edges, so it needs a range bound in addition to the mean. Average preserves the mean well and preserves the total only when the cell areas divide evenly, so the choice between an extensive and an intensive assertion depends on the quantity being resampled. METHOD PRESERVES SO ASSERT Nearest neighbour the only categorical option the exact set of values; only the counts change set equality — exact Bilinear weighted neighbours the mean, approximately; intermediates appear mean within a relative epsilon Cubic wider kernel the mean; can overshoot beyond the input range mean + an explicit range bound Average area-weighted mean the mean well; the total only when areas divide evenly mean, or total if extensive

Parameter reference

Assertion Form Right when
Value set unchanged Exact set equality Nearest neighbour, or categorical data
Mean preserved ∣μ′−μ∣/μ≤ε\lvert \mu' - \mu \rvert / \mu \le \varepsilon Any interpolating method, intensive quantity
Total preserved Relative bound on the sum Extensive quantity, area-weighted method
Range not exceeded min⁡′≥min⁡\min' \ge \min, max⁡′≤max⁡\max' \le \max Cubic, or any physically bounded band
Per-cell difference 95th or 99th percentile bound A cell-level guarantee is genuinely needed
Valid-cell count Exact, derived arithmetically Always — the mask changes shape too
Extent preserved Within a fraction of a cell Always

Step-by-step implementation

The assertions target rasterio 1.3+ and NumPy.

Step 1 — Assert the geometry before the values

A resampled raster has a new shape and a new geotransform, and both should follow arithmetically from the requested scale. Getting those wrong makes every value assertion meaningless.

import numpy as np

def assert_resampled_geometry(src_profile, dst_profile, scale: float):
    expected_w = round(src_profile["width"] * scale)
    expected_h = round(src_profile["height"] * scale)
    assert (dst_profile["width"], dst_profile["height"]) == (expected_w, expected_h), (
        f"resampled shape {dst_profile['width']}x{dst_profile['height']} "
        f"does not follow from scale {scale} (expected {expected_w}x{expected_h})"
    )
    # Pixel size must scale inversely; the origin must not move.
    assert np.isclose(dst_profile["transform"].a, src_profile["transform"].a / scale)
    assert np.isclose(dst_profile["transform"].c, src_profile["transform"].c)

Step 2 — Assert the invariant the method actually preserves

One helper per method keeps the choice explicit at the call site, which is where it should be visible.

def assert_value_set_preserved(before: np.ma.MaskedArray, after: np.ma.MaskedArray):
    """Nearest neighbour and mode only — no new value may appear."""
    new_values = set(np.unique(after.compressed())) - set(np.unique(before.compressed()))
    assert not new_values, (
        f"resampling introduced values absent from the source: {sorted(new_values)[:5]} — "
        "an interpolating method was used on categorical data"
    )

def assert_mean_preserved(before, after, rel_eps: float = 0.01):
    mu, mu2 = float(before.mean()), float(after.mean())
    assert abs(mu2 - mu) / abs(mu) <= rel_eps, (
        f"mean moved from {mu:.4f} to {mu2:.4f}, beyond {rel_eps:.1%}"
    )

Step 3 — Bound a percentile, not the maximum

Where a cell-level guarantee is genuinely required, the maximum difference after resampling is dominated by edge cells and cells adjacent to nodata, where interpolation has least information. Bounding it forces a threshold so loose it stops discriminating anywhere else.

def assert_cell_difference(before, after_upsampled, p: float = 95, bound: float = 0.5):
    """Compare only where both are valid; bound a percentile of |difference|."""
    both_valid = ~before.mask & ~after_upsampled.mask
    diff = np.abs(before[both_valid].astype("f8") - after_upsampled[both_valid].astype("f8"))
    measured = float(np.percentile(diff, p))
    assert measured <= bound, (
        f"p{p:.0f} cell difference {measured:.4f} exceeds {bound} "
        f"(max was {float(diff.max()):.4f}, dominated by edges)"
    )

Step 4 — Catch cubic overshoot explicitly

Cubic interpolation can produce values outside the range of its inputs. For a physically bounded quantity that is an impossible value, and no statistical assertion notices.

def assert_no_overshoot(before, after, lo=None, hi=None):
    lo = float(before.min()) if lo is None else lo
    hi = float(before.max()) if hi is None else hi
    assert float(after.min()) >= lo - 1e-9, f"undershoot: {after.min()} below {lo}"
    assert float(after.max()) <= hi + 1e-9, f"overshoot: {after.max()} above {hi}"
Cubic overshoot at a step edge A profile across a sharp edge in a source raster, resampled three ways. The source steps abruptly from a low value to a high value. The nearest-neighbour result reproduces the step exactly with no intermediate values. The bilinear result forms a straight ramp between the two levels and stays entirely within the input range. The cubic result forms a smooth curve that dips below the low value immediately before the step and rises above the high value immediately after it, leaving the input range in both directions. Dashed lines mark the input minimum and maximum, and the excursions beyond them are annotated as impossible values for any physically bounded quantity, invisible to a mean or percentile assertion. input max input min nearest — exact step bilinear — inside the range undershoot below the input min overshoot above the input max For a bounded quantity — a reflectance, a percentage, a probability — both excursions are impossible values. A mean assertion passes; a percentile assertion passes; only an explicit range bound catches it.

Verify the fix

Run the comparison against a deliberately wrong method to confirm each assertion discriminates:

pytest -q tests/test_resample.py -v

Resample a categorical fixture with bilinear and confirm the value-set assertion fails. Resample a step-edge fixture with cubic and confirm the range assertion fails while the mean assertion passes — that pairing is the clearest demonstration that the two checks cover different classes.

The mask changes shape too

Resampling does not only move values; it moves the boundary between valid and fill. An output cell that partially overlaps a masked region has to be either valid or not, and the rule differs by method — nearest neighbour takes the nearest source cell’s state, while an averaging method may treat any partial coverage as valid, as invalid, or as valid with a reduced weight.

The consequence is that the valid-cell count after resampling is not the source count scaled by the area ratio, and asserting that it is produces failures for correct behaviour. What can be asserted is that the count follows a stated rule and that the mask has no holes it did not have before — a fill region that becomes speckled after resampling indicates the mask was interpolated as if it were data, which is a real and common defect.

A mask edge, resampled correctly and incorrectly Three panels showing a mask boundary. In the source the fill region has a clean straight edge against the valid data. After nearest-neighbour resampling the edge remains clean and has simply moved to the nearest cell boundary in the new grid. After resampling with the mask interpolated as though it were data, the edge becomes a speckled band of partially-valid cells scattered along what should be a clean boundary, which is the recognisable signature of the defect. A caption states that the correct assertion is the absence of new holes and a connected edge, rather than a valid-cell count that scales with area. Source fill valid a clean edge Nearest neighbour fill valid edge moved to the nearest cell — still clean Mask interpolated a speckled band — the defect signature Assert that the mask gained no holes and its edge stays connected. Do not assert that the valid-cell count scales with area — it does not. The speckle is what an interpolated mask looks like, and it is invisible in any value statistic.

Recording which method ran

The single most useful thing a resampling test can record is the method that produced the output. Every assertion on this page depends on it — the same output is correct under one method and defective under another — and the method is exactly the parameter most likely to change without anyone intending it, because it is often a default rather than a choice.

Two places are worth writing it. In the run summary, alongside the scale factor and the resulting dimensions, so a later question about why a comparison behaved differently is answerable without re-running anything. And in the output’s own metadata, as a tag on the written file, so a downstream consumer inheriting the raster can tell whether the values it holds were interpolated or copied. A categorical raster that has been through a bilinear resampling is not recoverable, and the only defence is that somebody notices before it is used.

The assertion that follows is short and worth having: read the method back from the artefact and compare it against the one the pipeline configuration requested. A mismatch means a default was applied somewhere between the request and the write, which is a defect in the pipeline rather than in the data — and one that no examination of the pixels would ever reveal.

Failure modes and edge cases

  1. Comparing an upsampled output against its source cell by cell. Upsampling creates cells that had no source, so a direct comparison is undefined for most of them. Downsample the output back, or compare aggregates.
  2. A relative bound on a mean near zero. A difference raster has a mean close to zero by construction, so a relative epsilon becomes an absurdly tight bound. Use an absolute one whenever the expected value can approach zero.
  3. Interpolating categorical data. Bilinear on a land-cover raster produces class codes that do not exist, and the resulting map is meaningless in a way that looks smooth. The value-set assertion is the enforcement mechanism, not a nicety.
  4. Ignoring the nodata contribution. Interpolating a cell adjacent to fill mixes the sentinel into the result unless the method is mask-aware, producing a halo of wrong values along every mask edge. Resample with masking enabled and assert the halo is absent.
  5. Assuming the total is conserved. Only an area-weighted method conserves an extensive quantity, and only approximately. Asserting a conserved total after bilinear resampling fails for entirely correct behaviour.
  6. Bounding the maximum difference. Edge cells dominate it, so the threshold has to be loose enough to admit them, at which point it admits everything. Bound a percentile and record the maximum separately.

Conclusion

Comparing resampled rasters is an exercise in asserting the right invariant. The method determines what survives — a value set, a mean, a total — and asserting anything else produces failures for correct behaviour, which is how these checks come to be disabled. Pairing the invariant with an explicit range bound for interpolating methods, a percentile rather than a maximum for cell-level guarantees, and a mask-shape assertion for the boundary covers the resampling half of raster and grid assertions.