Geostatistical Modeling Spatial Statistics with Python

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
bash
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 λi\lambda_i for nn neighboring samples such that the estimator

Z^(s0)=∑i=1nλi Z(si)\hat{Z}(\mathbf{s}_0) = \sum_{i=1}^{n} \lambda_i \, Z(\mathbf{s}_i)

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

(C11T0)(λμ)=(c01)\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 Cij=C(si−sj)\mathbf{C}_{ij} = C(\mathbf{s}_i - \mathbf{s}_j) is the covariance between samples, c0\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:

σOK2(s0)=C(0)−λTc0−μ\sigma^2_{OK}(\mathbf{s}_0) = C(\mathbf{0}) - \boldsymbol{\lambda}^T \mathbf{c}_0 - \mu

where C(0)C(\mathbf{0}) is the process variance (the sill). Two properties of σOK2\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.

Relationship Between Variogram and Variance

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

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

so the covariance matrix C\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.

Kriging Variance and Sample Geometry Three zones illustrated left to right: dense sampling with low variance (blue shading), moderate sampling with medium variance (yellow), and sparse/edge zone with high variance (orange). Arrows annotate the variance level in each zone and a curve above shows variance rising with distance from nearest sample. Dense sampling Moderate sampling Sparse / edge zone σ² ≈ 0 – 0.2 × sill σ² ≈ 0.4 – 0.7 × sill σ² → sill σ²(s) Distance from nearest sample →

Annotated Implementation

Step 1 — Load and Project Input Data

python
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

python
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

python
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:

Parameter Definition Effect on variance surface
Nugget (C0C_0) Discontinuity at lag = 0; captures micro-scale noise and positioning error High nugget → variance stays elevated even adjacent to samples
Sill (C0+CC_0 + C) Total variance at which the variogram plateaus; equals C(0)C(\mathbf{0}) Sets the upper bound on kriging variance across the domain
Range (aa) Distance at which the variogram first reaches the sill Controls how far “influence” spreads; beyond the range, variance equals the sill
python
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.

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

Observation Likely cause
Variance uniformly high across the entire domain Variogram range fitted far too small; all cells treat every sample as beyond-range
Variance = 0 everywhere Variogram nugget forced to sill; no spatial structure fitted
Circular bullseye halos with sharp edges Grid resolution much finer than variogram range; appears correct but wastes compute
Variance increases toward known sample clusters Anisotropy unmodeled; the fitted isotropic range is perpendicular to the main correlation direction

Production Considerations

Memory and Scaling

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

Sample count Covariance matrix (float64) Factorisation cost
500 2 MB milliseconds
2,000 32 MB ~1 s
10,000 800 MB minutes
50,000 20 GB impractical

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 kk 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

python
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

python
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

Symptom Likely cause Fix
LinAlgError: singular matrix during kriging Duplicate or near-duplicate sample coordinates Enforce minimum separation distance; jitter by 1e-6 map units if data source forces duplicates
Variance surface is all NaN Grid falls outside sample extent; masked array not filled Verify grid bounds overlap sample bounding box; call np.ma.filled(ss, np.nan) explicitly
Variance near zero everywhere Variogram sill fitted near zero; likely wrong units Check that coordinates and values share compatible units; inspect krig.variogram_model_parameters
High nugget-to-sill ratio (>0.7) Measurement noise dominates; or positional error in GPS coordinates Review data collection methodology; consider nugget fixing via variogram_parameters
pykrige fitting silently falls back to linear model Insufficient sample count for nlags bins Reduce nlags or increase sample size; verify at least 30+ pairs per lag bin
Variance spikes at grid edges Asymmetric search neighborhood near boundary Add a buffer equal to the variogram range; mask the edge strip before reporting
Memory error on large grid n×nn \times n covariance matrix too large Switch to chunked execution with n_closest_points; consider gstools FFT backend
Anisotropic variance pattern Isotropic variogram imposed on directionally structured data Fit geometric anisotropy (range ratio + rotation angle) in skgstat or pykrige

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.

python
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

← Back to Kriging, Interpolation & Surface Generation Techniques

Related

Inverse Distance Weighting: Deterministic Spatial Interpolation in Python Ordinary & Universal Kriging in Python