Step-by-Step Ordinary Kriging with PyKrige

TL;DR #

Instantiate OrdinaryKriging(x, y, z, variogram_model='spherical', variogram_parameters=[sill, range, nugget]), then call z_pred, variance = ok.execute('grid', gridx, gridy). Coordinates must be 1D NumPy arrays in a projected (metric) CRS. The variogram_parameters list order is [sill, range, nugget]. Always check that z_pred contains no NaN values and inspect variance to identify data-sparse regions before using the surface downstream.

Why This Workflow Matters #

Ordinary kriging is the go-to spatial interpolator when the underlying field is stationary — meaning the mean is constant but unknown across the domain. It appears throughout environmental monitoring, soil science, air-quality mapping, and mining resource estimation because it is both statistically unbiased and delivers a spatially explicit uncertainty estimate alongside every prediction.

Within the broader Ordinary & Universal Kriging toolkit, ordinary kriging is the starting point: master it first, then extend to universal kriging when a measurable large-scale trend demands a drift term. The approach sits inside the larger Kriging, Interpolation & Surface Generation Techniques workflow where surface generation feeds directly into uncertainty and variance mapping downstream. Getting the variogram configuration and coordinate handling right here prevents compounding errors in every subsequent step.

Environment and Version Pinning #

pip install pykrige==1.7.2 numpy==1.26.4 scipy==1.13.0 matplotlib==3.9.0 pyproj==3.6.1

import numpy as np          # 1.26.4
from pykrige.ok import OrdinaryKriging  # pykrige 1.7.2
import matplotlib.pyplot as plt         # 3.9.0
from pyproj import Transformer          # pyproj 3.6.1

Python 3.10 or 3.11 is recommended. PyKrige 1.7.x requires NumPy ≥ 1.20 and SciPy ≥ 1.6. If your environment has NumPy 2.x, pin numpy<2 until upstream compatibility is confirmed.

Step-by-Step Implementation #

Step 1 — Prepare Coordinate and Value Arrays #

PyKrige requires strictly 1D NumPy arrays for all three inputs. Coordinates must be in a projected metric CRS — Euclidean distance is assumed throughout covariance matrix construction. If your source data is in geographic coordinates (EPSG:4326), reproject before passing to the kriging object.

from pyproj import Transformer
import numpy as np

# Example: reproject from WGS-84 lat/lon to UTM zone 32N (EPSG:32632)
transformer = Transformer.from_crs("EPSG:4326", "EPSG:32632", always_xy=True)

lon = np.array([8.45, 9.12, 10.03, 8.78, 9.55, 10.20, 8.90, 9.33])
lat = np.array([47.2,  47.8, 47.5,  48.1, 48.4, 48.0,  47.6, 48.2])
z   = np.array([12.4,  15.1, 13.8,  16.0, 14.5, 17.2,  13.1, 15.7])

x_m, y_m = transformer.transform(lon, lat)  # now in metres

Validate the arrays before proceeding:

assert x_m.ndim == 1 and y_m.ndim == 1 and z.ndim == 1, "All inputs must be 1D"
assert not np.any(np.isnan(z)), "z contains NaN — remove or impute before kriging"
assert len(x_m) == len(y_m) == len(z), "Array length mismatch"
print(f"n={len(z)}  z: min={z.min():.2f}, max={z.max():.2f}, mean={z.mean():.2f}")

Step 2 — Define the Target Grid #

The grid must cover the convex hull of your sample points. Extending significantly beyond that hull produces high-variance extrapolation zones. Use np.linspace for regular grids; match units to your projected CRS.

# Grid bounds match the sample extent in metres
gridx = np.linspace(x_m.min(), x_m.max(), 200)
gridy = np.linspace(y_m.min(), y_m.max(), 200)

Resolution trades off computation time against surface detail. For exact kriging, interpolation cost scales with the number of prediction points, not grid resolution itself, but denser grids require more memory for the output arrays.

Step 3 — Choose and Configure the Variogram Model #

The variogram defines the spatial autocorrelation structure and controls prediction accuracy more than almost any other parameter. Before supplying variogram_parameters, inspect the empirical semivariogram to choose a realistic range estimate.

from pykrige.ok import OrdinaryKriging

# Spherical model: realistic plateau cutoff at the range distance.
# variogram_parameters = [sill, range, nugget]
# - sill  : total variance at the plateau (~sample variance as a starting point)
# - range : distance at which sill is reached (inspect your empirical variogram)
# - nugget: y-intercept — micro-scale variability or measurement noise

