Geostatistical Modeling Spatial Statistics with Python

Ordinary & Universal Kriging in Python

Ordinary and Universal Kriging are the two most widely deployed geostatistical interpolators in production spatial-analysis pipelines, sitting within the broader Kriging, Interpolation & Surface Generation Techniques family. Ordinary Kriging (OK) handles stationary fields where the mean is constant but unknown; Universal Kriging (UK) extends this to fields that carry a deterministic spatial trend. Both methods return not just a prediction surface but also a spatially explicit variance surface — the property that separates them from purely deterministic approaches like Inverse Distance Weighting.

← Back to Kriging, Interpolation & Surface Generation Techniques

Prerequisites

  • Python ≥ 3.10
  • numpy >= 1.26, scipy >= 1.12
  • geopandas >= 0.14, pyproj >= 3.6
  • pykrige >= 1.7
  • scikit-gstat >= 1.0
  • Input data in a projected CRS (e.g., UTM, State Plane) — geographic coordinates (EPSG:4326) will produce incorrect distance calculations
  • No duplicate coordinate pairs; no NaN values in the target variable
  • Minimum ~30 sample points for reliable variogram estimation (80–200+ recommended for production)
python
# pip install "numpy>=1.26" "scipy>=1.12" "geopandas>=0.14" "pykrige>=1.7" "scikit-gstat>=1.0"

Mathematical Core

Ordinary Kriging

The OK predictor at an unsampled location s0\mathbf{s}_0 is a weighted sum of nn observations:

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

The weights λ\boldsymbol{\lambda} are found by solving the system

[C11⊤0][λμ]=[c01]\begin{bmatrix} \mathbf{C} & \mathbf{1} \\ \mathbf{1}^\top & 0 \end{bmatrix} \begin{bmatrix} \boldsymbol{\lambda} \\ \mu \end{bmatrix} = \begin{bmatrix} \mathbf{c}_0 \\ 1 \end{bmatrix}

where Cij=C(si,sj)\mathbf{C}_{ij} = C(\mathbf{s}_i, \mathbf{s}_j) is the n×nn \times n covariance matrix between observations, c0\mathbf{c}_0 is the nn-vector of covariances between each observation and the prediction location, μ\mu is a Lagrange multiplier enforcing the unbiasedness constraint ∑λi=1\sum \lambda_i = 1, and C(⋅)C(\cdot) is the theoretical covariance derived from the fitted variogram via C(h)=C(0)−γ(h)C(h) = C(0) - \gamma(h).

The kriging variance (minimised prediction error variance) is:

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

Universal Kriging

UK decomposes the random field as

Z(s)=∑k=0Kakfk(s)+ε(s)Z(\mathbf{s}) = \sum_{k=0}^{K} a_k f_k(\mathbf{s}) + \varepsilon(\mathbf{s})

where fk(s)f_k(\mathbf{s}) are known drift functions (e.g., f0=1f_0=1, f1=xf_1=x, f2=yf_2=y for a first-order linear trend), aka_k are unknown drift coefficients, and ε(s)\varepsilon(\mathbf{s}) is a zero-mean spatially correlated residual. The kriging system grows to accommodate the additional drift constraints:

[CFF⊤0][λμ]=[c0f0]\begin{bmatrix} \mathbf{C} & \mathbf{F} \\ \mathbf{F}^\top & \mathbf{0} \end{bmatrix} \begin{bmatrix} \boldsymbol{\lambda} \\ \boldsymbol{\mu} \end{bmatrix} = \begin{bmatrix} \mathbf{c}_0 \\ \mathbf{f}_0 \end{bmatrix}

where Fik=fk(si)\mathbf{F}_{ik} = f_k(\mathbf{s}_i) is the n×(K+1)n \times (K+1) drift matrix evaluated at the observation locations, and f0\mathbf{f}_0 is the drift vector at the prediction location. Checking for second-order stationarity in the residuals after drift removal confirms that UK is correctly specified.

