SciKit-GStat - Geostatistics

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

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

When to Use vs Alternatives

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

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