Spatial Pipeline EngineeringMRV Automation

Modeling Soil Organic Carbon with Digital Soil Mapping

Digital soil mapping turns a few hundred laboratory measurements into a wall-to-wall soil carbon surface using environmental covariates as the interpolating structure. Done carefully it is the only affordable way to map soil carbon over a project area; done carelessly it produces a beautiful raster whose stated accuracy is a fiction. This guide implements the careful version, within soil organic carbon modeling and validation in the spatial modeling and carbon stock validation stack.

The framing that keeps a pipeline honest is that the model is an interpolator over a covariate space, not a mechanistic account of soil formation. It can only be trusted where the covariate combinations it is asked to predict resemble those it was trained on. Most catastrophic digital soil maps are catastrophic in exactly those regions of covariate space where no core was ever collected — and nothing in a conventional accuracy report reveals this, which is why extrapolation detection is treated here as a first-class output rather than an optional extra.

Digital soil mapping pipeline with the two outputs that are usually missing A pipeline from a covariate stack and a core dataset through spatially blocked hyperparameter tuning to a fitted quantile regression forest. The model produces three outputs rather than one: a median prediction surface, a per-pixel prediction interval from the 5th and 95th quantiles, and an extrapolation mask marking pixels whose covariate combination lies outside the convex region the cores covered. A note states that the second and third outputs are usually absent from published maps and are the two a verifier asks for first. Covariate stack terrain · climate · land use bare-soil composite Cores SOC % · one lab method georeferenced, dated Blocked tuning GroupKFold over blocks never a random split Quantile forest full conditional distribution per pixel Median surface the map everyone ships Prediction interval q05 / q95 per pixel Extrapolation mask outside covariate hull usually missing the two a verifier asks for first

Root Cause Analysis

Three properties of soil data break the assumptions that make ordinary supervised learning trustworthy, and each needs an explicit countermeasure.

Samples are clustered, so the effective sample size is far below the row count. Four cores from one field corner are, for covariate purposes, close to one observation. A model trained on 400 such cores may have an effective sample size nearer 80, which changes both the achievable complexity and the honest error estimate. Random cross-validation hides this completely because it splits the cluster; blocked cross-validation exposes it, which is why the blocked score is often shocking on first sight and is nonetheless the true one.

The covariate stack must be temporally coherent with the sampling. A bare-soil composite built from acquisitions spanning three years, joined to cores collected in one autumn, embeds a mismatch between what the covariate describes and what the laboratory measured. Terrain is static and safe; climate normals are safe; anything derived from imagery must be windowed to the sampling campaign, and the window must be recorded. A model that appears to improve when you widen the imagery window is usually improving by fitting noise that happens to correlate with sampling location.

Extrapolation is invisible in the accuracy report. A random forest asked to predict at a covariate combination far outside its training data returns a confident value — the mean of whichever leaves it lands in — with no signal that it is guessing. Since project areas frequently include terrain, land-use, or climate combinations that no core sampled, a meaningful fraction of the map may be extrapolated. The dissimilarity index below makes that fraction explicit, and the honest treatment is to mask it or to widen the interval there rather than to present it as prediction.

Diagnostic Pipeline / Pre-Flight Validation

The pre-flight checks covariate coherence, computes the block size from the data rather than assuming it, and measures how much of the prediction area is extrapolated.

from dataclasses import dataclass

import numpy as np
import structlog
from scipy.spatial import cKDTree

log = structlog.get_logger()

MAX_EXTRAPOLATION_FRACTION = 0.25   # above this, the map is mostly guesswork


@dataclass(frozen=True)
class StackHealth:
    n_cores: int
    n_covariates: int
    imagery_window_days: int
    sampling_span_days: int
    block_size_m: float
    extrapolated_fraction: float
    usable: bool
    reason: str | None


