Validating Raster Predictions and Accuracy Assessment

Accuracy assessment is the part of a raster modelling project that other people will actually rely on. A number reported without a sampling design behind it is worthless; a number reported without a confidence interval is worse, because it implies precision it does not have. This topic — part of Raster Machine Learning & Model Inference — covers the mechanics of turning a prediction raster and some reference data into figures that hold up.

Two ideas do most of the work. The first is that validation samples must be spatially independent of training samples, which is a property of the design, not of the metric. The second is that map area and estimated area are different quantities: pixel counts inherit the map’s errors, and the confusion matrix is what converts them back into an unbiased estimate.

---

Prerequisites

pip install "rasterio>=1.3.0" "numpy>=1.23" "scikit-learn>=1.3" "geopandas>=0.14" "scipy>=1.11"

You also need a prediction raster that is pixel-aligned with the reference, which is the invariant established in running model inference over large rasters, and reference data that was never seen during training — the fold discipline from building training datasets from satellite imagery.

---

Step-by-Step Workflow

Step 1 — Put the pair on one grid

Comparing two rasters that differ by half a pixel produces an accuracy figure dominated by registration error. Assert the grids match, or resample the reference onto the prediction grid with nearest-neighbour — never the other way round, because resampling the prediction changes the thing being assessed.

import numpy as np
import rasterio
from rasterio.enums import Resampling
from rasterio.warp import reproject


def load_aligned_pair(pred_path: str, ref_path: str) -> tuple[np.ndarray, np.ndarray]:
    with rasterio.open(pred_path) as pred, rasterio.open(ref_path) as ref:
        p = pred.read(1)
        if (ref.crs, ref.transform, ref.shape) == (pred.crs, pred.transform, pred.shape):
            r = ref.read(1)
        else:
            r = np.zeros(pred.shape, dtype=ref.dtypes[0])
            reproject(
                source=rasterio.band(ref, 1), destination=r,
                src_transform=ref.transform, src_crs=ref.crs,
                dst_transform=pred.transform, dst_crs=pred.crs,
                resampling=Resampling.nearest,   # class codes must not be averaged
            )
    return p, r

Step 2 — Sample, do not census

Using every reference pixel looks rigorous and is not. Reference maps have their own errors, their coverage is rarely uniform, and a wall-to-wall comparison weights whichever region happens to be largest. A stratified random sample of a few hundred points per class gives tighter, more honest intervals from far less data.

import numpy as np


def stratified_sample(pred: np.ndarray, n_per_class: int = 375,
                      seed: int = 0) -> tuple[np.ndarray, np.ndarray]:
    """Draw equal numbers of pixel indices from each mapped class."""
    rng = np.random.default_rng(seed)
    rows, cols = [], []
    for cls in np.unique(pred[pred > 0]):
        idx = np.flatnonzero(pred.reshape(-1) == cls)
        take = rng.choice(idx, size=min(n_per_class, idx.size), replace=False)
        r, c = np.unravel_index(take, pred.shape)
        rows.append(r)
        cols.append(c)
    return np.concatenate(rows), np.concatenate(cols)

Equal allocation per class is the design that makes rare classes reportable. It does mean the raw sample is not representative of the map, which is exactly why the area estimator in step 5 re-weights each stratum by its mapped area.

Step 3 — Build the confusion matrix

With a sample in hand the matrix is a cross-tabulation. Writing it with np.bincount rather than a library call keeps it fast enough to run over millions of pixels when you do want a census.

import numpy as np


def confusion(pred: np.ndarray, ref: np.ndarray, n_classes: int) -> np.ndarray:
    """Rows = predicted class, columns = reference class, 1-based codes."""
    valid = (pred > 0) & (ref > 0)
    flat = (pred[valid].astype(np.int64) - 1) * n_classes + (ref[valid].astype(np.int64) - 1)
    return np.bincount(flat, minlength=n_classes ** 2).reshape(n_classes, n_classes)

