Spatial Pipeline EngineeringMRV Automation

Propagating Spatial Autocorrelation into Uncertainty Budgets

Aggregating a carbon estimate over an area is the last arithmetic step in most MRV pipelines and the one most likely to be wrong by a large factor. The reason is that the standard error of a mean shrinks as the square root of the sample size only when the samples are independent, and pixels in a carbon map are emphatically not. Neighbouring pixels share terrain, share climate, share management history, and — critically — share the model that produced them, so a million-pixel project does not contain a million independent observations. This guide shows how to account for that, within emission factor uncertainty mapping in the spatial modeling and carbon stock validation stack.

The practical consequence is stark. A project reporting a total stock with a confidence interval computed from independent pixel variance will typically report an interval one to two orders of magnitude too narrow. That interval then feeds an uncertainty deduction, and the deduction comes out near zero — which is exactly the outcome that makes a reviewer look closely, and exactly the outcome the arithmetic guarantees regardless of the data.

Independent pixels versus correlated pixels, and what happens to the interval Two panels side by side over the same project area. On the left, the independence assumption: one hundred thousand pixels are treated as one hundred thousand observations, the standard error divides by the square root of that number, and the resulting confidence interval on total stock is plus or minus a fraction of a percent. On the right, the same area with a correlation range of eight hundred metres: the number of effectively independent observations falls to a few hundred blocks, the standard error divides by the square root of that much smaller number, and the interval widens by a factor of roughly twenty. A panel notes that the data did not change, only the assumption about it, and that the wide interval is the honest one. Same map, same pixels, two assumptions Only one of them is defensible, and it is the one that widens the interval. Assumed independent n = 100,000 observations SE = σ / √100,000 interval ≈ ±0.4% Nothing in the pipeline objects. The uncertainty deduction comes out near zero, and the reviewer starts asking questions. Correlated, range 800 m n_eff ≈ 300 blocks SE = σ / √300 interval ≈ ±7% The data did not change. Only the assumption about how much independent information it contains.

Root Cause Analysis

Spatial autocorrelation in a carbon map has three distinct sources, and they need separating because they behave differently under aggregation.

The landscape itself is autocorrelated. Soil type, elevation, rainfall, and disturbance history all vary smoothly, so two nearby hectares genuinely resemble each other more than two distant ones. This is real structure in the world and it is the component a variogram is designed to describe. It has a finite range: beyond some distance, typically hundreds of metres to a few kilometres depending on the biome, the resemblance disappears.

The model that produced the map is autocorrelated in its errors. This is the component most often forgotten. A biomass model that systematically underpredicts in dense canopy produces errors that cluster wherever canopy is dense, and that clustering has whatever spatial pattern the canopy has. Crucially, model error correlation does not necessarily decay with distance at all — a bias affecting an entire forest type is perfectly correlated across every pixel of that type, however far apart they are. Averaging over more area does not reduce it.

The inputs are shared. Every pixel derived from a single satellite scene inherits that scene’s atmospheric correction, its calibration, and its geolocation. A single emission factor applied across a stratum makes every pixel in that stratum share one number, and the uncertainty in that number is fully correlated across all of them. This component is often the largest and it is the easiest to handle, because it does not require a variogram — a factor applied to the whole area propagates as a proportional uncertainty on the total, unreduced by area.

The failure this produces is systematic in one direction. Every one of the three sources makes the effective sample size smaller than the pixel count, so treating pixels as independent always understates the interval and never overstates it. That one-sided bias is why the assumption is not a neutral simplification.

Diagnostic Pipeline / Pre-Flight Validation

Before computing any interval, establish the correlation structure empirically and check that the aggregation is not operating in a regime where the arithmetic breaks down. The empirical variogram is the standard tool and it is worth computing even when a model-based approach will follow, because its shape is diagnostic.

import math
from dataclasses import dataclass

import numpy as np
import structlog

log = structlog.get_logger()


@dataclass(frozen=True)
class VariogramPoint:
    lag_m: float
    semivariance: float
    n_pairs: int


@dataclass(frozen=True)
class CorrelationStructure:
    """Fitted spherical variogram parameters, in the units of the variable."""
    nugget: float
    sill: float
    range_m: float
    fit_rmse: float

    @property
    def nugget_fraction(self) -> float:
        """Share of total variance that is uncorrelated at the shortest lag."""
        return self.nugget / self.sill if self.sill > 0 else 1.0


