The Variogram Explained: How Far Does a Measurement Reach

Problem statement

Every interpolation method rests on one assumption: nearby places are more alike than distant ones. The variogram is how you check that assumption and, if it holds, measure it.

It answers three questions with numbers rather than intuition:

• Is there spatial structure at all? If not, no interpolator can beat the mean.
• How far does one measurement tell you anything? That is the range, and it sets how far you can defensibly interpolate.
• How much of the variation is unexplainable at any distance? That is the nugget, and it caps how good any surface can be.

Fitted to 500 elevation samples over a 10 km window in Snowdonia:

nugget      0 m²
sill   44,096 m²
range   3,959 m

Quick answer

The variogram is half the mean squared difference between pairs of points, binned by the distance between them:

import numpy as np


def empirical_variogram(points, values, bin_width=250, max_lag=5000,
                        min_pairs=30):
    """Semivariance against separation distance."""
    d = np.sqrt(((points[:, None, :] - points[None, :, :]) ** 2).sum(-1))
    gamma = 0.5 * (values[:, None] - values[None, :]) ** 2

    upper = np.triu_indices(len(points), 1)          # each pair once
    d, gamma = d[upper], gamma[upper]

    edges = np.arange(0, max_lag + bin_width, bin_width)
    lags, semivariance, counts = [], [], []
    for lo, hi in zip(edges[:-1], edges[1:]):
        in_bin = (d >= lo) & (d < hi)
        if in_bin.sum() >= min_pairs:
            lags.append((lo + hi) / 2)
            semivariance.append(gamma[in_bin].mean())
            counts.append(int(in_bin.sum()))
    return np.array(lags), np.array(semivariance), np.array(counts)


lags, semivariance, counts = empirical_variogram(points, values)
print("  lag (m)   semivariance (m²)   pairs")
for i in (0, 1, 2, 4, 7, 17, 19):
    print(f"{lags[i]:9,.0f} {semivariance[i]:17,.0f} {counts[i]:10,}")

  lag (m)   semivariance (m²)   pairs
      125               954        254
      375             4,118        719
      625            10,637      1,059
    1,125            19,834      1,833
    1,875            27,842      2,727
    4,375            44,198      4,238
    4,875            45,546      4,375

Semivariance rising with distance is spatial structure. A flat line is not.

Empirical curve, nugget, sill and range, a model

1. Compute the empirical variogram

Take every pair of samples, compute half the squared difference in value, and average within distance bins. With 500 points that is 124,750 pairs — instant. With 50,000 points it is 1.2 billion, and you need a subsample or a spatial index.

Two choices matter:

• Bin width. Too narrow and each bin is noisy; too wide and the shape is smoothed away. Aim for at least 30 pairs per bin, ideally hundreds.
• Maximum lag. Do not go beyond about half the study area's diameter. Beyond that, few pairs exist and they are all edge-to-edge, which biases the estimate.

2. Read the three numbers off the curve

Nugget — the intercept as distance goes to zero. Two measurements at the same place should agree perfectly, so a non-zero nugget is measurement error plus variation at scales finer than your closest pair. It is the floor on how well any interpolator can do.

Sill — the plateau. Here the curve flattens: 44,198 m² at 4,375 m, 43,527 at 4,625 m, 45,546 at 4,875 m. The fitted sill was 44,096 m².

Textbooks say the sill should equal the variance of the data. Against the samples it does not — their variance is 38,896 m², so the sill is 13% higher. Against the whole DEM, whose variance is 42,804 m², the gap is 3%. That is not an error; it is what happens when the range is a large fraction of the study area, and it is explained below.

Range — the distance at which the curve flattens. Beyond it, two points are no more alike than two points chosen at random. The fitted range here was 3,959 m, inside the largest observed lag of 4,875 m.

3. Fit a model

Kriging needs a continuous function, not bins. The spherical model is the usual default:

def spherical(h, nugget, sill, rng):
    inside = nugget + sill * (1.5 * h / rng - 0.5 * (h / rng) ** 3)
    return np.where(h == 0, 0, np.where(h <= rng, inside, nugget + sill))