def empirical_block_size(xy: np.ndarray, residuals: np.ndarray,
                         max_lag_m: float = 20_000.0, bins: int = 20) -> float:
    """Block size from the residual semivariogram range, not a guess.

    Blocks smaller than the autocorrelation range still leak a neighbour's answer
    across the fold boundary, which is the whole failure blocking exists to stop.
    """
    tree = cKDTree(xy)
    pairs = tree.query_pairs(max_lag_m, output_type="ndarray")
    if len(pairs) < 100:
        log.warning("dsm.variogram.sparse", pairs=len(pairs))
        return max_lag_m / 4.0

    d = np.linalg.norm(xy[pairs[:, 0]] - xy[pairs[:, 1]], axis=1)
    gamma = 0.5 * (residuals[pairs[:, 0]] - residuals[pairs[:, 1]]) ** 2

    edges = np.linspace(0, max_lag_m, bins + 1)
    idx = np.digitize(d, edges) - 1
    means = np.array([gamma[idx == b].mean() if (idx == b).any() else np.nan
                      for b in range(bins)])
    sill = np.nanmean(means[-5:])
    reached = np.where(means >= 0.95 * sill)[0]
    rng = float(edges[reached[0] + 1]) if reached.size else max_lag_m / 4.0

    log.info("dsm.variogram", range_m=round(rng, 1), sill=round(float(sill), 4),
             pairs=len(pairs))
    return rng


def dissimilarity_index(train: np.ndarray, predict: np.ndarray) -> np.ndarray:
    """Per-pixel distance to the nearest training point in scaled covariate space,
    normalised by the mean nearest-neighbour distance WITHIN the training set.

    Values above 1 mean the pixel is further from any core than cores typically are
    from each other — the model is extrapolating, however confident it sounds.
    """
    mu, sd = train.mean(axis=0), train.std(axis=0) + 1e-9
    train_z, predict_z = (train - mu) / sd, (predict - mu) / sd

    tree = cKDTree(train_z)
    within, _ = tree.query(train_z, k=2)           # k=2: skip the point itself
    scale = float(within[:, 1].mean())

    distance, _ = tree.query(predict_z, k=1)
    return distance / max(scale, 1e-9)


def preflight(cores_xy: np.ndarray, core_values: np.ndarray, train_cov: np.ndarray,
              predict_cov: np.ndarray, imagery_window_days: int,
              sampling_span_days: int) -> StackHealth:
    block = empirical_block_size(cores_xy, core_values - core_values.mean())
    di = dissimilarity_index(train_cov, predict_cov)
    extrapolated = float((di > 1.0).mean())

    reason = None
    if imagery_window_days > sampling_span_days * 3:
        # A wide imagery window describes a different landscape from the one sampled.
        reason = "imagery_window_incoherent_with_sampling"
    elif extrapolated > MAX_EXTRAPOLATION_FRACTION:
        reason = "excessive_extrapolation"

    health = StackHealth(
        n_cores=len(cores_xy), n_covariates=train_cov.shape[1],
        imagery_window_days=imagery_window_days, sampling_span_days=sampling_span_days,
        block_size_m=round(block, 1), extrapolated_fraction=round(extrapolated, 3),
        usable=reason is None, reason=reason,
    )
    log.info("dsm.preflight", **health.__dict__)
    return health
Where a soil model is interpolating and where it is guessing A scatter of training cores in a two-dimensional covariate space defined by slope on the horizontal axis and mean annual rainfall on the vertical axis. The cores occupy a diagonal band. A shaded region marks the area within one nearest-neighbour distance of a core, labelled interpolation. Two clusters of prediction pixels lie outside it: a steep, dry corner and a flat, very wet corner, both labelled extrapolation. A side panel gives the arithmetic: 78 percent of the project area is inside the covered region with a prediction interval of plus or minus 6.4 tonnes of carbon dioxide equivalent per hectare, while 22 percent is outside with an interval of plus or minus 19.1, and notes that the accuracy statistic quoted for the map as a whole applies only to the first group. The map's stated accuracy applies to part of the map Covariate space, two of nine dimensions shown. Cores define where prediction is interpolation. slope (°) → rainfall (mm yr⁻¹) → interpolation within one core-spacing steep & dry flat & very wet Project area, split honestly inside covered region 78% prediction interval ±6.4 tCO₂e ha⁻¹ outside covered region 22% prediction interval ±19.1 tCO₂e ha⁻¹ The headline accuracy describes the first row only.

