Uncertainty & Variance Mapping in Geostatistics

Kriging solves two problems simultaneously: it estimates the value at each unsampled location, and it computes the statistical confidence of that estimate. The second output — the kriging variance surface — is what separates geostatistical methods from deterministic approaches and what makes them indispensable for regulatory compliance, resource allocation, and any decision that depends on knowing how much to trust a map. This guide is part of the broader Kriging, Interpolation & Surface Generation Techniques framework and covers the full variance-mapping workflow: mathematical derivation, annotated Python implementation, output diagnostics, and production performance patterns.

Prerequisites #

• Python 3.9 or later
• numpy >= 1.24, geopandas >= 0.13, pykrige >= 1.7, scipy >= 1.11, matplotlib >= 3.7, scikit-learn >= 1.3, rasterio >= 1.3
• Input point dataset in a projected CRS (UTM, State Plane, or equivalent) — geographic coordinates (EPSG:4326) will produce variogram distances in decimal degrees, making range and sill values meaningless
• No duplicate coordinates; no NaN values in the target variable
• Familiarity with semivariogram structure (nugget, sill, range) — see Ordinary & Universal Kriging for a variogram primer

pip install "numpy>=1.24" "geopandas>=0.13" "pykrige>=1.7" \
            "scipy>=1.11" "matplotlib>=3.7" "scikit-learn>=1.3" \
            "rasterio>=1.3"

Mathematical Core #

The Kriging System of Equations #

Ordinary Kriging (OK) finds a set of weights $\lambda_i$ for $n$ neighboring samples such that the estimator

$$\hat{Z}(\mathbf{s}_0) = \sum_{i=1}^{n} \lambda_i \, Z(\mathbf{s}_i)$$

is unbiased ($\sum_i \lambda_i = 1$) and minimizes the estimation variance. Solving the augmented system

$$\begin{pmatrix} \mathbf{C} & \mathbf{1} \\ \mathbf{1}^T & 0 \end{pmatrix} \begin{pmatrix} \boldsymbol{\lambda} \\ \mu \end{pmatrix} = \begin{pmatrix} \mathbf{c}_0 \\ 1 \end{pmatrix}$$

where $\mathbf{C}_{ij} = C(\mathbf{s}_i - \mathbf{s}_j)$ is the covariance between samples, $\mathbf{c}_0$ is the vector of covariances between samples and the prediction point, and $\mu$ is the Lagrange multiplier, yields the optimal weights. The kriging variance follows directly from the solution:

$$\sigma^2_{OK}(\mathbf{s}_0) = C(\mathbf{0}) - \boldsymbol{\lambda}^T \mathbf{c}_0 - \mu$$

where $C(\mathbf{0})$ is the process variance (the sill). Two properties of $\sigma^2_{OK}$ are immediately useful in practice:

1. Variance is zero at observed locations when the nugget effect is absent — the prediction is exact.
2. Variance grows monotonically with distance from the nearest samples until it plateaus at the sill — beyond the variogram range, all new observations are equally (un)informative.

A third property is less obvious and catches most newcomers: $\sigma^2_{OK}$ contains no term in $Z$. The right-hand side is assembled entirely from $C(\mathbf{0})$, $\mathbf{c}_0$, $\boldsymbol{\lambda}$ and $\mu$, every one of which depends only on the sample geometry and the fitted variogram — never on the measured values. Multiply every observation by ten and the variance surface does not move by a single unit. The Lagrange multiplier $\mu$ is the price paid for the unbiasedness constraint: it is near zero when the prediction point sits inside a well-distributed cloud of samples and grows when the estimator has to extrapolate, which is why $\sigma^2_{OK}$ climbs past the sill outside the sample hull rather than stopping at it. Read correctly, then, a kriging variance map is a sampling-configuration map. It answers “how well is this cell surrounded by data?”, not “how variable is the phenomenon here?”. Two consequences follow. First, the surface can be computed before a single measurement exists, which is exactly what makes kriging usable for sampling-network design — score candidate station layouts by the variance they would produce and add stations where the reduction is largest. Second, a genuinely heteroscedastic process — assay grades that scatter more where the mean grade is high, say — is handed the same uncertainty in its quiet regions as in its violent ones. When that matters, the fix is not a better variogram but a variance-stabilising transform (log or normal-score) applied before kriging, with the confidence limits back-transformed rather than $\sigma^2$ itself.

Relationship Between Variogram and Variance #

The semivariogram $\gamma(\mathbf{h})$ and the covariance function $C(\mathbf{h})$ are linked by