Ordinary vs Universal Kriging system structure Diagram showing how the OK kriging matrix extends to the UK system by adding drift function columns and rows. Left panel: OK system (C | 1 / 1T | 0). Right panel: UK system (C | F / FT | 0). An arrow labelled "add drift functions" points from OK to UK. Ordinary Kriging C n × n covariance 1 1ᵀ 0 Unbiasedness: Σλᵢ = 1 add drift functions Universal Kriging C n × n covariance F n × (K+1) drift matrix Fᵀ 0 f₀=1, f₁=x, f₂=y (linear drift) new in UK Both systems solve for weights λ and Lagrange multipliers μ that minimise prediction variance. UK requires the drift matrix F to be full-rank; collinear drift terms will produce a singular system. Shaded blocks are the additions UK makes to the OK system.

Data Preparation

Before fitting any variogram, the input data must be clean and in a suitable coordinate system. Sampling bias mitigation at this stage — for example, applying declustering weights before variogram computation — directly affects the accuracy of the fitted model.

python
import numpy as np
import geopandas as gpd

def validate_and_prepare(
    gdf: gpd.GeoDataFrame,
    value_col: str,
    target_crs: str = "EPSG:32633"
) -> tuple[np.ndarray, np.ndarray]:
    """
    Project to metric CRS, remove duplicate geometries, strip NaNs,
    and return (coords, values) arrays ready for variogram fitting.

    Parameters
    ----------
    gdf        : GeoDataFrame with Point geometry and a numeric value column.
    value_col  : Name of the column to interpolate.
    target_crs : EPSG code for the projected output CRS (default: UTM zone 33N).

    Returns
    -------
    coords : (n, 2) float64 array of (x, y) coordinates in metres.
    values : (n,)   float64 array of observed values.
    """
    if gdf.crs is None:
        raise ValueError("Input GeoDataFrame must have a defined CRS.")
    if not gdf.crs.is_projected:
        gdf = gdf.to_crs(target_crs)

    # Drop NaNs first, then exact-geometry duplicates (keep last measurement)
    clean = (
        gdf.dropna(subset=[value_col])
           .drop_duplicates(subset=["geometry"], keep="last")
    )
    coords = np.column_stack((clean.geometry.x, clean.geometry.y))
    values = clean[value_col].values.astype(np.float64)
    return coords, values

Variogram Fitting

The empirical semivariogram is the backbone of kriging. Every parameter — nugget, sill, range, and model family — must be chosen deliberately rather than left to library defaults.

Parameter definitions:

  • Nugget (C0C_0): discontinuity at the origin; captures measurement error and micro-scale variability below the sampling resolution.
  • Sill (C0+CC_0 + C): total variance where the variogram levels off; approximately equal to the marginal variance of the process.
  • Range (aa): lag distance at which the sill is effectively reached; observations separated by more than the range are spatially uncorrelated.

The Cressie–Hawkins robust estimator is preferred over the classical Matheron estimator when the dataset contains outliers:

2γ^(h)=(1∣N(h)∣∑N(h)∣Z(si)−Z(sj)∣1/2)40.457+0.494/∣N(h)∣2\hat{\gamma}(h) = \frac{\left(\frac{1}{|N(h)|}\sum_{N(h)}|Z(\mathbf{s}_i)-Z(\mathbf{s}_j)|^{1/2}\right)^4}{0.457 + 0.494/|N(h)|}

python
from skgstat import Variogram

def fit_variogram(
    coords: np.ndarray,
    values: np.ndarray,
    model_type: str = "spherical",
    n_lags: int = 20,
    maxlag: float | None = None
) -> Variogram:
    """
    Fit a theoretical variogram to empirical semivariogram estimates.

    Parameters
    ----------
    coords     : (n, 2) coordinate array in projected metres.
    values     : (n,)   observed values array.
    model_type : Theoretical model — 'spherical', 'exponential', or 'gaussian'.
    n_lags     : Number of lag bins for the empirical variogram (15–30 is typical).
    maxlag     : Upper distance limit for variogram fitting; defaults to half the
                 diagonal extent of the domain if None.

    Notes
    -----
    Use model_type='exponential' for geological data with asymptotic approach
    to the sill; use 'gaussian' for smooth, continuously differentiable fields.
    Always inspect V.plot() before using fitted parameters in kriging.
    """
    kwargs = dict(estimator="cressie", n_lags=n_lags, fit_method="trf")
    if maxlag is not None:
        kwargs["maxlag"] = maxlag
    V = Variogram(coords, values, **kwargs)
    V.model = model_type
    return V


