Computing Confusion Matrices from Raster Pairs

Encode the predicted and reference class of every valid pixel as one integer, then count them in a single pass:

import numpy as np

valid = (pred > 0) & (ref > 0)
codes = (pred[valid].astype(np.int64) - 1) * n_classes + (ref[valid].astype(np.int64) - 1)
cm = np.bincount(codes, minlength=n_classes ** 2).reshape(n_classes, n_classes)
# rows = predicted, columns = reference

That is roughly twenty times faster than a library cross-tabulation on tens of millions of pixels, and it makes the row and column convention explicit rather than implied. This page belongs to validating raster predictions and accuracy assessment in Raster Machine Learning & Model Inference.

Reading the Matrix Correctly

The matrix contains more information than any summary derived from it, and most of that information is in the off-diagonal cells — they name which classes are being confused with which, which is what tells you what to fix.

Anatomy of a confusion matrix A four by four matrix with predictions as rows and reference as columns. The diagonal holds correct pixels. Reading across a row shows what the map called that class when it was something else, which is commission error. Reading down a column shows where the true class ended up instead, which is omission error. Rows: what the map said · Columns: what it really was crops 812 41 18 4 forest 63 740 9 2 built 128 22 604 11 water 3 6 14 683 crops forest built water the orange cell is the story 128 pixels mapped as built-up are really cropland — bare fields read as bright impervious surface Built-up user's accuracy is 604/765 = 0.79; its producer's accuracy is 604/645 = 0.94. The map finds nearly all the built-up there is, and calls a lot of bare soil built-up as well.

That asymmetry is invisible in an overall accuracy figure and decisive for anyone using the product. A planner who needs every settlement found is well served; an analyst measuring impervious area is being handed a systematic overestimate.

Environment & Setup

Package Version pin Used for
rasterio >=1.3.0 Reading both rasters and checking alignment
numpy >=1.23 The bincount tabulation and metric arithmetic
pandas >=2.0 Rendering the matrix with class names 
pip install "rasterio>=1.3.0" "numpy>=1.23" "pandas>=2.0"

Complete Working Example

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


def load_aligned(pred_path: str, ref_path: str) -> tuple[np.ndarray, np.ndarray]:
    """Return prediction and reference on the PREDICTION grid."""
    with rasterio.open(pred_path) as pred, rasterio.open(ref_path) as ref:
        p = pred.read(1)
        same = ((ref.crs, ref.transform, ref.shape)
                == (pred.crs, pred.transform, pred.shape))
        if same:
            r = ref.read(1)
        else:
            r = np.zeros(pred.shape, dtype=ref.dtypes[0])
            reproject(rasterio.band(ref, 1), r,
                      src_transform=ref.transform, src_crs=ref.crs,
                      dst_transform=pred.transform, dst_crs=pred.crs,
                      resampling=Resampling.nearest)
    return p, r


def confusion_matrix(pred: np.ndarray, ref: np.ndarray,
                     class_codes: list[int]) -> pd.DataFrame:
    """Cross-tabulation with rows = predicted, columns = reference."""
    n = len(class_codes)
    lookup = np.zeros(max(class_codes) + 1, dtype=np.int64) - 1
    for i, code in enumerate(class_codes):
        lookup[code] = i

    valid = np.isin(pred, class_codes) & np.isin(ref, class_codes)
    pi = lookup[pred[valid].astype(np.int64)]
    ri = lookup[ref[valid].astype(np.int64)]

    counts = np.bincount(pi * n + ri, minlength=n * n).reshape(n, n)
    return pd.DataFrame(counts, index=class_codes, columns=class_codes)


def metrics(cm: pd.DataFrame) -> pd.DataFrame:
    m = cm.to_numpy().astype("float64")
    diag, row, col = np.diag(m), m.sum(axis=1), m.sum(axis=0)
    with np.errstate(divide="ignore", invalid="ignore"):
        users = np.divide(diag, row, out=np.full_like(diag, np.nan), where=row > 0)
        prods = np.divide(diag, col, out=np.full_like(diag, np.nan), where=col > 0)
    return pd.DataFrame({
        "n_predicted": row.astype(int),
        "n_reference": col.astype(int),
        "users_accuracy": users.round(3),
        "producers_accuracy": prods.round(3),
        "f1": (2 * users * prods / (users + prods)).round(3),
    }, index=cm.index)

