Empirical Variogram Estimation in Python

The empirical (or experimental) variogram is the first quantitative product of any geostatistical study: it turns a cloud of scattered samples into a curve that shows how spatial similarity decays with distance. Everything downstream — the fitted model, the kriging weights, the prediction intervals — is only as trustworthy as this estimate. This guide covers how to compute it correctly in Python, within the broader context of variogram modeling and semivariance analysis, using both the classical Matheron estimator and the robust Cressie-Hawkins estimator.

Prerequisites #

• Python 3.10+
• scikit-gstat>=0.6, numpy>=1.22, geopandas>=0.14, matplotlib>=3.6
• A GeoDataFrame of point observations in a projected (metric) CRS — UTM or a regional projection, never EPSG:4326
• Values that are already detrended and declustered, so the intrinsic stationarity assumption holds
• At least ~100 well-distributed samples so that each lag bin can hold 30–50 pairs

Mathematical Core #

Two estimators dominate practice. The classical Matheron estimator averages half the squared differences of all pairs whose separation falls in the lag class $N(h)$:

$$\hat{\gamma}(h) = \frac{1}{2\,|N(h)|} \sum_{(i,j) \in N(h)} \left[ z(s_i) - z(s_j) \right]^{2}$$

Here $z(s_i)$ is the observed value at location $s_i$, $N(h)$ is the set of pairs separated by a distance within the tolerance of lag $h$, and $|N(h)|$ is the number of such pairs. Because the differences are squared, a single outlying pair can dominate a bin.

The Cressie-Hawkins robust estimator replaces the squared differences with fourth roots of absolute differences and applies a bias correction:

$$\bar{\gamma}(h) = \frac{1}{2}\cdot\frac{\left(\dfrac{1}{|N(h)|}\displaystyle\sum_{(i,j)\in N(h)} |z(s_i) - z(s_j)|^{1/2}\right)^{4}}{0.457 + \dfrac{0.494}{|N(h)|} + \dfrac{0.045}{|N(h)|^{2}}}$$

The denominator constant $0.457 + 0.494/|N(h)| + 0.045/|N(h)|^2$ makes the estimator approximately unbiased for a Gaussian field. Use Cressie-Hawkins whenever the value distribution has a heavy upper tail; use Matheron for clean, near-Gaussian data where its lower variance is an advantage.

The fourth-root construction is doing something specific. Raising each difference to the power $\tfrac{1}{2}$ compresses the upper tail before averaging, so a pair whose difference is ten times the typical one contributes about 3.2 times as much to the mean rather than the 100 times it would contribute after squaring. Raising the compressed mean back to the fourth power restores variance units, and the constant $0.457$ undoes the bias that this compress-then-expand round trip introduces for a Gaussian field. The two remaining terms matter only for small bins: at $|N(h)| = 30$ they add about 3.5% to the denominator, while at $|N(h)| = 6$ they add roughly 15%. This is also where the estimator degrades. The correction is derived under Gaussianity, so on a strongly skewed field Cressie-Hawkins stays robust to individual outlying pairs but is no longer unbiased for the true semivariance. If the field is lognormal, transform first and estimate the variogram of the logs rather than asking the robust estimator to absorb the skew on its own.

Annotated Implementation #

scikit-gstat centres on one class, Variogram, whose signature exposes every estimation choice. The essential call is Variogram(coordinates, values, ...), where coordinates is an (n, 2) array of projected coordinates and values is a length-n array of observations.

1. Extract coordinates and values #

import geopandas as gpd
import numpy as np

# Point samples in a projected, metric CRS (UTM Zone 33N here)
gdf = gpd.read_file("soil_zinc.gpkg").to_crs("EPSG:32633")

# scikit-gstat wants a plain (n, 2) coordinate array and a length-n value array
coordinates = np.column_stack([gdf.geometry.x.values, gdf.geometry.y.values])
values = gdf["zinc"].to_numpy(dtype=float)

assert coordinates.shape[0] == values.shape[0]
assert not np.isnan(values).any(), "Remove or impute NaNs before estimation"

The to_crs call is not optional. Variogram computes Euclidean distances straight from the coordinate array, so degrees would produce meaningless lag bins.

2. Instantiate the Variogram #

from skgstat import Variogram

V = Variogram(
    coordinates,
    values,
    estimator="matheron",   # classical estimator; use "cressie" for robust
    model="spherical",      # theoretical model fitted to the binned points
    bin_func="even",        # equal-width lag classes
    n_lags=15,              # number of lag bins
    maxlag=0.5,             # cap at half the maximum pairwise distance
    normalize=False,        # keep lags and semivariance in data units
)