ok = OrdinaryKriging(
    x_m, y_m, z,
    variogram_model='spherical',
    variogram_parameters=[4.5, 45000.0, 0.2],  # units: variance, metres, variance
    nlags=15,           # number of lag bins for empirical variogram construction
    weight=True,        # weight lag bins by number of pairs (reduces outlier influence)
    verbose=False,
    coordinates_type='euclidean'
)

Inspect the fitted model before proceeding to interpolation:

# Access fitted parameters to verify reasonableness
params = ok.variogram_model_parameters
print(f"Fitted sill={params[0]:.3f}, range={params[1]:.1f}m, nugget={params[2]:.3f}")

# Nugget-to-sill ratio > 0.3 suggests strong measurement noise or very short-range structure
nugget_ratio = params[2] / params[0]
print(f"Nugget/sill ratio: {nugget_ratio:.2f}  {'⚠ high noise' if nugget_ratio > 0.3 else 'OK'}")

Step 4 — Execute Grid Interpolation #

# Returns: z_pred (predicted surface), ss (kriging variance)
z_pred, ss = ok.execute('grid', gridx, gridy)

print(f"Output shape : {z_pred.shape}")
print(f"Predicted    : min={z_pred.min():.3f}, max={z_pred.max():.3f}")
print(f"Variance     : min={ss.min():.4f},  max={ss.max():.4f}")
print(f"NaN in output: {np.isnan(z_pred).sum()}")

The execute method also accepts 'points' mode for irregular prediction locations:

# Predict at a list of arbitrary (x, y) pairs instead of a regular grid
x_new = np.array([x_m.mean(), x_m.mean() + 5000])
y_new = np.array([y_m.mean(), y_m.mean() - 3000])
z_pts, ss_pts = ok.execute('points', x_new, y_new)

Step 5 — Visualise Predictions and Variance #

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

im0 = axes[0].pcolormesh(gridx, gridy, z_pred, cmap='viridis', shading='auto')
axes[0].scatter(x_m, y_m, c=z, edgecolors='white', linewidth=0.5, s=40, zorder=5)
fig.colorbar(im0, ax=axes[0], label='Predicted value')
axes[0].set_title('Ordinary Kriging — Predicted Surface')
axes[0].set_xlabel('Easting (m)')
axes[0].set_ylabel('Northing (m)')

im1 = axes[1].pcolormesh(gridx, gridy, ss, cmap='plasma', shading='auto')
fig.colorbar(im1, ax=axes[1], label='Kriging variance')
axes[1].scatter(x_m, y_m, c='white', edgecolors='black', linewidth=0.5, s=40, zorder=5)
axes[1].set_title('Kriging Variance — Uncertainty Map')
axes[1].set_xlabel('Easting (m)')

plt.tight_layout()
plt.savefig('kriging_output.png', dpi=150, bbox_inches='tight')
plt.show()

Interpreting the Output #

Predicted surface (z_pred): A 2D array shaped (len(gridy), len(gridx)). Values should stay within the range of observed z for well-configured models. Extreme extrapolation values beyond the sample range signal variogram misconfiguration — usually an overfitted range or a near-zero nugget on noisy data.

Kriging variance (ss): Values are in the same units as z². Variance near zero occurs at or very close to sampled locations; it rises monotonically with distance from samples, reaching the sill in data-sparse regions. Kriging variance is independent of the observed values — it depends only on sample geometry and the variogram parameters. This makes it a reliable structural uncertainty indicator but not a substitute for cross-validated prediction error.

Nugget-to-sill ratio: A ratio below 0.1 suggests negligible measurement error and high spatial continuity. A ratio above 0.4 indicates that much of the total variance is unexplained at the smallest observable lag — either genuine micro-scale variability or instrument noise. High-nugget surfaces appear smoother (the kriging weights become more uniform) and variance maps lose their contrast near samples.

Inline SVG — Kriging Workflow Diagram #

Critical Best Practices #

Always Project to a Metric CRS First #

Geographic (lat/lon) coordinates violate the Euclidean distance assumption baked into the covariance matrix. A one-degree difference in longitude shrinks as latitude increases, so a variogram fitted in degree-space will have a physically meaningless range parameter. Project to UTM, State Plane, or another metre-based CRS using pyproj or geopandas before passing any coordinate to PyKrige. Refer to stationarity and trend analysis concepts to verify stationarity after projection — a trend that appears in unprojected space may dissolve in a properly projected one, or vice versa.