Deterministic Transformation Logic

A quantile regression forest is the workhorse here because it returns the full conditional distribution rather than a point estimate, which gives per-pixel prediction intervals at no extra fitting cost. The implementation tunes with blocked folds, calibrates the interval against held-out data, and emits three surfaces.

import numpy as np
import structlog
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import GroupKFold, GridSearchCV

log = structlog.get_logger()


class QuantileForest(RandomForestRegressor):
    """Random forest that retains leaf memberships so arbitrary quantiles of the
    conditional distribution can be read off after fitting."""

    def fit(self, X, y, **kwargs):
        super().fit(X, y, **kwargs)
        self._y = np.asarray(y)
        self._leaf_index = [
            {leaf: np.where(tree.apply(X) == leaf)[0] for leaf in np.unique(tree.apply(X))}
            for tree in self.estimators_
        ]
        return self

    def predict_quantiles(self, X, quantiles=(0.05, 0.5, 0.95)) -> np.ndarray:
        out = np.empty((X.shape[0], len(quantiles)), dtype="float64")
        for i in range(X.shape[0]):
            pooled = []
            for tree, index in zip(self.estimators_, self._leaf_index):
                leaf = tree.apply(X[i:i + 1])[0]
                pooled.append(self._y[index[leaf]])
            pooled = np.concatenate(pooled)
            out[i] = np.quantile(pooled, quantiles)
        return out


def tune_blocked(X: np.ndarray, y: np.ndarray, groups: np.ndarray) -> dict:
    """Tune with GroupKFold over spatial blocks.

    Tuning with random folds selects the hyperparameters that memorise location
    best — usually deep trees with tiny leaves — and those are exactly the wrong
    ones for prediction at unsampled sites.
    """
    grid = {"n_estimators": [500], "min_samples_leaf": [2, 4, 8, 16],
            "max_features": ["sqrt", 0.3, 0.5]}
    search = GridSearchCV(
        RandomForestRegressor(random_state=0), grid,
        cv=GroupKFold(n_splits=min(5, len(set(groups)))), scoring="neg_root_mean_squared_error",
        n_jobs=-1,
    )
    search.fit(X, y, groups=groups)
    log.info("dsm.tuning", best=search.best_params_,
             blocked_rmse=round(-float(search.best_score_), 3))
    return search.best_params_


def calibrate_interval(model: QuantileForest, X_hold: np.ndarray, y_hold: np.ndarray,
                       target: float = 0.90) -> float:
    """Empirical coverage of the nominal interval, and the factor that fixes it.

    A quantile forest's nominal 90% interval routinely covers 70-80% on held-out
    soil data. Reporting the nominal figure without checking is the quiet way an
    uncertainty statement becomes untrue.
    """
    q = model.predict_quantiles(X_hold, quantiles=(0.05, 0.5, 0.95))
    covered = float(((y_hold >= q[:, 0]) & (y_hold <= q[:, 2])).mean())

    # Widen (or narrow) symmetrically about the median until coverage matches.
    factor = 1.0
    for _ in range(40):
        low = q[:, 1] - (q[:, 1] - q[:, 0]) * factor
        high = q[:, 1] + (q[:, 2] - q[:, 1]) * factor
        got = float(((y_hold >= low) & (y_hold <= high)).mean())
        if abs(got - target) < 0.01:
            break
        factor *= 1.05 if got < target else 0.97

    log.info("dsm.interval.calibration", nominal=target, raw_coverage=round(covered, 3),
             widening_factor=round(factor, 3))
    return factor