Fixing the orientation once and documenting it prevents an entire genre of reporting error: with rows as predictions, dividing a row by its sum gives user’s accuracy (of the pixels the map calls class k, what fraction really are), and dividing a column by its sum gives producer’s accuracy (of the pixels that really are class k, what fraction the map found).

Step 4 — Derive metrics per class

import numpy as np


def class_metrics(cm: np.ndarray) -> dict[str, np.ndarray]:
    diag = np.diag(cm).astype("float64")
    row = cm.sum(axis=1).astype("float64")     # predicted totals
    col = cm.sum(axis=0).astype("float64")     # reference totals
    with np.errstate(divide="ignore", invalid="ignore"):
        users = np.divide(diag, row, out=np.full_like(diag, np.nan), where=row > 0)
        producers = np.divide(diag, col, out=np.full_like(diag, np.nan), where=col > 0)
        f1 = 2 * users * producers / (users + producers)
    return {
        "users_accuracy": users,
        "producers_accuracy": producers,
        "f1": f1,
        "overall": diag.sum() / cm.sum(),
    }

Report all three per class. A model with 0.95 user’s accuracy and 0.40 producer’s accuracy for “water” is not a good water map — it finds almost no water, but what it finds is water. The single overall figure hides exactly that.

Step 5 — Estimate area, not pixel count

The stratified estimator treats each mapped class as a stratum, uses the sample to estimate what each stratum really contains, and weights by the stratum’s mapped area.

import numpy as np


def area_estimates(cm: np.ndarray, mapped_px: np.ndarray, px_area_ha: float):
    """Stratified area estimate and standard error from a per-stratum sample."""
    w = mapped_px / mapped_px.sum()                     # stratum weights
    n = cm.sum(axis=1).astype("float64")                # sample size per stratum
    p = w[:, None] * cm / np.where(n[:, None] > 0, n[:, None], 1)   # cell proportions

    prop = p.sum(axis=0)                                # estimated class proportions
    var = ((w[:, None] ** 2) * (cm / np.where(n[:, None] > 0, n[:, None], 1))
           * (1 - cm / np.where(n[:, None] > 0, n[:, None], 1))
           / np.where(n[:, None] > 1, n[:, None] - 1, 1)).sum(axis=0)

    total_ha = mapped_px.sum() * px_area_ha
    return prop * total_ha, 1.96 * np.sqrt(var) * total_ha

The returned interval is what turns a map into a measurement. “Built-up area is 48,200 ha” is a claim about a picture; “built-up area is 41,300 ± 5,800 ha” is a claim about the ground.

---

Cross-Validation Designs That Match the Question

Every validation design answers a specific question, and choosing the wrong one produces a number that answers a question nobody asked. Before splitting anything, decide which of three questions you actually need answered: how well does the model fill gaps inside the area it was trained on, how well does it work in an unvisited part of the same landscape, or how well does it transfer to a different region entirely.

A random split answers the first question, and only the first — and even then it is optimistic, because gaps in practice are rarely single pixels. A blocked split answers the second, provided the blocks are larger than the range over which the target is autocorrelated. Leave-one-region-out answers the third, and it is the only design whose number you should quote when someone asks whether the model will work in their country.

The buffered variant deserves more use than it gets. Blocked splits still put training cells immediately adjacent to validation cells along every block boundary, and with a 5×5 focal feature those two cells literally share input pixels. Discarding a one-block strip between roles costs a modest fraction of the data and removes the leak entirely:

import numpy as np
from sklearn.model_selection import GroupKFold