def extract_variogram_params(V: Variogram) -> dict:
    """Return nugget, sill, and range from a fitted skgstat Variogram."""
    return {
        "nugget": V.parameters[2],   # index 2 in (range, sill, nugget)
        "sill":   V.parameters[1],
        "range":  V.parameters[0],
    }

Annotated Implementation

Ordinary Kriging

OK is the correct default when exploratory analysis (scatter of values vs. x, vs. y) reveals no systematic trend and the stationarity tests pass.

python
from pykrige.ok import OrdinaryKriging

def run_ordinary_kriging(
    coords: np.ndarray,
    values: np.ndarray,
    grid_x: np.ndarray,
    grid_y: np.ndarray,
    variogram_params: dict,
    variogram_model: str = "spherical"
) -> tuple[np.ndarray, np.ndarray]:
    """
    Execute Ordinary Kriging on a regular prediction grid.

    Parameters
    ----------
    coords           : (n, 2) observation coordinates (metres).
    values           : (n,)   observed values.
    grid_x, grid_y   : 1-D arrays defining the prediction grid axes.
    variogram_params : Dict with keys 'nugget', 'sill', 'range' from fit_variogram().
    variogram_model  : Must match the model used during variogram fitting.

    Returns
    -------
    z  : (len(grid_y), len(grid_x)) masked array of kriging predictions.
    ss : (len(grid_y), len(grid_x)) masked array of kriging variance.

    Notes
    -----
    Setting exact_values=False adds a small nugget perturbation that prevents
    singular covariance matrices when observations are very closely spaced.
    """
    ok = OrdinaryKriging(
        coords[:, 0],
        coords[:, 1],
        values,
        variogram_model=variogram_model,
        # Pass fitted parameters directly — never rely on PyKrige's auto-fit
        # when you have a validated variogram from scikit-gstat.
        variogram_parameters={
            "nugget": variogram_params["nugget"],
            "psill":  variogram_params["sill"] - variogram_params["nugget"],
            "range":  variogram_params["range"],
        },
        nlags=20,
        verbose=False,
        enable_plotting=False,
        exact_values=True,
    )
    z, ss = ok.execute("grid", grid_x, grid_y)
    return z, ss

Universal Kriging with Linear Drift

When a first-order spatial trend is confirmed, supply drift_terms=["regional_linear"]. For external covariate drift (e.g., elevation as a predictor of precipitation), pass custom drift functions — but ensure the drift matrix F\mathbf{F} is full-rank before proceeding.

python
from pykrige.uk import UniversalKriging

def run_universal_kriging(
    coords: np.ndarray,
    values: np.ndarray,
    grid_x: np.ndarray,
    grid_y: np.ndarray,
    variogram_params: dict,
    variogram_model: str = "spherical",
    drift_terms: list[str] | None = None
) -> tuple[np.ndarray, np.ndarray]:
    """
    Execute Universal Kriging with polynomial or external drift.

    Parameters
    ----------
    drift_terms : List of PyKrige drift keywords.
                  'regional_linear'  — first-order (x, y) trend.
                  'point_log'        — radially symmetric log trend.
                  Pass None to default to first-order linear.

    Notes
    -----
    For external covariate drift supply a custom callable via
    specified_drift in the UniversalKriging constructor instead.
    The drift matrix must not contain collinear columns — check
    np.linalg.matrix_rank(F) == F.shape[1] before fitting.
    """
    if drift_terms is None:
        drift_terms = ["regional_linear"]

    uk = UniversalKriging(
        coords[:, 0],
        coords[:, 1],
        values,
        variogram_model=variogram_model,
        variogram_parameters={
            "nugget": variogram_params["nugget"],
            "psill":  variogram_params["sill"] - variogram_params["nugget"],
            "range":  variogram_params["range"],
        },
        drift_terms=drift_terms,
        verbose=False,
    )
    z, ss = uk.execute("grid", grid_x, grid_y)
    return z, ss

