Back
Kriging Interpolation
Overview
Kriging is a geostatistical interpolation method that provides the Best Linear Unbiased Estimate (BLUE) of a spatially distributed variable. Unlike simpler methods like IDW, Kriging takes into account the spatial correlation structure of the data through a variogram model, and provides both an estimate and an associated variance (uncertainty) at each prediction point.
Variogram Model
Variogram Types and Parameters
// Available variogram models
enum class VariogramModel {
Spherical, // Reaches sill at the range distance
Exponential, // Approaches sill asymptotically
Gaussian // Smooth, parabolic behavior near origin
};
// Variogram parameters
struct VariogramParams {
double nugget = 0.0; // Nugget effect (discontinuity at origin)
double sill = 1.0; // Sill (total variance)
double range = 1.0; // Range (correlation distance)
VariogramModel model = VariogramModel::Spherical;
};
The variogram describes how spatial correlation decreases with distance:
- Nugget: Represents measurement error or micro-scale variation. A non-zero nugget means the variogram does not pass through the origin.
- Sill: The variance plateau. The variogram levels off at this value.
- Range: The distance at which the variogram reaches (or approaches) the sill. Beyond this distance, data points are considered uncorrelated.
Variogram Functions
Variogram Functions
// Evaluate a variogram model at a given distance
double variogram_model(double distance, const VariogramParams ¶ms);
// Calculate the experimental variogram from data
// Returns pairs of (distance, semivariance)
std::vector<std::pair<double, double>> calculate_experimental_variogram(
const Serie<Vector<2>> &points,
const Serie<double> &values,
size_t num_bins = 20
);
Experimental Variogram Example
// Compute the experimental variogram from observed data
df::Serie<Vector<2>> obs_pts = /* observation locations */;
df::Serie<double> obs_vals = /* measured values */;
auto exp_variogram = df::geo::calculate_experimental_variogram(
obs_pts, obs_vals, 15
);
// Plot or inspect the experimental variogram
for (const auto &[dist, semivar] : exp_variogram) {
std::cout << "Distance: " << dist
<< " Semivariance: " << semivar << std::endl;
}
// Fit variogram parameters (manual in this example)
df::VariogramParams params;
params.nugget = 0.1;
params.sill = 1.5;
params.range = 50.0;
params.model = df::VariogramModel::Spherical;
Ordinary Kriging
Ordinary Kriging Function
// Perform ordinary kriging interpolation
// Returns a pair: {estimated_values, kriging_variances}
std::pair<Serie<double>, Serie<double>> ordinary_kriging(
const Serie<Vector<2>> &points, // Known data locations
const Serie<double> &values, // Known data values
const Serie<Vector<2>> &targets, // Prediction locations
const VariogramParams ¶ms // Fitted variogram parameters
);
Ordinary Kriging Example
// Known measurement points
df::Serie<Vector<2>> stations{
{10, 20}, {30, 10}, {50, 40}, {20, 50}, {40, 30}
};
df::Serie<double> measurements{5.2, 3.8, 7.1, 6.5, 4.9};
// Variogram parameters (from experimental variogram fitting)
df::VariogramParams vparams;
vparams.nugget = 0.2;
vparams.sill = 4.0;
vparams.range = 30.0;
vparams.model = df::VariogramModel::Exponential;
// Target grid points
auto grid = df::grid::cartesian::from_dims<2>(
{50, 50}, {30.0, 30.0}, {60.0, 60.0}
);
// Perform kriging
auto [estimates, variances] = df::geo::ordinary_kriging(
stations, measurements, grid, vparams
);
// estimates: predicted values at grid points
// variances: kriging variance (uncertainty) at each point
// Low variance = high confidence, high variance = low confidence
for (size_t i = 0; i < 5; ++i) {
std::cout << "Point " << i
<< ": estimate=" << estimates[i]
<< " +/- " << std::sqrt(variances[i]) << std::endl;
}
Complete Example
Kriging Workflow
#include <dataframe/Serie.h>
#include <dataframe/geo/kriging.h>
#include <dataframe/geo/grid/from_dims.h>
#include <dataframe/math.h>
#include <iostream>
#include <cmath>
int main() {
// Weather station data
df::Serie<Vector<2>> stations{
{0, 0}, {100, 0}, {0, 100}, {100, 100}, {50, 50},
{25, 75}, {75, 25}, {30, 30}, {70, 70}
};
df::Serie<double> rainfall{12.0, 8.5, 15.2, 10.1, 11.8,
14.0, 9.2, 12.5, 11.0};
// Step 1: Compute experimental variogram
auto exp_vario = df::geo::calculate_experimental_variogram(
stations, rainfall, 10
);
// Step 2: Fit a variogram model
df::VariogramParams params;
params.nugget = 0.5;
params.sill = 6.0;
params.range = 60.0;
params.model = df::VariogramModel::Gaussian;
// Step 3: Create prediction grid
auto grid = df::grid::cartesian::from_dims<2>(
{100, 100}, {50.0, 50.0}, {100.0, 100.0}
);
// Step 4: Kriging interpolation
auto [pred, var] = df::geo::ordinary_kriging(
stations, rainfall, grid, params
);
// Step 5: Analyze results
auto [min_pred, max_pred] = df::math::bounds(pred);
std::cout << "Rainfall predictions: [" << min_pred
<< ", " << max_pred << "] mm\n";
auto [min_var, max_var] = df::math::bounds(var);
std::cout << "Kriging std dev: [" << std::sqrt(min_var)
<< ", " << std::sqrt(max_var) << "] mm\n";
return 0;
}