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 &params);

// 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    &params     // 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;
}