Estimating Area with Stratified Random Sampling

Use the mapped classes as strata, weight each by its mapped area, and let the reference sample correct the map’s errors:

import numpy as np

w = mapped_px / mapped_px.sum()                       # stratum weights
n = cm.sum(axis=1)                                    # samples per stratum
p = w[:, None] * cm / n[:, None]                      # estimated cell proportions
area_ha = p.sum(axis=0) * mapped_px.sum() * px_area_ha

The result differs from a pixel count by exactly the map’s commission and omission errors, and it comes with a variance. This page belongs to validating raster predictions and accuracy assessment in Raster Machine Learning & Model Inference.

Why the Pixel Count Is Always Wrong

A classified map makes two kinds of error for every class: it claims pixels that belong to something else, and it misses pixels that belong to it. Those two do not cancel. If a class is easy to over-predict — bright bare soil called built-up, say — its commission error exceeds its omission error and the pixel count is biased high, systematically and unavoidably.

What the pixel count actually counts The mapped extent of a class is made of correctly identified pixels plus commission errors. The true extent is the correct pixels plus the omission errors that the map missed. Only when commission and omission happen to be equal do the two totals agree, which is not a property any classifier guarantees. Mapped extent and true extent are different quantities mapped correct — 604 px commission — 161 true correct — 604 px omission — 41 Counting pixels gives 765; the ground holds 645 — an overestimate of 19%. The bias is the difference between the red bar and the green one.

The estimator below recovers the true extent from the confusion matrix, which is why the sample has to be a probability sample: every pixel in a stratum must have had a known chance of being selected, or the weights mean nothing.

Environment & Setup

Package Version pin Used for
numpy >=1.23 The estimator and its variance
rasterio >=1.3.0 Stratum pixel counts and pixel area from the transform
pandas >=2.0 Presenting estimates with intervals 
pip install "numpy>=1.23" "rasterio>=1.3.0" "pandas>=2.0"

Complete Working Example

import numpy as np
import pandas as pd
import rasterio


def stratum_counts(path: str, class_codes: list[int]) -> tuple[np.ndarray, float]:
    """Mapped pixel count per class and the area of one pixel in hectares."""
    with rasterio.open(path) as src:
        arr = src.read(1)
        px_area_ha = abs(src.transform.a * src.transform.e) / 10_000
    return np.array([(arr == c).sum() for c in class_codes], dtype="float64"), px_area_ha


def stratified_area(cm: np.ndarray, mapped_px: np.ndarray,
                    px_area_ha: float, class_codes: list[int],
                    z: float = 1.96) -> pd.DataFrame:
    """Unbiased area estimates with confidence intervals.

    cm: rows = mapped class (the stratum), columns = reference class.
    """
    w = mapped_px / mapped_px.sum()
    n = cm.sum(axis=1).astype("float64")
    safe_n = np.where(n > 0, n, 1.0)

    p_hat = cm / safe_n[:, None]                 # within-stratum proportions
    p = w[:, None] * p_hat                       # weighted cell proportions
    prop = p.sum(axis=0)                         # estimated class proportions

    var = ((w[:, None] ** 2) * p_hat * (1 - p_hat)
           / np.where(n > 1, n - 1, 1.0)[:, None]).sum(axis=0)
    total_ha = mapped_px.sum() * px_area_ha

    return pd.DataFrame({
        "mapped_ha": (mapped_px * px_area_ha).round(0),
        "estimated_ha": (prop * total_ha).round(0),
        "ci95_ha": (z * np.sqrt(var) * total_ha).round(0),
        "bias_pct": (100 * (mapped_px * px_area_ha - prop * total_ha)
                     / np.maximum(prop * total_ha, 1)).round(1),
    }, index=class_codes)

The bias_pct column is the one to look at first. A class whose mapped area exceeds its estimated area by more than about ten per cent has a commission problem worth investigating, and the confusion matrix row will say exactly which class its extra pixels were taken from.

Allocating the Sample

Proportional versus equal allocation Proportional allocation puts most points in the dominant class and leaves the rare class with a handful, so its interval is uselessly wide. Equal allocation gives every class enough points to be estimated, at a small cost to the precision of the overall figure. A compromise allocation, weighted but with a floor, is usually the right answer. proportional equal with a floor crops 780 forest 450 built 120 water 32 — interval is useless crops 420 forest 360 built 300 water 300 — usable The estimator re-weights by mapped area regardless of allocation, so unequal allocation is not a bias. It is purely a choice about which class gets the precision.

Because the weights in the estimator come from the map rather than from the sample, allocation does not bias anything — it only decides where the precision goes. That is a genuinely useful property: you can add a hundred more points to a troublesome class after the first round and simply recompute, without invalidating what came before.

A floor of roughly 250 to 300 points per stratum, with the remainder allocated proportionally, is a reasonable default. Below about 100 points a stratum’s cell proportions are so noisy that the interval swamps the estimate, and the honest presentation is to report the class as unmeasured rather than to publish a number with an interval wider than itself.

Verification

import numpy as np

assert np.isclose(prop.sum(), 1.0, atol=1e-6), "proportions must sum to one"
assert (var >= 0).all(), "negative variance: a stratum has one sample"
assert np.isclose((mapped_px * px_area_ha).sum(),
                  (prop * total_ha).sum(), rtol=1e-6), "total area must be preserved"
The estimator moves area between classes, never creates it Two stacked bars of equal total length represent the mapped and estimated composition of the study area. Between them, built-up shrinks and cropland grows, but the overall length is identical because every pixel belongs to exactly one class in both accountings. Total area is conserved; composition is not mapped crops forest built estimated crops forest built Both bars end on the dashed line — if yours do not, a stratum weight is wrong.

The last assertion is the most useful and the least obvious: the estimator redistributes area between classes but must conserve the total, because every pixel is in exactly one class in both the map and the reference. It is worth asserting rather than eyeballing: the failure is a quiet few per cent, not an obvious factor of two, and a few per cent is exactly the size of the effect the estimator exists to measure. If the totals disagree, a stratum weight is wrong — usually because the pixel counts were taken from a different file or a different nodata convention from the one the sample was drawn against.

Common Errors

Confidence intervals wider than the estimate

The stratum has too few samples, or the class is genuinely rare and heterogeneous. Add points to that stratum specifically; the estimator handles unequal allocation natively.

Estimated areas sum to more than the study area

Stratum weights were computed from a different raster than the one the sample came from, or nodata pixels were included in one count and not the other. Recompute both from the same file.

The estimate barely differs from the pixel count

That is a good sign, not a bug — it means commission and omission are balanced for that class. Check the confusion matrix to confirm both are small rather than large and offsetting.

Frequently Asked Questions

Q: Why not simply multiply pixel count by pixel area? Because that inherits every commission error in the map. A class the model over-predicts gains pixels it should not have, and no amount of careful pixel counting removes that bias. The stratified estimator uses the reference sample to correct it and returns a standard error as well.

Q: How many reference points do I need? A few hundred per stratum for a useful per-class interval, weighted toward the rare classes that matter. Overall precision improves slowly with sample size, so doubling from 500 to 1000 narrows the interval by only about thirty per cent.

Q: Can I use the same sample for accuracy and for area? Yes, and you should — they come from the same confusion matrix. The accuracy figures describe the map’s quality and the area figures describe the ground, and both follow from the same stratified sample without extra fieldwork.