Fitted to the data above it gives a range of 3,959 m and a sill of 44,096 m². Refitted on 15 equal bins out to half the largest separation — the binning the other interpolation pages here use — the same samples give 3,762 m and 42,519 m² (Example 4). The range is a model parameter, and it moves with the binning by about 5%; quote it with the binning.

Always check the fitted range against your largest lag. If the range comes back beyond it, the optimiser placed the plateau somewhere the data cannot see, and the only honest reading is "the range is at least max_lag".

4. Use the range to set your interpolation radius

The range is the practical answer to "how far can I interpolate?" Beyond it, the nearest sample carries no information about the target, and the best estimate is the local mean.

That gives you a defensible max_distance for masking a surface, and a defensible search radius for IDW and kriging neighbourhoods.

5. Check for anisotropy before trusting one curve

A single variogram assumes structure is the same in every direction. Rainfall along a valley, contamination down a hydraulic gradient and dune fields are all directional.

Compute variograms in four direction bands (0°, 45°, 90°, 135°, with a ±22.5° tolerance). If the ranges differ by more than about 50%, the field is anisotropic and an isotropic model will over-smooth along the long axis and under-smooth across it.

Fitting, directional variograms and the cloud

The data used on this page

The surface is real: the Copernicus DEM GLO-30 over a 10 × 10 km window centred on the summit of Snowdon (Yr Wyddfa), warped to a 30 m grid in EPSG:32630 (x 422,880–432,870, y 5,875,440–5,885,430; 333 × 333 cells, 99.8 km²). The 500 samples are cell centres drawn at random with a fixed seed; 12,000 further cells, none of them samples, are held out as test_points/test_values. grid_points/truth are every cell centre and its elevation, and truth_fn reads the DEM bilinearly at any location, which the other interpolation guides score against. Here only the 500 samples (points, values) and the DEM itself are used.

import os

import numpy as np
import rasterio
from rasterio.enums import Resampling
from rasterio.transform import from_origin
from rasterio.warp import reproject
from scipy.interpolate import RegularGridInterpolator

bounds = (422_880.0, 5_875_440.0, 432_870.0, 5_885_430.0)
transform = from_origin(bounds[0], bounds[3], 30, 30)
shape = (333, 333)
if not os.path.exists("snowdon_dem_utm30.tif"):
    COP = ("https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N53_00_{0}_00_DEM/"
           "Copernicus_DSM_COG_10_N53_00_{0}_00_DEM.tif")
    dem = np.full(shape, np.nan, "float32")
    for tile in ("W005", "W004"):
        part = np.full(shape, np.nan, "float32")
        with rasterio.open(COP.format(tile)) as src:
            reproject(rasterio.band(src, 1), part, dst_transform=transform, dst_crs="EPSG:32630",
                      dst_nodata=np.nan, resampling=Resampling.bilinear)
        dem = np.where(np.isnan(dem), part, dem)
    with rasterio.open("snowdon_dem_utm30.tif", "w", driver="GTiff", width=shape[1], height=shape[0], count=1,
                       dtype="float32", crs="EPSG:32630", transform=transform, nodata=np.nan) as dst:
        dst.write(dem, 1)
dem = rasterio.open("snowdon_dem_utm30.tif").read(1)

xs = bounds[0] + (np.arange(shape[1]) + 0.5) * 30
ys = bounds[3] - (np.arange(shape[0]) + 0.5) * 30
gx, gy = np.meshgrid(xs, ys)
grid_points = np.column_stack([gx.ravel(), gy.ravel()])
truth = dem.ravel().astype("float64")

rng = np.random.default_rng(0)
order = rng.permutation(truth.size)
sample_idx, test_idx = order[:500], order[500:12_500]
points, values = grid_points[sample_idx], truth[sample_idx]
test_points, test_values = grid_points[test_idx], truth[test_idx]
_bilinear = RegularGridInterpolator((ys[::-1], xs), dem[::-1], bounds_error=False, fill_value=None)
truth_fn = lambda xy: _bilinear(np.column_stack([xy[:, 1], xy[:, 0]]))
print(f"DEM {np.nanmin(dem):.0f} to {np.nanmax(dem):.0f} m; {len(points)} samples, {len(test_points):,} held out")

DEM 52 to 1072 m; 500 samples, 12,000 held out