def fit_and_predict(
    X_train: np.ndarray, y_train: np.ndarray, groups: np.ndarray,
    X_hold: np.ndarray, y_hold: np.ndarray, X_predict: np.ndarray,
    dissimilarity: np.ndarray,
) -> dict:
    params = tune_blocked(X_train, y_train, groups)
    model = QuantileForest(random_state=0, **params).fit(X_train, y_train)
    factor = calibrate_interval(model, X_hold, y_hold)

    q = model.predict_quantiles(X_predict)
    median = q[:, 1]
    low = median - (median - q[:, 0]) * factor
    high = median + (q[:, 2] - median) * factor

    extrapolated = dissimilarity > 1.0
    # Outside the covered region the interval is not merely wide, it is unverified.
    # Mark it rather than pretending the calibration transfers.
    low[extrapolated] = np.nan
    high[extrapolated] = np.nan

    log.info("dsm.predict", pixels=len(median),
             extrapolated_fraction=round(float(extrapolated.mean()), 3),
             median_mean=round(float(np.nanmean(median)), 3),
             mean_interval_width=round(float(np.nanmean(high - low)), 3))
    return {"median": median, "low": low, "high": high,
            "extrapolated": extrapolated, "widening_factor": factor,
            "hyperparameters": params}

The interval calibration step is the one most pipelines omit and the one that most changes the reported figure. A quantile forest’s nominal 90% interval commonly achieves 70–80% empirical coverage on soil data, because the forest underestimates variance at sparsely sampled covariate combinations. Measuring coverage on held-out data and widening until it matches turns an interval that sounds rigorous into one that is.

Nominal against achieved coverage for an uncalibrated and a calibrated prediction interval A calibration plot with nominal interval coverage from 0.5 to 0.95 on the horizontal axis and achieved coverage measured on held-out cores on the vertical axis. A diagonal line marks perfect calibration. The uncalibrated quantile forest sits well below the diagonal throughout, achieving 0.38 at a nominal 0.50 and only 0.76 at a nominal 0.95. After widening by a factor of 1.34 the calibrated curve tracks the diagonal within two percentage points across the range. An annotation states that the uncalibrated 90 percent interval delivers 72 percent coverage, so roughly one core in four falls outside an interval described to the reader as containing nine in ten. An interval is a claim — measure whether it is true Achieved coverage on 84 held-out cores that took no part in fitting or tuning. perfect 1.0 0.8 0.6 0.4 0.50 0.65 0.80 0.95 nominal coverage achieved coverage uncalibrated calibrated ×1.34 At the nominal 90% interval 72% achieved roughly one core in four falls outside an interval described as holding nine in ten

Compliance Gating & Audit Trail Generation

Four artefacts make the map auditable. The blocked cross-validation score with its block size and the rationale for that size, alongside the random-fold score so the optimism gap is visible rather than hidden. The calibrated coverage of the reported interval, measured on held-out data. The extrapolation mask and the fraction of the project area it covers. And the covariate manifest: every layer, its source, its temporal window, and its resolution, versioned so a re-run resolves the same inputs.

Where the extrapolated fraction is material, the honest treatments are, in order of preference: collect additional cores targeting the uncovered covariate space — a small, well-targeted campaign often removes most of the extrapolation; restrict the crediting area to the covered region; or apply a conservativeness deduction sized to the wider interval. What is not defensible is reporting one accuracy figure over a map of which a fifth is extrapolated.

The stock conversion and reporting basis are handled separately, and the equivalent-soil-mass requirement described in the parent topic applies to whatever this model predicts. Route the outputs through the schema contract in the MRV data schema reference, and chain the covariate manifest into MRV data lineage and provenance tracking so a re-run five years later resolves the same imagery.

Production Integration

  1. Build the covariate stack on the canonical equal-area grid, windowing every imagery-derived layer to the sampling campaign and recording the window.
  2. Join cores to covariates with an explicit CRS on both sides, rejecting any core lacking a recorded laboratory method.
  3. Pre-flight: derive the block size from the residual variogram, compute the dissimilarity index, and fail on incoherent windows or excessive extrapolation.
  4. Tune and fit with blocked folds, then calibrate the interval against a held-out set that took no part in tuning.
  5. Predict three surfaces — median, calibrated interval, extrapolation mask — and mask the interval where the model extrapolates.
  6. Emit the manifest with both cross-validation scores, the block size, coverage, and the extrapolated fraction, and hand off to the design-based validation step.