print(V)                    # prints estimator, model, fitted nugget/sill/range

The constructor immediately estimates the empirical variogram and fits the named model. To work with the empirical estimate alone, read V.bins (lag-class upper edges) and V.experimental (the estimated semivariance per bin) before trusting the fit.

3. Switch to the robust estimator #

# Rebuild with the Cressie-Hawkins robust estimator for heavy-tailed data
V_robust = Variogram(
    coordinates, values,
    estimator="cressie",
    model="spherical",
    n_lags=15,
    maxlag=0.5,
)

# Compare the two empirical curves directly
for h, g_m, g_r in zip(V.bins, V.experimental, V_robust.experimental):
    print(f"lag {h:8.1f} m   matheron {g_m:8.2f}   cressie {g_r:8.2f}")

If the two curves diverge sharply at particular lags, those lags contain influential outlier pairs and the robust estimate should be preferred.

4. Inspect the variogram cloud and pair counts #

# Per-bin pair counts — the single most important diagnostic
counts = np.fromiter((len(g) for g in V.lag_classes()), dtype=int)
for h, c in zip(V.bins, counts):
    flag = "  <-- under-supported" if c < 30 else ""
    print(f"lag {h:8.1f} m : {c:5d} pairs{flag}")

# The raw cloud: one half-squared-difference per pair, and each pair's distance
distances = V.distance          # condensed pairwise distances
cloud = 0.5 * np.square(V.pairwise_diffs)   # per-pair semivariances
print(f"cloud points: {cloud.size}, expected n*(n-1)/2 = "
      f"{values.size * (values.size - 1) // 2}")

The cloud is the un-aggregated scatter behind the binned curve; scanning it reveals the outlier pairs that motivate the robust estimator.

Diagnostic Configuration #

Three parameters control the quality of the empirical estimate. Tune them deliberately rather than accepting defaults.

Number of lag bins (n_lags). Too few bins over-smooth and hide the range; too many leave each bin under-populated and noisy. Ten to twenty bins is typical. Sweep the value and watch the pair counts.

Maximum lag (maxlag). Restrict to roughly half the domain diameter. scikit-gstat accepts a float (absolute distance), a fraction in $(0, 1)$ of the maximum pairwise distance, or the string "median" for the median pairwise distance:

# Absolute cap at 1500 m; often clearer than a fraction of the max distance
V_capped = Variogram(coordinates, values, maxlag=1500.0, n_lags=12)

# Equal-count binning: every bin holds the same number of pairs, ideal for
# unevenly sampled data where equal-width bins leave sparse tails.
V_uniform = Variogram(coordinates, values, bin_func="uniform", n_lags=12)

Minimum pairs per bin. scikit-gstat does not hard-enforce a floor, so verify it yourself. Bins with fewer than 30 pairs should be dropped or absorbed by widening bins or lowering maxlag:

def flag_thin_bins(V, min_pairs=30):
    """Return the lag distances whose bins are under-supported."""
    counts = np.fromiter((len(g) for g in V.lag_classes()), dtype=int)
    thin = V.bins[counts < min_pairs]
    if thin.size:
        print(f"WARNING: {thin.size} bins below {min_pairs} pairs: {thin}")
    return thin

flag_thin_bins(V)

Empirical Cloud versus Binned Variogram #

The cloud on the left is what V.distance and V.pairwise_diffs expose; the binned curve on the right is V.bins against V.experimental. Model fitting operates only on the binned points, which is why their quality — driven by pair counts — matters more than the raw cloud.

Output Interpretation #

Read the binned empirical variogram before fitting anything:

• Rise from a positive intercept: the intercept at the shortest lag approximates the nugget — measurement error plus sub-sample-spacing variability. A large intercept relative to the plateau means a noisy or under-sampled process.
• Plateau (sill): where the curve flattens marks the total variance and, by $\gamma(h) = C(0) - C(h)$, the point beyond which pairs are uncorrelated. A curve that never flattens signals residual trend — stop and revisit stationarity and trend analysis before proceeding.
• Range: the lag at which the plateau is reached is the correlation range, the key input to kriging and to distance-band weights.
• Noisy tail: erratic points at the widest lags are the tell-tale of thin bins. They should be excluded by capping maxlag, not fitted through.

Good signs are a smooth monotone rise to a clear plateau with every bin holding ample pairs. Warning signs are a sawtooth curve (too many bins, too few pairs), an unbounded climb (trend), or a flat line at the sill from lag zero (pure nugget — no resolvable spatial structure).

