Geostatistical Modeling Spatial Statistics with Python

Fitting Empirical Variograms with SciKit-GStat

TL;DR: Pass your coordinate array and value array straight into Variogram(coordinates, values, n_lags=15, model="spherical"). The object lazily computes everything: V.experimental holds the semivariances, V.bins the lag distances, V.plot() draws the cloud with the fitted curve, and V.parameters returns [range, sill, nugget]. Call V.describe() for the full numeric summary.

Why This Matters

The empirical variogram is the measured fingerprint of spatial continuity in your data, and every kriging surface you later produce inherits its range, sill, and nugget. Getting the estimation right is the whole point of empirical variogram estimation, the first stage of the broader variogram modeling and semivariance analysis workflow. scikit-gstat collapses the pairwise distance calculation, binning, estimator choice, and model fitting into a single Variogram object, so you can move from raw observations to a defensible model in a handful of lines.

The experimental semivariance for a lag class hh is the classical Matheron estimator:

γ^(h)=12N(h)∑(i,j)∈N(h)(z(xi)−z(xj))2\hat{\gamma}(h) = \frac{1}{2N(h)} \sum_{(i,j) \in N(h)} \left(z(x_i) - z(x_j)\right)^2

where N(h)N(h) is the set of point pairs separated by a distance falling inside lag class hh, and z(xi)z(x_i) is the observed value at location xix_i. scikit-gstat evaluates this for every lag class and stores the results in V.experimental.

Environment

bash
pip install "scikit-gstat>=1.0" "numpy>=1.23" "matplotlib>=3.7"
python
import numpy as np
import matplotlib.pyplot as plt
from skgstat import Variogram

scikit-gstat pulls in scipy and numpy automatically. matplotlib is only needed for V.plot(); the numeric attributes work without it.

Step-by-Step Implementation

Step 1 — Assemble coordinates and values

python
rng = np.random.default_rng(7)

# 220 sample locations in a 1000 m x 1000 m block (projected, metric CRS)
coords = rng.uniform(0, 1000, size=(220, 2))

# A regionalized variable with a smooth trend plus short-range noise
field = (
    8.0
    + 0.004 * coords[:, 0]                       # gentle east-west gradient
    + 3.0 * np.sin(coords[:, 1] / 180.0)         # medium-scale structure
    + rng.normal(0, 0.8, size=len(coords))       # nugget-scale noise
)

coordinates must be an (N, 2) array of metric coordinates and values a matching 1-D array of length N. Use a projected CRS: distances in degrees are not isotropic and will distort every lag.

Step 2 — Construct the Variogram

python
V = Variogram(
    coordinates=coords,
    values=field,
    n_lags=15,            # number of lag classes
    model="spherical",    # theoretical model fitted to the cloud
    estimator="matheron", # classical semivariance estimator
    maxlag="median",      # cap the largest lag at the median pairwise distance
    normalize=False,
)
print(V)

Constructing the object does not eagerly compute anything expensive; the pairwise distances, binning, and fit are all evaluated lazily the first time you access a result attribute. Printing V triggers the fit and shows a compact summary.

Step 3 — Inspect the experimental variogram

python
lag_distances = V.bins          # upper edge of each lag class (metres)
semivariances = V.experimental  # semivariance for each lag class

for d, g in zip(lag_distances, semivariances):
    print(f"lag <= {d:7.1f} m   gamma = {g:6.3f}")

V.bins and V.experimental are aligned arrays of length n_lags. They are the raw material of the plot and of the model fit; everything downstream is derived from them.

Step 4 — Check how many pairs fall in each lag

python
pairs_per_lag = np.bincount(V.lag_groups(), minlength=V.n_lags)
print("pairs per lag:", pairs_per_lag)

V.lag_groups() returns the lag-class index of every point pair, so np.bincount gives the pair count per bin. Bins with fewer than roughly 30 pairs produce unstable semivariances; if the tail is sparse, lower maxlag or reduce n_lags.

Step 5 — Plot the cloud and fitted model

python
fig = V.plot(hist=True, show=False)
fig.savefig("empirical_variogram.png", dpi=150, bbox_inches="tight")

V.plot() overlays the fitted theoretical model on the experimental points, and hist=True adds the pairs-per-lag histogram beneath. A clean fit rises from near the nugget, bends through the experimental points, and levels off at the sill.

Step 6 — Extract range, sill and nugget

