Variogram fitting

Kriging, IDW, and SpatialFilter use a fitted variogram when no radius is provided. A variogram describes how dissimilar values become as the distance between points increases. Variogram builds it in two steps: estimate the pattern from the data, then fit a theoretical model to that estimate.

Empirical semivariogram

For every pair of source points i,ji, j separated by distance h=dijh = d_{ij}, the squared half-difference 12(zi−zj)2\frac{1}{2}(z_i - z_j)^2 is a noisy sample of how dissimilar values tend to be at that separation. Averaging these over all pairs whose distance falls in a lag bin around hh gives the classical (Matheron) empirical semivariogram estimator (Matheron, 1971; Cressie, 1993):

γ^(h)=12 ∣N(h)∣∑(i,j)∈N(h)(zi−zj)2\hat{\gamma}(h) = \frac{1}{2\,|N(h)|} \sum_{(i,j) \in N(h)} \left(z_i - z_j\right)^2

where N(h)N(h) is the set of point pairs whose separation falls in the lag bin centered on hh. n_lags controls how many such bins the pairwise distance range is divided into; max_pairs bounds how many point pairs are sampled when there are more than that available; max_distance bounds the maximum separation considered.

Theoretical models

The empirical estimate γ^(h)\hat\gamma(h) is a scatter of noisy points, one per lag bin, rather than a function that can be evaluated at any distance. Variogram fits three candidate theoretical models to it, each parameterized by a nugget c0c_0 (discontinuity at h=0h=0, representing measurement noise or sub-lag-scale variability), a partial sill cc (additional variance contributed by spatial structure), and a range aa (the distance beyond which points are treated as spatially uncorrelated). Writing r=h/ar = h / a:

Spherical

γ(h)={c0+c(1.5 r−0.5 r3)h<ac0+ch≥a\gamma(h) = \begin{cases} c_0 + c\left(1.5\,r - 0.5\,r^{3}\right) & h < a \\ c_0 + c & h \geq a \end{cases}

Exponential

γ(h)=c0+c(1−e−3r)\gamma(h) = c_0 + c\left(1 - e^{-3r}\right)

Gaussian

γ(h)=c0+c(1−e−3r2)\gamma(h) = c_0 + c\left(1 - e^{-3r^{2}}\right)

The factor of 3 in the exponential and gaussian models is the practical range convention. It rescales aa so the model reaches about 95% of the sill at h=ah = a, matching the spherical model's behavior of reaching the sill exactly at h=ah = a (Isaaks and Srivastava, 1989). This keeps range comparable across the three models.

Fitting and model selection

For each candidate model, c0c_0, cc, and aa are chosen to minimize the (weighted) squared error between the model curve γ(h)\gamma(h) and the binned empirical estimate γ^(h)\hat\gamma(h). Variogram.fit returns the spherical, exponential, or gaussian model with the lowest fitting error, so the caller does not need to choose one first.

Variogram.evaluate(distances) returns the fitted model's semivariance at arbitrary distances. These values feed directly into the kriging linear system.