Example 1 — variogram, fit and sanity checks together

import numpy as np
from scipy.optimize import curve_fit


def spherical(h, nugget, sill, rng):
    h = np.asarray(h, dtype=float)
    inside = nugget + sill * (1.5 * h / rng - 0.5 * (h / rng) ** 3)
    return np.where(h == 0, 0.0, np.where(h <= rng, inside, nugget + sill))


def fit_variogram(lags, semivariance, values, extent, model=spherical):
    """Fit, then say plainly whether the fit is supported by the data.

    extent: the diameter of the study area, used to flag the case where the
    range is large enough that the sample variance stops being a fair check.
    """
    p0 = [semivariance[0] * 0.1, semivariance.max(), lags.max() / 2]
    bounds = ([0, 0, lags[0]], [semivariance.max(), semivariance.max() * 3,
                                lags.max() * 3])
    params, _ = curve_fit(model, lags, semivariance, p0=p0,
                          bounds=bounds, maxfev=20000)
    nugget, sill, rng = params

    data_variance = float(np.var(values))
    print(f"  nugget {nugget:10,.0f}   ({nugget / (nugget + sill):.1%} of total)")
    print(f"  sill   {sill:10,.0f}   data variance {data_variance:,.0f} "
          f"({sill / data_variance - 1:+.0%} vs data)")
    print(f"  range  {rng:10,.0f} m  last observed lag {lags.max():,.0f} m")

    if rng > lags.max():
        print("  ! the fitted range is beyond the largest lag — treat it as "
              "a lower bound and widen max_lag if the extent allows")
    if rng > 0.25 * extent:
        print("  ! the range is a large fraction of the study extent — the "
              "sample variance underestimates the sill, so do not use "
              "'sill == variance' as a check here")
    return params


extent = float(np.hypot(bounds[2] - bounds[0], bounds[3] - bounds[1]))
nugget, sill, rng = fit_variogram(lags, semivariance, values, extent)

  nugget          0   (0.0% of total)
  sill       44,096   data variance 38,896 (+13% vs data)
  range       3,959 m  last observed lag 4,875 m
  ! the range is a large fraction of the study extent — the sample variance underestimates the sill, so do not use 'sill == variance' as a check here

Printing the warning is the point of the function. A fit that silently returns three numbers invites them to be quoted as measurements.

Example 2 — directional variograms

import numpy as np


def directional_variogram(points, values, angle_deg, tolerance_deg=22.5,
                          bin_width=250, max_lag=5000):
    """The variogram restricted to pairs aligned with one direction."""
    diff = points[:, None, :] - points[None, :, :]
    d = np.sqrt((diff ** 2).sum(-1))
    gamma = 0.5 * (values[:, None] - values[None, :]) ** 2

    with np.errstate(invalid="ignore"):
        bearing = np.degrees(np.arctan2(diff[..., 1], diff[..., 0])) % 180
    target = angle_deg % 180
    offset = np.minimum(np.abs(bearing - target), 180 - np.abs(bearing - target))

    upper = np.triu_indices(len(points), 1)
    d, gamma, offset = d[upper], gamma[upper], offset[upper]
    aligned = offset <= tolerance_deg

    edges = np.arange(0, max_lag + bin_width, bin_width)
    out = []
    for lo, hi in zip(edges[:-1], edges[1:]):
        m = aligned & (d >= lo) & (d < hi)
        if m.sum() >= 30:
            out.append((float((lo + hi) / 2), float(gamma[m].mean()), int(m.sum())))
    return out

Run it for 0, 45, 90 and 135 degrees and compare where each curve flattens. Different ranges in different directions means anisotropy, which most kriging implementations can model once you give them the major axis and the ratio.

Example 3 — a variogram cloud, for finding the outlier that ruined the fit

import numpy as np