python
effective_range, sill, nugget = V.parameters
print(f"effective range : {effective_range:8.2f} m")
print(f"sill            : {sill:8.3f}")
print(f"nugget          : {nugget:8.3f}")
print(V.describe())   # full dict: model, estimator, fit metrics, parameters

Interpreting the Output

V.parameters returns [effective_range, sill, nugget]. The effective range is the separation distance at which spatial correlation effectively vanishes, so points farther apart than the range are treated as independent. The sill is the semivariance plateau, approximately equal to the sample variance of a stationary variable. The nugget is the intercept at zero lag; a nugget greater than zero encodes measurement error plus variation at scales finer than your sampling spacing.

Attribute Meaning Practical read
V.experimental Semivariance per lag class Should rise then flatten toward the sill
V.bins Lag distances (upper edges) x-axis of the variogram
V.parameters [range, sill, nugget] The three numbers you carry into kriging
V.rmse Fit residual against experimental points Lower is a tighter fit
V.describe() Full summary dictionary Model name, estimator, and all fit metrics

If the experimental points never flatten, the variable is likely non-stationary and carries a trend that must be removed first, a diagnostic covered under stationarity and trend analysis.

Critical Best Practices

Cap maxlag at roughly half the domain

Semivariances at large separation distances rest on very few pairs and swing wildly. Setting maxlag="median" or an explicit value near half the domain diameter keeps the fit anchored to the well-sampled short and medium lags. Never let the variogram extend to the full diagonal of the study area.

Choose the estimator to match your data

estimator="matheron" is the classical choice but is sensitive to outliers because it squares differences. For skewed or contaminated data, estimator="cressie" (the Cressie-Hawkins robust estimator) or estimator="dowd" down-weights extreme pairs and yields a more stable sill.

Verify pair counts before trusting the fit

A beautiful-looking curve fitted through bins of five pairs each is fiction. Always compute np.bincount(V.lag_groups()) and confirm each lag you fit to holds enough pairs. Sparse tails are the single most common cause of an inflated range.

Keep coordinates in a projected CRS

scikit-gstat computes Euclidean distances directly on the coordinate array. Latitude and longitude degrees are not metric and are anisotropic away from the equator, so reproject to UTM or a national grid first. See GeoPandas data preparation for the reprojection workflow.

Re-fit rather than rebuild when tuning

Changing V.model, V.n_lags, or V.estimator on an existing object invalidates the cache and re-fits automatically. You do not need to reconstruct the Variogram from scratch when experimenting, which keeps parameter sweeps fast.

Troubleshooting

Symptom Likely cause Fix
ValueError on construction coordinates not shaped (N, 2) or length mismatch with values Reshape coordinates and confirm len(values) == len(coords)
Experimental variogram never flattens Trend or non-stationarity in the data Detrend before estimation; revisit stationarity diagnostics
Very noisy far lags Too few pairs per bin Lower maxlag or reduce n_lags
V.parameters range implausibly large Fit driven by sparse tail bins Cap maxlag, check np.bincount(V.lag_groups())
Sill far above sample variance Outliers inflating Matheron estimator Switch to estimator="cressie" or "dowd"
V.plot() raises on a headless server No display backend Call matplotlib.use("Agg") before importing pyplot

Next Steps

With a clean experimental variogram in hand, tune the lag geometry that controls it in choosing lag bins and bandwidth for variograms, then compare candidate theoretical curves in fitting spherical, exponential and Gaussian variogram models.

Related

← Back to Empirical Variogram Estimation

Frequently Asked Questions

What do V.bins and V.experimental return in scikit-gstat?

V.bins returns the upper edge of each lag class as an array of separation distances, and V.experimental returns the matching semivariance value for each of those lag classes. Together they are the x and y coordinates of the experimental variogram points that the theoretical model is fitted to.

Does scikit-gstat report the effective range or the model range parameter?

V.parameters[0] is the effective range, the distance at which the model reaches about 95 percent of the sill. scikit-gstat rescales the internal shape parameter of exponential and Gaussian models so that the reported range is directly comparable across model types.

Why does my empirical variogram look noisy at large lags?

Large separation distances contain few point pairs, so their semivariance estimates have high variance. Restrict maxlag to about half the domain diameter and require roughly 30 or more pairs per bin so the far lags do not distort the fit.

Related

Choosing Lag Bins and Bandwidth for Variograms