Interpolating Sparse Yield Monitor Data with Kriging

TL;DR: Project yield points to a metric UTM CRS, fit an Ordinary Kriging model with a spherical variogram using pykrige, execute on a 10 m grid, and write a dual-band GeoTIFF — Band 1 carries predicted yield, Band 2 carries prediction variance for uncertainty masking.

Why Kriging Arises in Yield-Data Workflows

Combine harvesters log yield at 1–2 Hz, but GPS drift, header-width jumps, and equipment lag create irregular spacing and systematic gaps — particularly on headlands, around obstacles, and where the harvester decelerates. Straight-line interpolation methods such as IDW treat every sample point equally regardless of spatial structure, which causes concentric bulls-eye artefacts around dense clusters and degrades prescription accuracy by 10–25% compared to geostatistical methods on typical Midwest corn-belt datasets.

Kriging solves this by first modelling the spatial autocorrelation structure of the yield variable — how similar values are as a function of separation distance — through a variogram, and then using that model as a distance-weighting kernel. The result is a best-linear-unbiased prediction at each grid node, together with a variance surface that quantifies where predictions are unreliable. Without this uncertainty band, agronomists have no signal about which zones of the interpolated raster to trust when feeding the output into downstream management zone classification.

The diagram below shows how raw sparse monitor points are transformed through variogram fitting and grid prediction into the final raster pair.

Kriging interpolation pipeline for sparse yield monitor data Four-stage flow: sparse yield points, CRS projection and thinning, variogram fitting and Ordinary Kriging, dual-band GeoTIFF output (yield + variance). Sparse yield points CSV / GeoJSON (WGS84) Project + thin UTM / grid decimation Variogram fit + Ordinary Kriging Dual-band GeoTIFF Band 1: yield · Band 2: σ² 1. Ingest 2. Prepare 3. Interpolate 4. Export

Prerequisites

These differ from the Spatial Interpolation for Yield Data cluster page only in one package:

TEXT 
pykrige==1.7.2
geopandas==0.14.4
rasterio==1.3.10
numpy==1.26.4
scipy==1.13.0

Install with:

BASH 
pip install pykrige==1.7.2 geopandas==0.14.4 rasterio==1.3.10 numpy==1.26.4 scipy==1.13.0

Input requirements:

  • Yield points as a CSV or GeoJSON file with lon, lat, and yield_bu_ac columns (or equivalent)
  • Points must number at least 30 after cleaning; below that, leave-one-out cross-validation is unreliable
  • No pre-existing CRS projection required — the script handles the WGS84 → UTM step explicitly

Windows users: install pykrige via conda install -c conda-forge pykrige to avoid GDAL build conflicts with the pip wheel.

Step-by-Step

Step 1 — Load and reproject to UTM

Kriging’s covariance function measures Euclidean distance in the units of the coordinate system. Using geographic coordinates means that 1° of longitude narrows towards the poles, producing anisotropic variograms and wrong range estimates. Always reproject to the field’s UTM zone before fitting.

Reproject to your target UTM zone using the approach described in understanding CRS in precision agriculture. Replace EPSG:32615 (UTM zone 15 N) with the code for your farm’s location.

Step 2 — Define the output grid

A 10 m cell size matches the practical spatial resolution of most yield monitor GPS logs (post-filtering) and aligns with the grid spacing that variable-rate controllers can act on. Use numpy.arange rather than numpy.linspace so that the grid origin is pinned to xmin and the spacing is exact.

Step 3 — Fit variogram and run Ordinary Kriging

OrdinaryKriging auto-fits a spherical variogram using weighted least-squares when weight=True. Set nlags=12 for moderate-density fields; increase to 20 for dense urban test plots. The execute("grid", grid_x, grid_y) call returns both the prediction array z_interp and the variance surface ss.

Step 4 — Export dual-band GeoTIFF

Writing Band 2 (variance) alongside the prediction is essential. Variance values above the sample variance indicate regions where the model is extrapolating rather than interpolating — those zones should be masked before the raster is passed to prescription generation software. Use nodata=np.nan so that downstream tools such as gdalwarp and rasterio.mask handle edge clipping correctly.

Step 5 — Cross-validate with leave-one-out RMSE

The complete, directly runnable script:

PYTHON 
import numpy as np
import geopandas as gpd
import rasterio
from rasterio.transform import from_origin
from pykrige.ok import OrdinaryKriging

# ── 1. Load and clean ──────────────────────────────────────────────────────
gdf = gpd.read_file("sparse_yield_points.geojson")
# Drop rows with missing yield values before any spatial operation
gdf = gdf.dropna(subset=["yield_bu_ac"])
assert len(gdf) >= 30, f"Too few points after cleaning: {len(gdf)}"

# Ensure WGS84 before reprojection
if gdf.crs is None:
    gdf = gdf.set_crs("EPSG:4326")
elif gdf.crs.to_epsg() != 4326:
    gdf = gdf.to_crs("EPSG:4326")