def variogram_cloud(points, values, max_lag=2000, top=10):
    """Every pair as a point, so a single bad sample becomes visible."""
    d = np.sqrt(((points[:, None, :] - points[None, :, :]) ** 2).sum(-1))
    gamma = 0.5 * (values[:, None] - values[None, :]) ** 2
    i, j = np.triu_indices(len(points), 1)
    d, gamma = d[i, j], gamma[i, j]

    close = d < max_lag
    order = np.argsort(-gamma[close])[:top]
    idx_i, idx_j = i[close][order], j[close][order]

    print(f"  {top} most discordant close pairs:")
    for a, b in zip(idx_i, idx_j):
        print(f"    samples {a:4d} and {b:4d}: {np.hypot(*(points[a] - points[b])):6.0f} m "
              f"apart, values {values[a]:8.1f} and {values[b]:8.1f}")

    from collections import Counter
    culprits = Counter(np.concatenate([idx_i, idx_j]).tolist())
    print(f"  samples appearing most often: {culprits.most_common(3)}")
    return d[close], gamma[close]

A single sample with a transcription error — a decimal point in the wrong place — appears in many of the most discordant close pairs. The binned variogram hides it in an average; the cloud does not.

Example 4 — the other numbers on this page

print(f"pairs: {len(points) * (len(points) - 1) // 2:,}")
print("plateau:", ", ".join(f"{g:,.0f} m² at {l:,.0f} m" for l, g in zip(lags[-3:], semivariance[-3:])))
print(f"variance: samples {values.var():,.0f} m², every DEM cell {np.nanvar(dem):,.0f} m²; "
      f"sill / sample variance {sill / values.var():.2f}")

# the same fit on 15 equal bins out to half the largest separation, as the other interpolation pages do
from scipy.spatial.distance import pdist
d, g = pdist(points), 0.5 * pdist(values[:, None], "sqeuclidean")
edges = np.linspace(0, d.max() / 2, 16)
lag15 = np.array([d[(d >= a) & (d < b)].mean() for a, b in zip(edges[:-1], edges[1:])])
sv15 = np.array([g[(d >= a) & (d < b)].mean() for a, b in zip(edges[:-1], edges[1:])])
print("15 equal bins:")
_ = fit_variogram(lag15, sv15, values, extent)

# does a regional trend explain the high sill?
A = np.column_stack([np.ones(len(points)), (points - points.mean(0)) / 1000])
coef, *_ = np.linalg.lstsq(A, values, rcond=None)
residual = values - A @ coef
print(f"plane: {coef[1]:+.1f} m/km east, {coef[2]:+.1f} m/km north, "
      f"explains {1 - residual.var() / values.var():.1%} of the variance")
r_lags, r_sv, _ = empirical_variogram(points, residual)
print("detrended:")
_ = fit_variogram(r_lags, r_sv, residual, extent)

pairs: 124,750
plateau: 44,198 m² at 4,375 m, 43,527 m² at 4,625 m, 45,546 m² at 4,875 m
variance: samples 38,896 m², every DEM cell 42,804 m²; sill / sample variance 1.13
15 equal bins:
  nugget          0   (0.0% of total)
  sill       42,519   data variance 38,896 (+9% vs data)
  range       3,762 m  last observed lag 6,569 m
  ! the range is a large fraction of the study extent — the sample variance underestimates the sill, so do not use 'sill == variance' as a check here
plane: +8.0 m/km east, +26.1 m/km north, explains 16.6% of the variance
detrended:
  nugget          0   (0.0% of total)
  sill       40,952   data variance 32,455 (+26% vs data)
  range       3,690 m  last observed lag 4,875 m
  ! the range is a large fraction of the study extent — the sample variance underestimates the sill, so do not use 'sill == variance' as a check here

Semivariance, the nugget and a shifted sill

Why half the squared difference

The factor of one half makes the semivariance equal the variance of the difference between two points divided by two, which under the standard stationarity assumption equals C(0) − C(h) — the covariance at zero lag minus the covariance at lag h.

That identity is what lets kriging use the variogram directly in its weighting system. It also means the sill should equal the variance of the data, which is the check in Example 1.

Why the nugget is not always error

A nugget of zero, as measured here, says elevation is spatially continuous down to the closest sample separation, which is what you expect from a DEM: it is itself a smooth interpolated product with no measurement noise at 30 m.

A large nugget in field data means one of two things, and the variogram cannot tell them apart:

• Measurement error — repeat samples at the same spot would disagree.
• Micro-scale variation — real variation at distances shorter than your closest pair.

The distinction matters. Measurement error should be smoothed through; micro-scale variation is real signal you have failed to sample. The way to separate them is to add a few duplicate or very-close samples to the design.