Building the Prediction Grid

python
def make_prediction_grid(
    coords: np.ndarray,
    resolution: float,
    buffer: float = 0.0
) -> tuple[np.ndarray, np.ndarray]:
    """
    Build a regular grid spanning the convex hull of the observation points.

    Parameters
    ----------
    resolution : Cell size in the same units as coords (usually metres).
    buffer     : Optional margin added around the bounding box.
    """
    x_min, y_min = coords.min(axis=0) - buffer
    x_max, y_max = coords.max(axis=0) + 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

Output Interpretation

Reading the Prediction Surface

z (returned from ok.execute or uk.execute) is a masked NumPy array. Values outside the search neighbourhood or flagged as extrapolation are automatically masked. Before downstream analysis:

python
# Convert masked array to a plain float array (fill masked cells with NaN)
z_array  = np.where(z.mask, np.nan, z.data)
ss_array = np.where(ss.mask, np.nan, ss.data)

Reading the Variance Surface

ss is the kriging variance σOK2(s0)\sigma^2_{OK}(\mathbf{s}_0), not the standard deviation. A 95% prediction interval at each grid cell is:

Z^(s0)±1.96σOK2(s0)\hat{Z}(\mathbf{s}_0) \pm 1.96 \sqrt{\sigma^2_{OK}(\mathbf{s}_0)}

What good looks like:

  • Variance is lowest at and near observation points (it equals the nugget at exact data locations when exact_values=True).
  • Variance rises smoothly toward data-sparse zones and domain boundaries.
  • The spatial pattern of variance mirrors the sampling density — dense clusters have low local variance; voids have high variance.

Warning signs:

  • Flat variance surface: variogram range is too large; the covariance matrix is effectively constant and weights are near-uniform.
  • Negative kriging variance: numerical precision issue; add a small nugget (0.01× sill minimum) or check for near-duplicate coordinates.
  • Extreme variance spikes at observation locations: exact_values=False may be needed, or duplicates were not removed.

Integrating this variance surface into your pipeline for stakeholder-facing outputs is covered in Uncertainty & Variance Mapping.

Leave-One-Out Cross-Validation

python
from sklearn.model_selection import LeaveOneOut
from sklearn.metrics import mean_squared_error
import numpy as np

def loocv_kriging(coords, values, variogram_params, model="spherical"):
    """
    Leave-one-out cross-validation for Ordinary Kriging.

    Returns RMSE, MAE, and mean standardised residual
    (should be ~0 with std ~1 for a well-calibrated variogram).
    """
    from pykrige.ok import OrdinaryKriging
    loo = LeaveOneOut()
    predictions, actuals = [], []

    for train_idx, test_idx in loo.split(coords):
        ok = OrdinaryKriging(
            coords[train_idx, 0], coords[train_idx, 1], values[train_idx],
            variogram_model=model,
            variogram_parameters={
                "nugget": variogram_params["nugget"],
                "psill":  variogram_params["sill"] - variogram_params["nugget"],
                "range":  variogram_params["range"],
            },
            verbose=False,
            enable_plotting=False,
        )
        z_pred, _ = ok.execute("points",
                               coords[test_idx, 0],
                               coords[test_idx, 1])
        predictions.append(float(z_pred))
        actuals.append(values[test_idx[0]])

    predictions = np.array(predictions)
    actuals     = np.array(actuals)
    residuals   = actuals - predictions
    rmse = np.sqrt(mean_squared_error(actuals, predictions))
    mae  = np.mean(np.abs(residuals))
    return {"rmse": rmse, "mae": mae, "mean_residual": residuals.mean()}