The ratio of nugget to sill is worth reading as a number rather than only as a shape. A nugget below about 25% of the sill indicates strong spatial dependence, and kriging will meaningfully outperform a global mean; between 25% and 75% the dependence is moderate; above 75% the process is close to spatially random at the sampled resolution, and kriged predictions collapse toward the mean with prediction variances near the sill almost everywhere. The classification is a convention rather than a test, but it is a fast way to judge whether a full kriging workflow earns its cost. Bear in mind that the nugget you read off the plot is an extrapolation: the shortest lag bin has a positive centre distance, so the intercept at $h = 0$ is inferred from the fitted model, and a bin structure that starts too coarse will report as nugget what is really unresolved short-range structure.

Production Considerations #

The variogram cloud has $n(n-1)/2$ pairs, so pairwise-distance computation is $O(n^2)$ in both time and memory. At $n = 5{,}000$ that is ~12.5 million pairs — manageable — but at $n = 50{,}000$ the condensed distance array alone approaches 10 GB and naive estimation becomes infeasible. Three mitigations apply:

# 1. Subsample for exploratory estimation, then confirm on the full set.
rng = np.random.default_rng(42)
idx = rng.choice(values.size, size=3000, replace=False)
V_sample = Variogram(coordinates[idx], values[idx], n_lags=15, maxlag=0.5)

# 2. Cap maxlag so distant pairs are never materialised into bins.
# 3. For very large n, estimate per spatial tile and pool bin sums,
#    or move directional/large-scale estimation to gstools, which
#    supports binned estimators over chunked inputs.

Because the estimate is a sum over pairs, tiled estimation is exact if you accumulate the per-bin pair sums and counts across tiles rather than averaging tile-level curves. Keep the random seed fixed when subsampling so the exploratory variogram is reproducible.

Troubleshooting #

Next Steps #

Once the empirical curve is trustworthy, fit a theoretical model to it: the mechanics of weighted least-squares fitting in scikit-gstat are covered in fitting empirical variograms with scikit-gstat, and the binning decisions that shape the curve are detailed in choosing lag bins and bandwidth for variograms. From there, move to the sibling guides on theoretical variogram models and anisotropy and directional variograms.

Related

• Fitting Empirical Variograms with SciKit-GStat — weighted least-squares fitting of a model to the binned points
• Choosing Lag Bins and Bandwidth for Variograms — binning schemes, bin count, and maximum-lag selection
• Theoretical Variogram Models — the model families you fit to this empirical estimate
• Stationarity & Trend Analysis — the assumption a bounded variogram depends on

← Back to Variogram Modeling & Semivariance Analysis

Frequently Asked Questions #

What is the difference between the variogram cloud and the empirical variogram?

The variogram cloud is the raw scatter of one half-squared-difference per sample pair, plotted against each pair’s separation distance. The empirical variogram aggregates that cloud into lag bins and reports one averaged semivariance per bin. The cloud reveals outlier pairs and the spread within a lag; the binned curve is what you fit a theoretical model to.

How do I set the maximum lag distance in scikit-gstat?

Pass maxlag to the Variogram. Use a float for an absolute distance, a fraction between 0 and 1 to take that share of the maximum pairwise distance, or the string "median" to use the median pairwise distance. A robust default is about half the domain diameter, because bins beyond that contain too few pairs to estimate reliably.

Why are the last few points of my empirical variogram so noisy?

The widest lag bins contain the fewest pairs, because only observations on opposite sides of the study area fall into them. With few pairs the averaged semivariance has high variance and scatters wildly. Cap maxlag near half the domain diameter and require at least 30 pairs per bin so the unreliable tail is excluded from fitting.

Does scikit-gstat need projected coordinates?

Yes. The Variogram computes Euclidean lag distances directly from the coordinate array, so those coordinates must be in a projected metric CRS. Passing longitude and latitude in degrees produces lag distances in degrees, which correspond to different ground distances at different latitudes and distort the entire curve. Reproject with GeoPandas to a UTM zone first.

How many samples do I need for a reliable empirical variogram?

Around 100 well-distributed observations is the practical floor, because 15 lag bins each need 30 to 50 pairs and the pairs are not spread evenly across lags. Below roughly 50 points the fitted range commonly varies by a factor of two between random resamples of the same field, so the curve describes the sample rather than the process. Between 100 and 200 points the sill is usually resolved but the nugget stays poorly constrained, because the shortest lag holds the fewest pairs. Above about 500 points the binned curve stabilises and the choice of estimator matters more than the sample size.