Why the sill does not equal the data variance here

"The sill should equal the variance of the data" is the standard check, and these samples fail it by 13%: sill 44,096 m² against a sample variance of 38,896 m². The variance of every cell in the window is 42,804 m², within 3% of the sill — so most of the gap is 500 samples under-representing the window's spread.

The check is not wrong; its precondition is. It holds when the range is small compared with the study area, so that the domain contains many independent patches. Here the range is 4.0 km inside a 10 km window — the domain holds about two independent patches across, and the sample variance of such a domain systematically underestimates the true sill, because the sample never sees the full spread the process is capable of.

It is worth ruling out the obvious alternative. A regional trend also inflates a variogram, so fitting on detrended residuals should fix it if a trend is to blame. Removing a plane (+8.0 m/km east, +26.1 m/km north) explained only 16.6% of the variance and made the ratio worse, 1.13 to 1.26. The trend is not the cause; the domain size is.

The practical rule: use the sill-equals-variance check only when the fitted range is under about a quarter of the study extent. Above that, expect the sill to exceed the variance and do not "correct" it.

Why a shifted sill matters for kriging variance

Kriging predictions are insensitive to a scaled sill — multiply the whole variogram by a constant and the weights are unchanged. Kriging variances scale directly with it.

Kriging variances are also set mostly by the model's shape at short lags, which is where each prediction's nearest samples sit. In My kriging fit fails or returns a flat surface the kriging variance overstated the true error by 64% on these samples, because the fitted curve sat 1.8 times above the empirical semivariance at the first lag. Use kriging variance for relative comparison across a map; do not read it as a confidence interval unless the variogram is well constrained by lags well short of the domain size.

Trends, pairs per bin and flat variograms

• Use a projected CRS. Lag distances in degrees are not distances.
• Trend first, variogram second. A strong regional trend makes the variogram rise without limit. Remove the trend, fit on the residuals, add the trend back.
• At least 30 pairs per bin, and ideally 100+. Noisy bins produce arbitrary fits.
• Do not exceed half the study extent as a maximum lag.
• The number of pairs grows as n². Above about 5,000 samples, subsample or bin spatially.
• A flat variogram means no spatial structure. Stop; interpolation cannot help.
• A sill above the data variance is expected when the range approaches the domain size. Removing a trend will not fix it.
• A sill far below the data variance usually means a trend was removed too aggressively.
• Anisotropy is common and invisible in an omnidirectional variogram.

Internal links

• Spatial interpolation explained — what the range is for
• How to fit a variogram and krige a surface in Python — using the fitted model
• My kriging fit fails or returns a flat surface — when the fit goes wrong
• IDW, kriging, splines or TIN? Choosing an interpolator — when the variogram is worth the effort
• Spatial autocorrelation explained — the same structure measured a different way
• Sample design explained: where to measure — using the range to plan a survey
• How to choose a cell size for an interpolated surface — the range as an upper bound on useful detail
• How to cross-validate an interpolated surface — validating the model the variogram produced

FAQ

What is a variogram?

A plot of half the mean squared difference between pairs of samples against the distance separating them. It measures how quickly similarity decays with distance.

What are the nugget, sill and range?

The nugget is the intercept — variation at zero distance, from measurement error or micro-scale structure. The sill is the plateau, which should match the data variance. The range is the distance at which the plateau is reached.

What does a flat variogram mean?

No spatial structure. Nearby samples are no more alike than distant ones, so no interpolator will beat predicting the mean everywhere.

How far can I interpolate?

Not beyond the range. Past it, the nearest sample carries no information about the target location.

Why is my fitted sill larger than my data variance?

Usually because the range is a large fraction of the study area. Here the range was 4.0 km inside a 10 km window and the sill came out 13% above the sample variance; detrending made it worse (1.13 to 1.26), not better, so a trend is not the cause.

How many samples do I need for a variogram?

At least 100 for a rough shape, 200–500 for a usable fit, and enough that every distance bin holds 30 or more pairs.

Do I need a variogram for IDW?

Not to run it, but yes to justify it. The range tells you what search radius is defensible, and a flat variogram tells you not to bother at all.