A spatial k-fold cross-validation approach is preferable to random LOO for datasets with strong spatial autocorrelation, as random splits underestimate true generalisation error.

Production Considerations

Computational Scaling

Kriging solves an (n+1)×(n+1)(n+1) \times (n+1) linear system (OK) or (n+K+1)×(n+K+1)(n + K + 1) \times (n + K + 1) system (UK) per prediction location when using a global neighbourhood. This scales as O(n3)\mathcal{O}(n^3) for factorisation. For n>500n > 500, always use a local search neighbourhood:

python
# Local neighbourhood: restrict to the 30 nearest neighbours
ok = OrdinaryKriging(
    ...,
    nlags=20,
)
# PyKrige searches within a neighbourhood automatically when you call
# ok.execute("grid", grid_x, grid_y, backend="C")
# The C backend is ~10x faster than the pure-Python backend for grids > 10^4 cells.

Memory Management for Large Grids

For grids exceeding 10610^6 cells, execute in spatial chunks to avoid MemoryError:

python
def chunked_kriging(ok_model, grid_x, grid_y, chunk_size=500):
    """Tile the y-axis into strips and concatenate results."""
    z_strips, ss_strips = [], []
    for start in range(0, len(grid_y), chunk_size):
        gy_chunk = grid_y[start : start + chunk_size]
        z_c, ss_c = ok_model.execute("grid", grid_x, gy_chunk)
        z_strips.append(z_c)
        ss_strips.append(ss_c)
    return np.ma.vstack(z_strips), np.ma.vstack(ss_strips)

Anisotropy

Geological and hydrological processes rarely exhibit isotropic spatial continuity. Detect anisotropy by computing directional variograms at 0°, 45°, 90°, and 135°; if the range differs significantly by azimuth, set:

python
ok = OrdinaryKriging(
    ...,
    anisotropy_scaling=2.5,   # major_range / minor_range
    anisotropy_angle=45.0,    # clockwise rotation from north (degrees)
)

Block Kriging

When the output must represent areal averages (census tracts, satellite pixel footprints), switch to block kriging by passing a pointcloud of sub-block integration points instead of a single prediction location. This matters when reporting regulatory compliance quantities — point kriging predictions are not equal to block averages in non-linear problems.

Troubleshooting

Symptom Likely cause Fix
LinAlgError: Singular matrix Duplicate coordinates, zero nugget Deduplicate with drop_duplicates(subset=["geometry"]); add nugget ≥ 0.01 × sill
Predictions far outside observed range Variogram range ≫ domain extent; weights nearly uniform Fit variogram with maxlag ≤ half the domain diagonal; verify V.plot()
MemoryError on large grids Full grid materialised in RAM Use chunked execution or backend="C" with a local neighbourhood
Kriging variance is negative Near-duplicate coordinates, precision error Enforce minimum nugget; use exact_values=False
UK drift matrix singular Collinear drift terms Check np.linalg.matrix_rank(F); remove redundant drift covariates
LOO RMSE much larger than variogram RMSE Variogram fitted on full dataset; overfitting Re-fit variogram on each training fold inside the CV loop
High standardised residuals at cluster boundaries Incorrect stationarity assumption Test for second-order stationarity and switch to UK if a trend is detected
Flat prediction surface with no local detail Search neighbourhood too large Reduce nlags and verify the variogram range is realistic for the domain

Next Steps

For a step-by-step walkthrough of setting up PyKrige, customising variogram models, and exporting prediction rasters to GeoTIFF, see the Step-by-Step Ordinary Kriging with PyKrige guide. Once your prediction and variance surfaces are ready, Uncertainty & Variance Mapping shows how to convert kriging variance into confidence intervals and risk layers for stakeholder reporting.

Related

Step-by-Step Ordinary Kriging with PyKrige Inverse Distance Weighting: Deterministic Spatial Interpolation in Python Uncertainty & Variance Mapping in Geostatistics