scikit-gstat

Quick Reference

import skgstat as skg
import numpy as np

# Create variogram
V = skg.Variogram(coordinates=coords, values=values, n_lags=15)

# Fit model
V.model = 'spherical'
print(f"Range: {V.parameters[0]:.2f}, Sill: {V.parameters[1]:.2f}")

# Kriging interpolation
ok = skg.OrdinaryKriging(V)
predictions = ok.transform(grid_coords)

Key Classes

| Class | Purpose | |-------|---------| | Variogram | Empirical and theoretical variograms | | OrdinaryKriging | Interpolation with spatial correlation | | DirectionalVariogram | Anisotropic variograms | | SpaceTimeVariogram | Spatio-temporal analysis |

Essential Operations

Create and Fit Variogram

import skgstat as skg

V = skg.Variogram(
    coordinates=coords,      # (n, 2) array of x, y
    values=values,           # (n,) array of measurements
    n_lags=15,
    maxlag='median'          # or specific distance
)

# Fit model: 'spherical', 'exponential', 'gaussian', 'matern', 'stable'
V.model = 'spherical'

# Get parameters
print(f"Range: {V.parameters[0]:.2f}")
print(f"Sill: {V.parameters[1]:.2f}")
print(f"Nugget: {V.parameters[2]:.2f}")
print(f"RMSE: {V.rmse:.4f}")

Ordinary Kriging

import skgstat as skg
import numpy as np

V = skg.Variogram(coords, values, model='spherical')
ok = skg.OrdinaryKriging(V)

# Create prediction grid
x = np.linspace(0, 100, 50)
y = np.linspace(0, 100, 50)
xx, yy = np.meshgrid(x, y)
grid_coords = np.column_stack([xx.ravel(), yy.ravel()])

# Predict
predictions = ok.transform(grid_coords)
Z = predictions.reshape(xx.shape)

# Get variance
ok.return_variance = True
predictions, variance = ok.transform(grid_coords)

Directional Variogram

import skgstat as skg

DV = skg.DirectionalVariogram(
    coordinates=coords,
    values=values,
    azimuth=45,          # Direction in degrees
    tolerance=22.5,      # Angular tolerance
    bandwidth='q33'      # Perpendicular bandwidth
)

# Check anisotropy
for az in [0, 45, 90, 135]:
    DV.azimuth = az
    print(f"Azimuth {az}: Range = {DV.parameters[0]:.2f}")

Cross-Validation

import skgstat as skg
from sklearn.model_selection import cross_val_score

V = skg.Variogram(coords, values, model='spherical')
ok = skg.OrdinaryKriging(V)

scores = cross_val_score(ok, coords, values, cv=5, scoring='neg_mean_squared_error')
print(f"CV RMSE: {np.sqrt(-scores.mean()):.4f}")

Robust Estimators

import skgstat as skg

# Use robust estimator for noisy data
V = skg.Variogram(
    coords, values,
    estimator='cressie'  # 'matheron', 'cressie', 'dowd', 'genton'
)

Quick Model Reference

| Model | Behavior | |-------|----------| | spherical | Most common, linear near origin | | exponential | Never reaches sill, gradual approach | | gaussian | Parabolic near origin, smooth | | matern | Flexible smoothness control |

When to Use vs Alternatives

| Use Case | Tool | Why | |----------|------|-----| | Variogram analysis + kriging | scikit-gstat | Modern API, sklearn-compatible | | GSLIB-style simulation (SGSIM) | GeostatsPy | Full GSLIB simulation engine | | Kriging with trend/drift | pykrige | Universal kriging, regression kriging | | Random field generation | gstools | Flexible covariance, SRF generation | | Spatio-temporal variograms | scikit-gstat | Built-in SpaceTimeVariogram | | Production geomodelling | SGeMS / Petrel | GUI, large-scale 3D models | | Robust variogram estimation | scikit-gstat | Cressie, Dowd, Genton estimators | | ML pipeline integration | scikit-gstat | sklearn fit/transform interface |

Choose scikit-gstat when: You want a Pythonic, scikit-learn-compatible API for variogram fitting and kriging. Best for exploratory geostatistical analysis with cross-validation and integration into ML pipelines.

Choose GeostatsPy when: You need GSLIB-compatible simulation workflows (SGSIM, SISIM) or are working with traditional geostatistical conventions.

Choose pykrige when: You need universal kriging with external drift variables or regression kriging combining geostatistics with machine learning predictions.

Common Workflows

Variogram Fitting and Ordinary Kriging

• [ ] Load spatial data as numpy arrays (coordinates and values)
• [ ] Create Variogram object with appropriate n_lags and maxlag
• [ ] Test estimators: matheron (default) vs cressie (robust) for noisy data
• [ ] Fit multiple models (spherical, exponential, gaussian) and compare RMSE
• [ ] Check anisotropy with DirectionalVariogram at 0, 45, 90, 135 degrees
• [ ] Select best model based on RMSE and visual fit
• [ ] Create OrdinaryKriging object from fitted variogram
• [ ] Define prediction grid and run ok.transform(grid_coords)
• [ ] Set ok.return_variance = True to get kriging variance
• [ ] Cross-validate with cross_val_score() to assess prediction quality
• [ ] Map predictions and kriging variance

Common Issues

| Issue | Solution | |-------|----------| | Variogram flat or erratic | Adjust n_lags and maxlag (try maxlag='median') | | Poor model fit (high RMSE) | Try different model types or nested structures | | Kriging too slow | Reduce number of conditioning points or grid resolution | | Nugget too large | May indicate measurement error; try robust estimators | | Anisotropy unclear | Use smaller angular tolerance in DirectionalVariogram |

Tips

1. Maxlag should be ~50% of study area diagonal
2. Use robust estimators (cressie, dowd) with noisy data
3. Test multiple models and compare RMSE
4. Check anisotropy with directional variograms before kriging

References

• Variogram Models - Model equations and parameters
• Kriging Methods - Kriging types and configuration

Scripts

• scripts/variogram_analysis.py - Complete variogram analysis workflow