# ── 2. Reproject to metric UTM (replace 32615 with your zone) ─────────────
UTM_EPSG = 32615
gdf_utm = gdf.to_crs(f"EPSG:{UTM_EPSG}")
assert gdf_utm.crs.is_projected, "CRS must be projected (metric) for Kriging"

x = gdf_utm.geometry.x.values
y = gdf_utm.geometry.y.values
z = gdf_utm["yield_bu_ac"].values.astype(np.float64)

# ── 3. Build output grid at 10 m resolution ───────────────────────────────
GRID_RES = 10.0
xmin, ymin, xmax, ymax = gdf_utm.total_bounds
grid_x = np.arange(xmin, xmax, GRID_RES)
grid_y = np.arange(ymin, ymax, GRID_RES)
cols, rows = len(grid_x), len(grid_y)

# ── 4. Fit spherical variogram and execute Ordinary Kriging ───────────────
ok = OrdinaryKriging(
    x, y, z,
    variogram_model="spherical",
    nlags=12,
    weight=True,       # weight lags by number of pairs (critical on sparse data)
    verbose=False,
    enable_plotting=False,
)
z_pred, z_var = ok.execute("grid", grid_x, grid_y)

# Flip so row 0 is northernmost (rasterio convention: origin = top-left)
z_pred = np.flipud(z_pred.data)
z_var  = np.flipud(z_var.data)

# ── 5. Leave-one-out cross-validation ─────────────────────────────────────
residuals = []
for i in range(len(x)):
    mask = np.ones(len(x), dtype=bool)
    mask[i] = False
    ok_loo = OrdinaryKriging(
        x[mask], y[mask], z[mask],
        variogram_model="spherical",
        nlags=12,
        weight=True,
        verbose=False,
        enable_plotting=False,
    )
    pred, _ = ok_loo.execute("points", np.array([x[i]]), np.array([y[i]]))
    residuals.append(z[i] - pred[0])

rmse = float(np.sqrt(np.mean(np.array(residuals) ** 2)))
print(f"LOO-CV RMSE: {rmse:.2f} bu/ac")
# Fail loudly if RMSE exceeds 20 % of mean yield — indicates a bad variogram fit
assert rmse / z.mean() < 0.20, f"RMSE too high ({rmse:.2f}); inspect variogram parameters"

# ── 6. Export dual-band GeoTIFF ───────────────────────────────────────────
transform = from_origin(xmin, ymax, GRID_RES, GRID_RES)
with rasterio.open(
    "yield_kriged.tif",
    "w",
    driver="GTiff",
    height=rows,
    width=cols,
    count=2,
    dtype="float32",
    crs=gdf_utm.crs,
    transform=transform,
    nodata=np.nan,
) as dst:
    dst.write(z_pred.astype("float32"), 1)  # Band 1: predicted yield (bu/ac)
    dst.write(z_var.astype("float32"),  2)  # Band 2: prediction variance
    dst.update_tags(
        1, description="Ordinary Kriging prediction — yield bu/ac"
    )
    dst.update_tags(
        2, description="Kriging variance — mask where > sample_variance"
    )

print(f"Output: yield_kriged.tif  ({cols}×{rows} px, {GRID_RES} m resolution)")

Inline verification: after running the script, open the output file and confirm the spatial extent matches expectations:

PYTHON 
with rasterio.open("yield_kriged.tif") as src:
    print(src.profile)
    assert src.count == 2, "Expected two bands"
    assert src.crs.to_epsg() == UTM_EPSG
    band1 = src.read(1, masked=True)
    print(f"Yield range: {band1.min():.1f} – {band1.max():.1f} bu/ac")

Gotchas and Edge Cases

  • Running Kriging on geographic coordinates (WGS84) is the single most common error. The variogram range will be expressed in degrees, the fitted sill will mismatch the sample variance, and predictions will be silently wrong. The assert gdf_utm.crs.is_projected guard in Step 2 catches this at runtime.

  • pykrige scales at O(n³) due to the dense covariance matrix inversion (scipy.linalg.solve). For fields with more than 8 000 cleaned points, spatially subsample with gdf.sample(frac=0.5, random_state=42) or switch to gstools which supports sparse solvers. The LOO-CV loop in Step 5 compounds this — consider running it only on a random 10 % subset for production jobs.

  • np.flipud is required before writing. pykrige returns the grid with grid_y[0] as the southernmost row, but rasterio.open expects the first row to correspond to ymax (north). Skipping this flip produces a vertically mirrored raster that passes all CRS checks but is spatially inverted — visible only when overlaid on a basemap.

  • Masking high-variance zones before prescription generation is not optional. Where yield monitor passes are more than 150 m apart, the kriging variance can exceed the sample variance of the whole field, meaning the model is pure extrapolation. Clip these pixels with np.where(z_var > z.var(), np.nan, z_pred) before passing the raster to downstream management zone workflows.

Parent Guide

This guide is part of Spatial Interpolation for Yield Data — see there for the full pipeline context including data cleaning, IDW fallback, and integration with field boundary masking.

Parent section Spatial Interpolation for Yield Data in Python