def empirical_variogram(
    x: np.ndarray,
    y: np.ndarray,
    values: np.ndarray,
    *,
    lag_m: float,
    n_lags: int,
    max_pairs: int = 2_000_000,
) -> list[VariogramPoint]:
    """Binned semivariance against separation distance.

    Subsamples when the full pair set would be prohibitive. Subsampling
    is unbiased for the variogram; it only widens the estimate's own
    uncertainty, which is reported through n_pairs.
    """
    n = len(values)
    if n * (n - 1) // 2 > max_pairs:
        rng = np.random.default_rng(seed=0)
        keep = rng.choice(n, size=int(math.sqrt(2 * max_pairs)), replace=False)
        x, y, values = x[keep], y[keep], values[keep]
        log.info("variogram.subsampled", kept=len(keep), original=n)

    dx = x[:, None] - x[None, :]
    dy = y[:, None] - y[None, :]
    dist = np.sqrt(dx * dx + dy * dy)
    diff_sq = (values[:, None] - values[None, :]) ** 2

    iu = np.triu_indices(len(values), k=1)
    dist, diff_sq = dist[iu], diff_sq[iu]

    points: list[VariogramPoint] = []
    for i in range(n_lags):
        lo, hi = i * lag_m, (i + 1) * lag_m
        sel = (dist >= lo) & (dist < hi)
        count = int(sel.sum())
        if count < 30:
            # Too few pairs for a stable estimate; report, do not fabricate.
            log.warning("variogram.sparse_lag", lag_m=(lo + hi) / 2, n_pairs=count)
            continue
        points.append(
            VariogramPoint(
                lag_m=(lo + hi) / 2,
                semivariance=float(diff_sq[sel].mean() / 2),
                n_pairs=count,
            )
        )
    return points


def assert_structure_usable(cs: CorrelationStructure, extent_m: float) -> None:
    """Refuse structures that cannot support an aggregation.

    Two regimes break the arithmetic: a range comparable to the extent,
    where the whole area is one correlated blob and there is effectively
    one observation; and a pure nugget, where the variogram found no
    structure at all and is probably being computed on residual noise.
    """
    if cs.range_m > extent_m / 3:
        raise ValueError(
            f"correlation range {cs.range_m:.0f} m exceeds a third of the "
            f"{extent_m:.0f} m extent — the area contains too few independent "
            "blocks for a mean to be meaningful; report the estimate at a "
            "coarser unit or widen the area"
        )

    if cs.nugget_fraction > 0.95:
        raise ValueError(
            f"nugget is {cs.nugget_fraction:.0%} of the sill — no spatial "
            "structure detected; check that the variable is the model output "
            "rather than its residual, and that coordinates are projected"
        )

    log.info(
        "variogram.accepted",
        range_m=cs.range_m,
        nugget_fraction=round(cs.nugget_fraction, 3),
        fit_rmse=round(cs.fit_rmse, 4),
    )

The check that fires most usefully is the first. A project whose correlation range is comparable to its own extent is not a project with a precise mean and a wide interval — it is a project with one observation, and the correct response is to say so rather than to report a number with a confidence interval that implies replication.

A variogram and what each of its three parameters costs the uncertainty budget A variogram curve rising from a non-zero intercept at the origin to a plateau. The intercept is labelled nugget and annotated as the uncorrelated component, which is the only part that averages away with more pixels. The plateau is labelled sill and annotated as the total variance. The distance at which the curve reaches the plateau is labelled range and annotated as the block size for effective sample size: pixels closer than this share information. A shaded region below the nugget line is labelled the only variance that shrinks with area, and a second annotation notes that everything above it persists no matter how many pixels are averaged. Only the nugget averages away The rest of the variance is shared between neighbours and survives aggregation. sill nugget range separation distance → semivariance total variance what one pixel varies by block size for n_eff pixels closer than this share their information shaded: shrinks with area everything above it does not

Deterministic Transformation Logic

The aggregation itself replaces the pixel count with an effective sample size and then adds back the components that do not decay with distance at all. Three terms, computed separately and combined in quadrature only when they are genuinely independent of one another.

import math
from dataclasses import dataclass


