🧭 Concept Interpolation & surfaces Spatial statistics Interpolation & Surfaces · step 5 of 13 Updated

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   51,039 m²
range   4,537 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)
  lag (m)   semivariance (m²)   pairs
      125             1,021        217
      375             4,590        692
      625             8,945      1,036
    1,125            19,938      1,813
    1,875            30,584      2,691
    4,375            50,987      4,263
    4,875            50,984      4,288

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

A variogram curve with the nugget at the intercept, the sill as the plateau and the range as the distance at which the plateau is reached.
Three numbers describe the curve, and each one answers a different practical question.

Step-by-step solution

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 cleanly: 50,987 m² at 4,375 m, 51,741 at 4,625 m, 50,984 at 4,875 m. The fitted sill was 51,039 m².

Textbooks say the sill should equal the variance of the data. This one does not — the DEM's variance is 41,950 m², so the sill is 22% higher. 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 4,537 m, comfortably 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 4,537 m and a sill of 51,039 m², and the empirical points sit on the curve out to the largest lag.

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.

A variogram plateau at 51,039 square metres against a data variance of 41,950, with the range at 4.5 km against a 10 km study window.
A range that is nearly half the study extent makes the sample variance a biased estimate of the sill. The check "sill should equal the variance" does not apply here.

Code examples

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:+.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
  nugget          0   (0.0% of total)
  sill       51,039   data variance 41,950 (+22% vs data)
  range       4,537 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.

Explanation

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 this dataset fails it by 22%: sill 51,039 m² against a variance of 41,950 m².

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.5 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 (+14.9 m/km east, +17.3 m/km north) explained only 10.7% of the variance and made the ratio slightly worse, 1.22 to 1.32. 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.

That is part of why kriging variance overstated the true error by 68% in How to fit a variogram and krige a surface. 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.

Five checks on a fitted variogram, ending with whether the range is small enough for the sill-equals-variance check to apply.
The last one is why this dataset fails the textbook check without anything being wrong.

Edge cases or notes

  • 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.

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.5 km inside a 10 km window and the sill came out 22% above the variance; detrending made it worse, not better, so the domain size is 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.