$$\gamma(\mathbf{h}) = C(\mathbf{0}) - C(\mathbf{h})$$

so the covariance matrix $\mathbf{C}$ in the kriging system can be computed from the fitted variogram model. Choosing a variogram model (spherical, exponential, Gaussian, Matérn) directly controls how quickly variance rises away from sampled locations. The spherical model, for instance, reaches its sill abruptly at the range, while the exponential model approaches it asymptotically — the choice matters most in data-sparse regions at the edge of the domain.

Annotated Implementation #

Step 1 — Load and Project Input Data #

import logging
import geopandas as gpd
import numpy as np
from pathlib import Path

logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")

def prepare_spatial_data(
    input_path: Path,
    value_col: str,
    target_crs: int = 32618,
) -> tuple[np.ndarray, np.ndarray]:
    """
    Load point features, enforce projected CRS, remove duplicates,
    and return clean coordinate and value arrays.
    """
    gdf = gpd.read_file(input_path)

    if gdf.crs is None:
        raise ValueError("Dataset has no CRS. Assign one before proceeding.")

    if not gdf.crs.is_projected:
        logging.warning(
            "Geographic CRS detected — reprojecting to EPSG:%d.", target_crs
        )
        gdf = gdf.to_crs(target_crs)

    # Drop rows with missing target values before removing duplicates
    gdf = gdf.dropna(subset=[value_col])

    dupes = gdf.duplicated(subset=["geometry"], keep="first")
    if dupes.any():
        logging.info("Removed %d duplicate coordinate(s).", dupes.sum())
        gdf = gdf[~dupes].reset_index(drop=True)

    coords = np.column_stack([gdf.geometry.x, gdf.geometry.y])
    values = gdf[value_col].to_numpy(dtype=np.float64)
    return coords, values

Step 2 — Build the Prediction Grid #

def build_prediction_grid(
    coords: np.ndarray,
    resolution: float,
    buffer: float = 0.0,
) -> tuple[np.ndarray, np.ndarray]:
    """
    Create a regular grid spanning the sample extent plus an optional buffer.

    Parameters
    ----------
    resolution : float
        Cell size in the same linear units as coords (metres for UTM).
    buffer : float
        Extra padding around the sample bounding box (metres). Useful for
        reducing edge-artifact masking; set to at least the variogram range.
    """
    x_min, x_max = coords[:, 0].min() - buffer, coords[:, 0].max() + buffer
    y_min, y_max = coords[:, 1].min() - buffer, coords[:, 1].max() + buffer

    grid_x = np.arange(x_min, x_max + resolution, resolution)
    grid_y = np.arange(y_min, y_max + resolution, resolution)
    return grid_x, grid_y

Step 3 — Run Ordinary Kriging and Extract Variance #

from pykrige.ok import OrdinaryKriging

def run_kriging_with_variance(
    coords: np.ndarray,
    values: np.ndarray,
    grid_x: np.ndarray,
    grid_y: np.ndarray,
    variogram_model: str = "spherical",
    nlags: int = 12,
) -> tuple[np.ndarray, np.ndarray]:
    """
    Fit a variogram, execute Ordinary Kriging, and return both the
    prediction surface and the kriging variance surface.

    Returns
    -------
    z : np.ndarray, shape (len(grid_y), len(grid_x))
        Predicted values at each grid node.
    sigma2 : np.ndarray, same shape
        Kriging variance (σ²) at each grid node — NOT standard error.
        Take np.sqrt(sigma2) to get σ in the same units as `values`.
    """
    x, y = coords[:, 0], coords[:, 1]

    krig = OrdinaryKriging(
        x, y, values,
        variogram_model=variogram_model,
        nlags=nlags,            # number of lag bins for empirical variogram
        weight=True,            # weight bins by pair count — reduces outlier influence
        enable_plotting=False,
        verbose=False,
    )

    # execute() returns (prediction, variance) — both are masked arrays
    z, sigma2 = krig.execute("grid", grid_x, grid_y)

    # Log fitted variogram parameters for reproducibility
    params = krig.variogram_model_parameters
    logging.info(
        "Variogram fitted | model=%s nugget=%.4f sill=%.4f range=%.1f",
        variogram_model, params[2], params[0] + params[2], params[1],
    )

    return np.ma.filled(z, np.nan), np.ma.filled(sigma2, np.nan)

Key decisions in this implementation:

• weight=True down-weights variogram lag bins with few pairs, improving stability with irregular or clustered sample layouts.
• sigma2 from execute() is the variance (σ²). Squaring it again is a common off-by-one error; take np.sqrt(sigma2) for standard error.
• Masked array fill to np.nan propagates naturally through NumPy operations and Matplotlib’s imshow.

Variogram Configuration and Diagnostics #

Variance surfaces are only as credible as the variogram model they derive from. Three parameters govern the spatial decay of confidence:

from skgstat import Variogram

def diagnose_variogram(coords: np.ndarray, values: np.ndarray) -> None:
    """
    Fit and print variogram diagnostics using scikit-gstat.
    Use these parameter estimates to inform pykrige model selection.
    """
    v = Variogram(
        coords,
        values,
        model="spherical",
        n_lags=12,
        maxlag=0.6,          # use at most 60 % of the max pairwise distance
        fit_method="trf",    # Trust Region Reflective — robust to outliers
    )
    print(v.describe())
    # Inspect the nugget-to-sill ratio:
    nsr = v.parameters[2] / (v.parameters[0] + v.parameters[2])
    if nsr > 0.7:
        logging.warning(
            "Nugget-to-sill ratio %.2f > 0.7: most variation is unexplained "
            "noise. Variance will remain high across the domain.", nsr
        )

A well-fitted variogram shows:

• Empirical bins that track the theoretical curve within reasonable scatter.
• A nugget that accounts for genuine measurement noise, not half the sill.
• A range that reflects a physically plausible correlation distance for the phenomenon.

One caveat sits underneath all of this. The kriging system treats the variogram as known exactly, but $C_0$, $C$ and $a$ are estimates fitted to a scatter of empirical lag points, and their own uncertainty never enters $\sigma^2_{OK}$. The variance surface is therefore systematically optimistic, and the shortfall widens as the fit thins out — barely noticeable with several hundred well-spread samples, material below roughly fifty, where two or three influential lag bins can move the fitted range by a factor of two. The effect compounds with weak stationarity: an unmodelled trend inflates the empirical variogram at long lags, so the fitted sill comes out too high and the range too short, and the resulting surface manages to overstate uncertainty in the sampled interior while understating the error made when predicting past the sample hull. That combination is the worst case, because the map looks conservative exactly where it is not. Cross-validation, covered below, is the practical remedy: its z-scores are measured against residuals the model actually produced rather than against the ones its assumptions predict.

Output Interpretation #

Reading the Variance Surface: What Good Looks Like #

A well-behaved variance surface has three recognizable features:

1. Low-variance halos around sample clusters. The minimum variance occurs at sample locations (zero when nugget = 0) and grows outward. If samples are reasonably spaced, these halos merge into a low-variance interior separated from a high-variance fringe near the grid boundary.
2. Variance plateau at the sill beyond the range. Cells more than one variogram range from any sample will all show variance values near the sill — the kriging estimator is effectively guessing the population mean there.
3. No isolated spikes inside the domain. Isolated high-variance cells embedded in a low-variance field signal either a duplicate removal failure (single nearby sample was silently dropped) or a numerical singularity in the covariance matrix inversion.

Warning Signs #

Production Considerations #

Memory and Scaling #

Ordinary Kriging constructs and inverts a dense $n \times n$ covariance matrix where $n$ is the number of sample points. This scales as $O(n^3)$ for factorisation and $O(n^2)$ for storage:

For large datasets, use one of three strategies:

• Moving-window kriging: partition the prediction grid into tiles and restrict the kriging search neighborhood to the nearest $k$ samples per tile. pykrige supports this via backend="loop" with a n_closest_points parameter.
• Sparse covariance approximations: gstools implements matrix-free kriging via FFT-based covariance computations, suitable for gridded or near-regular sample configurations.
• Block kriging: aggregate sample support areas into blocks; variance represents average uncertainty over each block rather than a point, substantially reducing matrix size.

Chunked Grid Processing #