Manually Specify Variogram Parameters for Production Use #

Auto-fitting (variogram_parameters omitted) uses weighted least-squares optimisation on the empirical semivariogram. This is fast for exploration but brittle with sparse, clustered, or anisotropic data. For repeatable deployments, fit parameters once using scikit-gstat or manual inspection, then hard-code them. Document the fitting rationale alongside the parameter values.

Validate with Cross-Validation Before Deployment #

Leave-one-out cross-validation (LOOCV) is straightforward with PyKrige’s built-in method:

# PyKrige 1.7+ cross_validate returns arrays of predicted vs actual
results = ok.get_kriging_points()  # inspect fitted point estimates vs observations

# Manual LOOCV for more control:
from sklearn.model_selection import LeaveOneOut
loo = LeaveOneOut()
errors = []
for train_idx, test_idx in loo.split(x_m):
    ok_cv = OrdinaryKriging(
        x_m[train_idx], y_m[train_idx], z[train_idx],
        variogram_model='spherical',
        variogram_parameters=[4.5, 45000.0, 0.2],
        verbose=False
    )
    z_cv, _ = ok_cv.execute('points', x_m[test_idx], y_m[test_idx])
    errors.append(float(z_cv) - z[test_idx][0])

errors = np.array(errors)
print(f"LOOCV RMSE : {np.sqrt(np.mean(errors**2)):.4f}")
print(f"Mean error : {errors.mean():.4f}  (near 0 = unbiased)")

A mean standardised error near zero and variance near one indicate a calibrated model. Pair this with the cross-validation strategies available for spatial data — particularly spatial k-fold, which avoids spatial autocorrelation leakage between training and test folds.

Account for Anisotropy in Directional Phenomena #

Wind-driven pollutant dispersion, river valley morphology, and glacial deposits often exhibit stronger spatial continuity along one axis than the perpendicular direction. Ignoring anisotropy inflates variance orthogonal to the dominant structure. PyKrige exposes anisotropy_scaling and anisotropy_angle parameters:

ok_aniso = OrdinaryKriging(
    x_m, y_m, z,
    variogram_model='spherical',
    variogram_parameters=[4.5, 45000.0, 0.2],
    anisotropy_scaling=2.5,   # ratio of major to minor range axis
    anisotropy_angle=35.0,    # angle of major axis (degrees clockwise from north)
    verbose=False
)

Determine the anisotropy angle and scaling from a directional variogram analysis in scikit-gstat or gstools before fixing these values.

Scale to Large Datasets with Moving-Window Kriging #

Exact kriging inverts an $N \times N$ covariance matrix, so computational cost scales as $\mathcal{O}(N^3)$. For datasets beyond roughly 5,000 points this becomes prohibitive in single-process Python. PyKrige provides a moving-window implementation that locally limits the covariance matrix size:

# Moving-window kriging limits neighbourhood to n_closest_points
ok_large = OrdinaryKriging(
    x_m, y_m, z,
    variogram_model='spherical',
    variogram_parameters=[4.5, 45000.0, 0.2],
    verbose=False
)
z_pred, ss = ok_large.execute(
    'grid', gridx, gridy,
    n_closest_points=50,   # use the 50 nearest samples per prediction point
    backend='loop'         # required for n_closest_points to take effect
)

For datasets in the hundreds of thousands of points, consider gstools with its FFT-based covariance estimation, which scales to $\mathcal{O}(N \log N)$.

Troubleshooting #

Next Steps #

For a broader treatment of when to choose ordinary versus universal kriging, plus the mathematical derivation of the kriging system, see Ordinary & Universal Kriging. Once your predicted surface is in place, quantify spatial uncertainty rigorously by pairing the variance map with the techniques in Uncertainty & Variance Mapping.

← Back to Ordinary & Universal Kriging

Related #

• Ordinary & Universal Kriging in Python — parent guide covering theory, model selection, and the full kriging system
• Uncertainty & Variance Mapping — turn kriging variance arrays into calibrated uncertainty surfaces
• Cross-Validation Strategies for Spatial Models — spatial k-fold and blocked CV to avoid autocorrelation leakage
• Testing for Second-Order Stationarity in Python — verify the constant-mean assumption before fitting ordinary kriging