@dataclass(frozen=True)
class UncertaintyBudget:
    """The three components of an aggregated uncertainty, kept separate.

    They are kept separate because they respond differently to a change in
    area, and a verifier asking 'what would halve this' needs to see which
    term dominates.
    """
    area_ha: float
    mean_estimate: float
    se_random: float        # decays with effective sample size
    se_model_bias: float    # does not decay with area at all
    se_factor: float        # proportional, shared across the whole stratum

    @property
    def se_total(self) -> float:
        return math.sqrt(
            self.se_random ** 2 + self.se_model_bias ** 2 + self.se_factor ** 2
        )

    @property
    def relative_half_width_95(self) -> float:
        return 1.96 * self.se_total / self.mean_estimate

    @property
    def dominant_term(self) -> str:
        terms = {
            "random": self.se_random,
            "model_bias": self.se_model_bias,
            "emission_factor": self.se_factor,
        }
        return max(terms, key=terms.get)


def effective_sample_size(
    n_pixels: int, pixel_m: float, range_m: float, nugget_fraction: float
) -> float:
    """Independent observations implied by the correlation structure.

    Pixels within one correlation range of each other contribute roughly one
    observation between them. The nugget fraction restores the share of
    variance that is genuinely independent at pixel scale, so a high-nugget
    field keeps more of its nominal sample size.
    """
    if range_m <= pixel_m:
        return float(n_pixels)

    pixels_per_block = (range_m / pixel_m) ** 2
    n_blocks = n_pixels / pixels_per_block

    # Blend: the nugget share behaves independently, the structured share
    # behaves at block scale.
    return nugget_fraction * n_pixels + (1 - nugget_fraction) * n_blocks


def aggregate_with_correlation(
    *,
    pixel_values: list[float],
    pixel_m: float,
    cs: CorrelationStructure,
    model_bias_rel: float,
    factor_rel: float,
) -> UncertaintyBudget:
    """Aggregate a carbon field to a total with a defensible interval.

    model_bias_rel and factor_rel are relative standard errors that do NOT
    shrink with area — a systematic model bias and a shared emission factor
    respectively. Passing zero for either is a claim that the model is
    unbiased and the factor exact, which no pipeline should assert silently.
    """
    if model_bias_rel <= 0 or factor_rel <= 0:
        raise ValueError(
            "model bias and emission factor uncertainty must be positive; "
            "zero asserts a perfect model and an exact factor"
        )

    n = len(pixel_values)
    mean = sum(pixel_values) / n
    var = sum((v - mean) ** 2 for v in pixel_values) / (n - 1)

    n_eff = effective_sample_size(n, pixel_m, cs.range_m, cs.nugget_fraction)
    area_ha = n * (pixel_m ** 2) / 10_000
    total = mean * area_ha

    return UncertaintyBudget(
        area_ha=area_ha,
        mean_estimate=total,
        se_random=math.sqrt(var / n_eff) * area_ha,
        se_model_bias=total * model_bias_rel,
        se_factor=total * factor_rel,
    )

The refusal to accept zero for the two non-decaying terms is deliberate and it changes behaviour more than the variogram does. In a large project the random term is often not the dominant one even after correction — model bias and emission factor uncertainty are, and neither of them cares how many pixels there are. A pipeline that models only the random term produces a budget that gets arbitrarily tight as the project grows, which is obviously wrong and yet is what most implementations do.

Compliance Gating & Audit Trail Generation

The audit record for an aggregated figure needs four things, and a budget object with the terms kept separate supplies all of them.

The variogram itself, with its fitted parameters and the number of pairs behind each lag. A range fitted from three sparse lags is not evidence, and a verifier who cannot see the pair counts cannot tell the difference between a fitted structure and a plausible one.

The effective sample size and the pixel count together. The ratio between them is the single number that says how much the correction mattered, and stating it pre-empts the question of whether autocorrelation was considered at all.

Which term dominates the budget. This is what a reviewer needs to assess whether the project’s proposed improvements are the ones that would help — buying more imagery does not reduce an emission-factor term, and refitting the model does not reduce a term driven by a factor from a national database.

Whether the correlation range was estimated on this project’s data or imported. Importing a range from a published study of a different biome is sometimes the only option, and it is acceptable when disclosed; it is not acceptable when a fitted-looking number in a report turns out on questioning to have come from a paper about a different continent.

Production Integration

The natural place for this logic is immediately before the reporting step and immediately after the per-pixel Monte Carlo, so that the pixel-level uncertainty from Monte Carlo uncertainty propagation for emission factors feeds the aggregation rather than being replaced by it. The two are complements: the Monte Carlo establishes what a pixel’s uncertainty is, and this step establishes how much of it survives averaging.