def process_grid_chunked(
    coords: np.ndarray,
    values: np.ndarray,
    grid_x: np.ndarray,
    grid_y: np.ndarray,
    chunk_size: int = 400,
    variogram_model: str = "spherical",
) -> np.ndarray:
    """
    Tile the prediction grid to keep per-tile memory within bounds.
    Instantiate the kriging object once and reuse it across all tiles.
    """
    x, y = coords[:, 0], coords[:, 1]
    krig = OrdinaryKriging(
        x, y, values,
        variogram_model=variogram_model,
        enable_plotting=False,
        verbose=False,
    )

    ny, nx = len(grid_y), len(grid_x)
    full_variance = np.full((ny, nx), np.nan, dtype=np.float64)

    for i in range(0, ny, chunk_size):
        for j in range(0, nx, chunk_size):
            sub_x = grid_x[j : j + chunk_size]
            sub_y = grid_y[i : i + chunk_size]
            _, sub_var = krig.execute("grid", sub_x, sub_y)
            full_variance[i : i + chunk_size, j : j + chunk_size] = (
                np.ma.filled(sub_var, np.nan)
            )
            logging.debug("Completed tile (%d:%d, %d:%d).", i, i+chunk_size, j, j+chunk_size)

    return full_variance

Writing Variance to GeoTIFF #

import rasterio
from rasterio.transform import from_bounds

def export_variance_geotiff(
    variance: np.ndarray,
    grid_x: np.ndarray,
    grid_y: np.ndarray,
    crs_epsg: int,
    output_path: str,
) -> None:
    """Export the variance surface as a single-band GeoTIFF."""
    transform = from_bounds(
        grid_x.min(), grid_y.min(), grid_x.max(), grid_y.max(),
        width=len(grid_x), height=len(grid_y),
    )
    with rasterio.open(
        output_path, "w",
        driver="GTiff",
        height=variance.shape[0],
        width=variance.shape[1],
        count=1,
        dtype=variance.dtype,
        crs=f"EPSG:{crs_epsg}",
        transform=transform,
        nodata=np.nan,
        compress="lzw",
    ) as dst:
        dst.write(variance, 1)
    logging.info("Variance raster written to %s", output_path)

Troubleshooting #

Cross-Validation of Variance Estimates #

A kriging variance surface is a model prediction, not an empirical measurement. Validate it against held-out data using leave-one-out or spatial k-fold cross-validation. Spatial cross-validation from the cross-validation strategies framework is critical here — standard random k-fold underestimates error because nearby samples are spatially autocorrelated.

from sklearn.model_selection import LeaveOneOut
import numpy as np

def loocv_variance_calibration(
    coords: np.ndarray,
    values: np.ndarray,
    variogram_model: str = "spherical",
) -> dict[str, float]:
    """
    Leave-one-out cross-validation to assess variance calibration.
    A well-calibrated model has mean(z_score^2) ≈ 1.0.
    """
    loo = LeaveOneOut()
    z_scores = []

    for train_idx, test_idx in loo.split(coords):
        krig = OrdinaryKriging(
            coords[train_idx, 0],
            coords[train_idx, 1],
            values[train_idx],
            variogram_model=variogram_model,
            enable_plotting=False,
            verbose=False,
        )
        z_pred, s2_pred = krig.execute(
            "points",
            coords[test_idx, 0],
            coords[test_idx, 1],
        )
        z_actual = values[test_idx]
        sigma = np.sqrt(np.ma.filled(s2_pred, np.nan))
        if sigma > 0:
            z_scores.append(float((z_actual - z_pred) / sigma))

    z_scores = np.array(z_scores)
    return {
        "mean_z_score": float(np.mean(z_scores)),
        "mean_squared_z_score": float(np.mean(z_scores**2)),
        "rmse": float(np.sqrt(np.mean((z_scores * np.std(values))**2))),
    }

Interpretation: mean_squared_z_score near 1.0 indicates well-calibrated variance — the model neither overstates nor understates prediction uncertainty. Values < 1.0 mean the variance is overestimated (conservative); values > 1.0 indicate underestimated variance (dangerous for risk-sensitive decisions).

Next Steps #

For step-by-step variogram fitting and parameter selection within the pykrige and scikit-gstat APIs, see Step-by-Step Ordinary Kriging with PyKrige. If your domain exhibits a systematic spatial trend that inflates kriging variance near the domain boundary, switch to Universal Kriging as covered in Ordinary & Universal Kriging before running variance extraction. Where second-order stationarity cannot be verified, revisit your variogram fitting strategy before relying on variance outputs for production decisions.

Related

• Ordinary & Universal Kriging — mathematical foundations and pykrige implementation patterns
• Step-by-Step Ordinary Kriging with PyKrige — end-to-end variogram fitting and grid execution
• Inverse Distance Weighting — deterministic baseline to compare against kriging variance outputs
• Cross-Validation Strategies — spatial k-fold and buffered leave-one-out for validating interpolation models
• Stationarity & Trend Analysis — pre-requisite checks before trusting kriging variance estimates

← Back to Kriging, Interpolation & Surface Generation Techniques