The lookup table rather than a dictionary comprehension is what keeps this fast on a full scene: it turns the class-code translation into a single vectorised index operation instead of forty million Python-level lookups.

Alignment Is the Precondition

A half-pixel offset is an accuracy tax When prediction and reference grids are offset by half a pixel, every boundary between two classes produces a line of disagreeing pixels. For a landscape of small fields this can account for several percentage points of apparent error that has nothing to do with the model. grids aligned half-pixel offset boundary pixels agree a line of false errors on every boundary In a landscape of one-hectare fields this alone can cost three to five points of overall accuracy.

Always resample the reference onto the prediction grid, never the reverse. Resampling the prediction changes the object under assessment, which means the accuracy figure no longer describes the file you are going to publish. Use nearest-neighbour, because any interpolation of class codes invents values.

If the reference is vector rather than raster, burn it onto the prediction grid directly with the prediction’s own transform, exactly as described in rasterizing labels into segmentation masks. That avoids the intermediate raster and its own alignment question entirely.

Verification

assert cm.to_numpy().sum() > 0, "no overlapping valid pixels"
assert cm.shape[0] == cm.shape[1] == len(class_codes)
print("overall accuracy:", np.diag(cm.to_numpy()).sum() / cm.to_numpy().sum())
print((cm.to_numpy().sum(axis=1) == 0).sum(), "classes never predicted")
The matrix covers the intersection, not the map The prediction footprint and the reference footprint overlap only partly. The confusion matrix is built from the intersection, so the accuracy figure describes that intersection rather than the whole prediction. Reporting the fraction of the prediction that the assessment actually covers is part of reporting the accuracy. What the accuracy figure actually covers prediction reference only assessed intersection Report it: "overall accuracy 0.83, assessed over 24% of the mapped area".

A class that is never predicted shows as an all-zero row, and it is worth failing a report over rather than quietly publishing a matrix with a missing class. The complementary check — an all-zero column — means the reference contains no examples of that class, so its accuracy is unmeasured rather than perfect, and any summary that averages over classes must exclude it explicitly.

One more check catches a subtle and common mistake: compare the total valid count against the number of pixels valid in the prediction alone. If the joint mask drops more than a few per cent, the reference footprint is smaller than the prediction footprint and the accuracy applies only to the overlap — which is a legitimate result, but one that must be stated rather than assumed away. The distinction matters most when the reference data was collected somewhere convenient — along roads, near towns, inside one project’s study area — because then the assessed subset is not a random sample of the map and the figure cannot be extrapolated to the rest of it without an argument.

Common Errors

The matrix diagonal is nearly empty

Rows and columns are transposed relative to the class codes, or the two rasters use different code conventions — one 0-based and one 1-based. Print the unique values of both arrays before tabulating.

IndexError in the lookup table

A class value larger than any declared code is present, usually a stray 255 from an ignore region. Mask ignore values out before tabulating rather than extending the matrix to hold them.

Accuracy changes when the window changes

The comparison is over a subset with a different class mix, which is expected — accuracy is not a property of the model alone but of the model and the landscape together. Report the area the figure covers, and if two windows must be compared, compare the per-class producer and user accuracies rather than the overall figure, since those are far less sensitive to how the class mix happens to fall.

Frequently Asked Questions

Q: Rows predicted or rows reference? Pick one and document it in the function docstring. With rows as predictions, a row normalised by its sum gives user’s accuracy and a column normalised by its sum gives producer’s accuracy. Swapping the convention halfway through an analysis is the commonest source of contradictory accuracy claims.

Q: Should I build the matrix over every pixel? Only for a map-to-map comparison. For accuracy against reference data, a stratified sample gives better per-class precision and an honest variance estimate, because a census weights whatever class happens to dominate the scene.

Q: How do I handle a class that appears in one raster only? Size the matrix from the declared class list rather than from the values present, so an absent class appears as an empty row or column rather than silently shifting every index after it.