One operational detail matters more than it looks. The variogram should be computed on the quantity being aggregated, not on its residuals against field data, unless the intent is specifically to characterise model error correlation — in which case that is a separate variogram feeding the model bias term rather than the random one. Conflating the two is common and it produces a range that describes neither.

How each uncertainty term responds when the project area grows tenfold Three horizontal bars showing the relative size of the three uncertainty terms at a small project area and at ten times that area. The random sampling term shrinks substantially when the area grows, because more independent blocks are averaged. The model bias term does not change at all, because a systematic bias affects every pixel in the same direction regardless of how many there are. The shared emission factor term also does not change, because one factor applies across the whole stratum. A panel notes that as a project grows the budget becomes dominated by the two terms that area cannot help, and that a pipeline modelling only the random term produces a total uncertainty that falsely approaches zero. What growing the project tenfold actually buys Two of the three terms do not move. Random sampling Model bias Emission factor small area 10× area ±6.0% ±1.9% — √n_eff helps here ±4.3% ±4.3% — unchanged ±3.4% ±3.4% — unchanged

Frequently Asked Questions

Is a variogram always necessary, or is a fixed block size acceptable?

A fixed block size is acceptable as a conservative stand-in provided it is chosen larger than any plausible correlation range for the biome and the choice is disclosed. It costs precision — the interval comes out wider than a fitted structure would give — but it errs in the safe direction and it removes a fitting step that can go wrong. What is not acceptable is a fixed block size chosen small enough to be convenient, because that reintroduces the original error in a form that looks like it was considered.

What if the variogram does not reach a sill within the project extent?

Then the project has no scale at which its pixels become independent, and the mean over the project is one observation rather than an average of many. This happens with strong regional gradients — a rainfall or elevation trend running across the site. The usual fix is to detrend first: fit and remove the large-scale trend, compute the variogram on the residual, and carry the trend’s own uncertainty as a separate term. Reporting a mean with a narrow interval over an undetrended gradient is the failure mode this avoids.

Does this apply to area-based activity data as well as to stock?

Yes, and it is often overlooked there because area feels like it is counted rather than estimated. A deforestation area derived from a classified map has classification errors that cluster — confusion between forest and shrubland concentrates in transitional zones — so the area estimate has a correlated error structure exactly like a stock estimate. The standard approach in that setting is a design-based estimator with a stratified reference sample, which handles the correlation implicitly by sampling rather than by modelling it.

How does this interact with the uncertainty deduction a methodology requires?

Directly, and usually in the project’s disfavour, which is why the arithmetic gets attention. Most crediting methodologies deduct a share of the estimate as a function of the relative half-width of the confidence interval, with a threshold below which no deduction applies. Correcting for autocorrelation typically moves a project from below that threshold to above it. The correct response is to reduce the uncertainty rather than the estimate of it — better calibration, project-specific factors, more field plots — and the budget’s dominant-term field says which of those would actually work.

Can more pixels ever compensate for correlation?

For the random term, yes, but with sharply diminishing returns: doubling the pixel count at fixed area does nothing at all once the pixel is smaller than the correlation range, because the extra pixels carry no new information. Increasing the area does help the random term, since it adds genuinely new blocks. Neither helps the model bias or emission factor terms. This is the practical reason to keep the terms separate — it makes visible that resampling a map to a finer grid, a common instinct, changes nothing.

Should model error correlation be a separate variogram or folded into the main one?

Separate, because they answer different questions and often have different ranges. The variogram of the predicted field describes how the landscape varies. The variogram of the residuals against field observations describes how the model’s errors cluster, and its range is frequently much longer — a bias tied to forest type persists across every stand of that type. Folding them together produces a single range that under-corrects for the model term and over-corrects for the landscape one.

What is a reasonable relative uncertainty to expect after correcting properly?

For a project-scale forest carbon stock with project-specific allometry and a decent field plot network, a 95% half-width in the region of eight to fifteen percent is typical and defensible. Below five percent is unusual and invites scrutiny about which terms were omitted. Above about twenty-five percent, methodological deductions usually become large enough that improving the estimate is cheaper than accepting the deduction. Those bands are heuristics from practice rather than requirements, and a specific methodology’s own thresholds always govern.