At scale the quantile prediction dominates runtime because it pools leaf memberships per pixel. Chunk the prediction grid and process chunks in parallel with the same fitted model, and precompute the leaf index once rather than per chunk — the same tile-partitioned pattern used for async satellite tile processing with Dask.

Frequently Asked Questions

Which covariates actually matter for soil carbon?

Terrain derivatives — elevation, slope, curvature, and a wetness index — are consistently the strongest, because they control water movement and deposition. Climate normals matter at regional scale but are nearly constant within a single project and can then contribute little. Land-use history is often the single most informative layer where it exists at adequate quality, since management dominates soil carbon at field scale. Bare-soil spectral composites help on exposed soil and contribute nothing under permanent cover. Start with terrain plus land-use history, add the rest, and let the blocked score decide — a covariate that improves the random score but not the blocked score is adding location, not information.

How many cores do I need for a usable model?

Fewer than most people expect for a rough map, far more than most expect for a defensible one. A hundred well-spread cores can produce a model with genuine skill; the constraint is usually coverage of covariate space rather than raw count. Judge sufficiency by the extrapolated fraction rather than by a number: if 20% of your project sits outside the covered region, twenty well-targeted cores in that region are worth more than two hundred more in the middle of the existing cluster.

Is kriging or machine learning better here?

They answer slightly different questions and the modern default combines them. Kriging exploits spatial autocorrelation directly and gives principled uncertainty, but it struggles with many covariates. Tree ensembles exploit covariates well and handle non-linearity, but ignore residual spatial structure. Regression kriging — a machine-learning trend plus kriged residuals — usually beats either alone where the residuals still carry spatial structure. Check for that structure explicitly with a residual variogram; if the residuals are spatially uncorrelated, the tree model has already captured what there was.

Why calibrate the prediction interval instead of trusting the quantiles?

Because the quantiles are computed from the training distribution within leaves, which understates variance where the training data is sparse — exactly where the interval matters most. Empirical coverage on held-out data is the only way to know what your nominal 90% interval actually delivers. Report both the nominal and the achieved coverage, and the widening factor you applied. Verifiers respond well to this because it is visibly a check rather than a claim.

How should the model handle depth?

Two approaches work and they are not equivalent. The simplest is to model each standard depth increment separately — 0–5, 5–15, 15–30 cm — which is easy to implement and lets each increment have its own covariate relationships, since surface carbon responds to management while deeper carbon responds to texture and drainage. The alternative is to fit a depth function, typically a spline, to each profile and then model the spline coefficients spatially, which enforces a physically sensible continuity with depth and handles cores sampled at inconsistent increments. Depth functions are preferable when your cores come from mixed campaigns with different sampling protocols, which is the common case for projects that inherit legacy data. Whichever you choose, model concentration by increment and convert to stock afterwards, so the equivalent-soil-mass correction remains available.

What should I do when the model performs poorly no matter what I try?

Accept it and change the claim, rather than searching for a configuration that scores well by accident. A blocked R² near 0.2 usually means one of three things: the covariates genuinely do not carry the signal at this scale, which is common in flat, uniform landscapes where terrain explains nothing; the sample is too small or too clustered to support any model; or the dominant driver is management history you do not have. All three are diagnosable. If terrain and climate are near-constant across the area, look for management data. If the effective sample size is small, target new cores at the covariate gaps rather than adding more of the same. And if the signal simply is not there, a design-based estimate of the area mean — no map at all — is a legitimate and defensible output that many methodologies accept, and it is far better than a map with an invented accuracy.

Should I predict the whole project area or only the covered region?

Predict everywhere, then mask. The full surface is useful for planning and for identifying where to sample next, while the mask is what governs what may be credited. Shipping only the covered region loses information that helps target the next campaign; shipping the full surface without the mask invites a downstream consumer to credit the extrapolated part. Both layers, clearly labelled, is the arrangement that survives review.