def buffered_groups(block_x: np.ndarray, block_y: np.ndarray,
                    n_folds: int = 5, buffer: int = 1) -> np.ndarray:
    """Group ids with a discarded buffer ring; -1 marks a sample to drop."""
    fold = (block_x // (buffer + 1) + block_y // (buffer + 1)) % n_folds
    on_buffer = ((block_x % (buffer + 1)) != 0) | ((block_y % (buffer + 1)) != 0)
    return np.where(on_buffer, -1, fold)

Whichever design you pick, report it alongside the number. “Overall accuracy 0.83 under five-fold leave-one-region-out with 20 km blocks” is a statement someone can evaluate; “overall accuracy 0.83” is not.

---

Parameter Reference

---

Verification & Testing

The first check is that the matrix and the map agree on how much of the scene each class covers.

import numpy as np

assert cm.sum() == sample_size, "matrix lost samples to nodata"
assert np.isclose(prop.sum(), 1.0, atol=1e-6), "proportions must sum to one"
assert (se >= 0).all(), "negative variance means a stratum had one sample"

The second is a null-model comparison. Compute the accuracy of a map that assigns every pixel to the most common class; if the real model is not comfortably above it, the model has learned the prior and nothing else. For a scene that is 52% crops, a majority-class map scores 0.52, and a model at 0.61 is much weaker than it sounds.

The third is spatial: map the errors. Rasterise the sample points with a correct/incorrect flag and look at where the failures concentrate. Errors clustered in one valley or on one date are a data problem; errors spread uniformly are a model problem. This diagnostic changes what you do next more often than any aggregate statistic does.

---

Troubleshooting

Overall accuracy is suspiciously high, above 0.97

Check the split before celebrating. Either validation pixels came from the same fields as training pixels, or the reference layer is the same product that generated the training labels. Both produce this signature. Re-run with a spatially blocked design from designing spatial cross-validation for raster models.

A class has producer’s accuracy of zero

The model never predicts it. Usually the class is too rare in training, or two classes are spectrally identical and the model has collapsed them. Check the confusion matrix row for where its pixels went — they are always somewhere specific.

RuntimeWarning: invalid value encountered in divide

A class has no samples in either the row or the column. Guard the divisions as the code above does, and report nan rather than silently propagating a zero that reads as “0% accurate”.

Estimated area is wildly different from mapped area

That is the estimator working, up to a point — but a factor-of-two difference means the map has a systematic commission problem for that class. Look at the off-diagonal cells: the pixels have been taken from somewhere identifiable.

Confidence intervals are enormous for one class

Its stratum had too few samples, or the class is genuinely rare and heterogeneous. Increase allocation for that stratum specifically rather than raising the sample everywhere; the estimator handles unequal allocation natively.

---

Frequently Asked Questions

Q: Why is counting classified pixels a biased area estimate? Because classification errors are asymmetric. A class that is over-predicted gains pixels it should not have and loses none, so its pixel count is biased upward by exactly the commission error. The stratified estimator corrects this using the confusion matrix, and it also gives a standard error, which a pixel count cannot.

Q: What overall accuracy should I expect from a land cover model? For a handful of spectrally distinct classes over a single scene, 85 to 92 per cent is routine. For ten or more classes across a region with several ecological zones, 70 to 80 per cent is a good result. Anything above 95 per cent on a multi-class problem is usually evidence of a leaky validation split rather than a strong model.

Q: Is kappa a useful metric for raster accuracy? Rarely. Kappa corrects overall accuracy for chance agreement, but the correction depends on class prevalence in a way that makes comparisons across maps meaningless. Report per-class producer and user accuracy with confidence intervals instead; they carry the information a reader needs.

Q: How do I assess a continuous prediction rather than classes? Use the same sampling design with regression metrics: mean absolute error and bias per stratum, plus a scatter of predicted against reference with a one-to-one line. Bias by stratum is the diagnostic that matters, because a model can have excellent overall error while being systematically low in the range you care about.

Q: What if my reference data is itself imperfect? Then the assessment measures agreement, not truth, and it should say so. Quantify the reference error where you can — double-interpret a subset of points and report the interpreter agreement rate — and treat the resulting accuracy as an upper bound on how far the map can be trusted. A map cannot score higher than its reference layer is reliable.

---

Related

• Computing Confusion Matrices from Raster Pairs — the aligned-pair cross-tabulation in full.
• Designing Spatial Cross-Validation for Raster Models — blocks, buffers and leave-one-region-out.
• Estimating Area with Stratified Random Sampling — the estimator, its assumptions, and how to report it.
• Measuring Boundary Accuracy with IoU and F1 — object-level metrics for segmentation output.
• Thresholding Change Maps and Removing Noise — the same assessment logic applied to change products.