Fitting a Nested Anisotropic Variogram Model

TL;DR: When one model family cannot span both scales, sum several: nugget + sph_reduced(u1, c1) + exp_reduced(u2, c2), where each u_k is the lag re-expressed in that structure’s own range, anisotropy ratio and axis angle. Fit it in stages, refine jointly with scipy.optimize.least_squares, prove admissibility with np.linalg.cholesky, and keep the extra parameters only if leave-one-out kriging improves.

Why This Matters #

A soil, groundwater or ore-grade field rarely has a single characteristic length. Local variation from sampling support, micro-topography and short-range mixing produces a steep rise over the first hundred metres or so; a regional gradient — a plume moving with the hydraulic head, a stratigraphic dip, a prevailing wind — produces a much slower climb that is strongly directional. A single spherical or exponential model has one range and one shape parameter, so fitting it to both leaves the near lags underestimated and the far lags overestimated at the same time. The kriging weights that follow are wrong in a systematic, direction-dependent way, and the kriging variance is confidently wrong, which is worse.

The fix is a nested model: a sum of admissible structures, each with its own sill, range and anisotropy. This page assumes you can already read a set of directional variograms and have seen geometric anisotropy detected and modelled in Python for a single structure; it extends that machinery to two. The single-structure ground work — spherical, exponential and Gaussian shapes, and how their range parameters are defined — is covered under theoretical variogram models.

Formally, write the nested variogram as

$$\gamma(\mathbf{h}) = c_0\,\mathbb{1}[\mathbf{h} \neq \mathbf{0}] + \sum_{k=1}^{K} c_k\, g_k\!\big(u_k(\mathbf{h})\big),$$

where $c_0$ is the nugget, $c_k \geq 0$ is the sill contribution of structure $k$, $g_k$ is a normalised admissible shape (spherical, exponential, Matérn) rising from 0 to 1, and $u_k$ is the lag vector re-expressed in that structure’s own metric:

$$u_k(\mathbf{h}) = \frac{1}{a_k}\sqrt{\big(h_x\cos\theta_k + h_y\sin\theta_k\big)^2 + \frac{\big(-h_x\sin\theta_k + h_y\cos\theta_k\big)^2}{\lambda_k^{2}}}.$$

Here $a_k$ is the range along the axis of maximum continuity, $\theta_k$ that axis’s azimuth, and $\lambda_k \in (0, 1]$ the ratio of minor to major range. Admissibility is inherited: each $c_k g_k(u_k)$ is a valid variogram because a non-singular linear transform of a positive-definite isotropic covariance is still positive-definite, and a non-negative sum of conditionally negative-definite functions is conditionally negative-definite. The one thing you may not do is let any $c_k$ go negative to buy a better fit.

Environment and Version Pinning #

pip install "numpy==1.26.4" "scipy==1.13.1" "gstools==1.5.2" \
            "scikit-gstat>=1.0,<2.0" "matplotlib==3.9.2"

import numpy as np
import gstools as gs
from scipy.optimize import least_squares

rng = np.random.default_rng(20260807)

Step-by-Step Implementation #

1. Build a field that genuinely has two structures #

Simulating the field you intend to recover is the only way to know whether the fitting procedure works. Two independent random fields are added, so the resulting covariance really is the sum of the two models plus a nugget.

N, SIDE = 600, 6000.0
x = rng.uniform(0.0, SIDE, N)
y = rng.uniform(0.0, SIDE, N)

# gstools: len_scale IS the range for Spherical; for Exponential the
# practical range (95% of sill) is about 3 * len_scale.
short = gs.Spherical(dim=2, var=0.55, len_scale=180.0)
longr = gs.Exponential(dim=2, var=1.30, len_scale=900.0,
                       anis=0.30, angles=np.deg2rad(35.0))

z = (gs.SRF(short, seed=11)((x, y))
     + gs.SRF(longr, seed=29)((x, y))
     + rng.normal(0.0, np.sqrt(0.15), N))

print(f"n = {N}, sample variance = {z.var(ddof=1):.3f}")

n = 600, sample variance = 1.938

The theoretical total sill is $0.15 + 0.55 + 1.30 = 2.00$, so a sample variance of 1.938 is about right for 600 points over six kilometres.

2. Compute directional empirical variograms #

Writing the estimator out makes the azimuth convention unambiguous, which matters more than it sounds — see the best practice on convention drift below.

def directional_variogram(x, y, z, azimuth_deg, tol_deg,
                          n_lags=18, maxlag=3000.0, min_pairs=30):
    """Matheron estimator restricted to a symmetric azimuth cone.

    Azimuth is measured anticlockwise from the +x axis, in degrees.
    A cone of 90 degrees gives the omnidirectional variogram.
    """
    P = np.column_stack([x, y])
    i, j = np.triu_indices(len(z), k=1)
    d = P[j] - P[i]
    h = np.hypot(d[:, 0], d[:, 1])
    ang = np.degrees(np.arctan2(d[:, 1], d[:, 0])) % 180.0
    dev = np.abs(ang - (azimuth_deg % 180.0))
    dev = np.minimum(dev, 180.0 - dev)

    keep = (dev <= tol_deg) & (h <= maxlag) & (h > 0.0)
    edges = np.linspace(0.0, maxlag, n_lags + 1)
    idx = np.digitize(h[keep], edges) - 1
    sq = (z[i[keep]] - z[j[keep]]) ** 2
    hk = h[keep]

    bins, gam, cnt = [], [], []
    for b in range(n_lags):
        m = idx == b
        if m.sum() >= min_pairs:
            bins.append(hk[m].mean())
            gam.append(0.5 * sq[m].mean())
            cnt.append(float(m.sum()))
    return np.array(bins), np.array(gam), np.array(cnt)


omni = directional_variogram(x, y, z, 0.0, 90.0)
major = directional_variogram(x, y, z, 35.0, 22.5)
minor = directional_variogram(x, y, z, 125.0, 22.5)

for tag, v in (("major", major), ("minor", minor)):
    k = np.argmin(np.abs(v[0] - 750.0))
    print(f"{tag}: gamma({v[0][k]:.0f} m) = {v[1][k]:.2f}, "
          f"gamma({v[0][-1]:.0f} m) = {v[1][-1]:.2f}")

major: gamma(750 m) = 1.44, gamma(2917 m) = 1.95
minor: gamma(750 m) = 1.89, gamma(2917 m) = 2.02

Both directions reach the same sill near 2.0, and they differ only in how fast they get there — the signature of geometric rather than zonal anisotropy. What the numbers do not show is the shape of the climb, and that is where the case for a second structure is made.

3. Fit the short structure to the near lags #

Restrict the first fit to lags the short structure can actually see. The long structure has already contributed something by 400 m, so this stage is biased — deliberately, and only as a starting point.

def sph_reduced(u, sill):
    t = np.clip(u, 0.0, 1.0)
    return sill * (1.5 * t - 0.5 * t ** 3)

def exp_reduced(u, sill):
    return sill * (1.0 - np.exp(-u))

def reduced_lag(h, az_deg, a, ratio, theta_deg):
    """Lag h at azimuth az_deg, expressed in one structure's own metric."""
    phi = np.deg2rad(az_deg - theta_deg)
    return h * np.hypot(np.cos(phi), np.sin(phi) / ratio) / a


ob, og, oc = omni
near = ob <= 400.0

def resid_short(p):
    nug, c1, a1 = p
    pred = nug + sph_reduced(ob[near] / a1, c1)
    return np.sqrt(oc[near]) * (og[near] - pred) / pred

fit1 = least_squares(resid_short, [0.2, 0.6, 200.0],
                     bounds=([0.0, 0.0, 20.0], [1.0, 2.0, 600.0]))
nug0, c1_0, a1_0 = fit1.x
print(f"stage 1  nugget={nug0:.3f}  sill={c1_0:.3f}  range={a1_0:.0f} m")

stage 1  nugget=0.155  sill=0.612  range=191 m

The sill comes out at 0.612 against a true 0.55, and the range at 191 m against 180 m, precisely because the long structure’s early contribution has been absorbed. That is expected and harmless.

4. Subtract, then fit the long anisotropic structure #

Now use every azimuth at once. Subtracting the stage-one model from each directional empirical point leaves a residual that the long structure alone should explain, and fitting all cones jointly is what identifies its ratio and angle.

AZ = np.arange(0.0, 180.0, 22.5)
emp = [(a, *directional_variogram(x, y, z, a, 22.5)) for a in AZ]

def resid_long(p):
    c2, ell2, ratio2, th2 = p
    out = []
    for az, b, g, c in emp:
        target = g - (nug0 + sph_reduced(b / a1_0, c1_0))
        pred = exp_reduced(reduced_lag(b, az, ell2, ratio2, th2), c2)
        out.append(np.sqrt(c) * (target - pred))
    return np.concatenate(out)

fit2 = least_squares(resid_long, [1.0, 800.0, 0.5, 30.0],
                     bounds=([0.0, 100.0, 0.05, -90.0],
                             [3.0, 4000.0, 1.00, 270.0]))
c2_0, ell2_0, r2_0, th2_0 = fit2.x
print(f"stage 2  sill={c2_0:.3f}  len_scale={ell2_0:.0f} m  "
      f"ratio={r2_0:.2f}  angle={th2_0 % 180.0:.1f} deg")

stage 2  sill=1.221  len_scale=1032 m  ratio=0.31  angle=32.8 deg

Note the plain sqrt(count) weighting here rather than the relative weighting of stage 3: the residual target passes through zero at short lag, so dividing by the model would blow up.

5. Attach the anisotropy to the structure, not to the model #

Each structure now carries four of its own parameters. reduced_lag is called twice per evaluation, once per structure, with different a, ratio and theta. That is the whole mechanism, and it is what a single global coordinate rotation cannot reproduce.

6. Refine all nine parameters jointly #

def nested_gamma(h, az, p):
    nug, c1, a1, r1, t1, c2, ell2, r2, t2 = p
    return (nug
            + sph_reduced(reduced_lag(h, az, a1, r1, t1), c1)
            + exp_reduced(reduced_lag(h, az, ell2, r2, t2), c2))

def resid_joint(p):
    out = []
    for az, b, g, c in emp:
        pred = nested_gamma(b, az, p)
        out.append(np.sqrt(c) * (g - pred) / pred)
    return np.concatenate(out)

p_start = [nug0, c1_0, a1_0, 1.0, 0.0, c2_0, ell2_0, r2_0, th2_0]
lo = [0.0, 0.0,  20.0, 0.10, -90.0, 0.0,  100.0, 0.05, -90.0]
hi = [1.0, 2.0, 900.0, 1.00, 270.0, 3.0, 5000.0, 1.00, 270.0]

fit3 = least_squares(resid_joint, p_start, bounds=(lo, hi), x_scale="jac")
p_hat = fit3.x
names = ["nugget", "c1", "a1", "ratio1", "angle1",
         "c2", "len2", "ratio2", "angle2"]
for n_, v in zip(names, p_hat):
    print(f"{n_:>7} = {v:9.3f}")
print(f"weighted SSR = {2 * fit3.cost:.2f}")

 nugget =     0.142
     c1 =     0.571
     a1 =   168.412
 ratio1 =     0.940
 angle1 =    18.700
     c2 =     1.284
   len2 =   947.106
 ratio2 =     0.281
 angle2 =    33.407
weighted SSR = 214.63

7. Prove the model is admissible #

Admissibility is guaranteed by construction only while every sill stays non-negative, so test it rather than assume it. Build the covariance matrix on a dense set of locations and try to factorise it; a model that is not positive-definite will fail. Note the explicit handling of the origin — the nugget makes $\gamma$ discontinuous there, and $C(\mathbf{0})$ is the full sill.

def nested_cov(dx, dy, p):
    h = np.hypot(dx, dy)
    az = np.degrees(np.arctan2(dy, dx))
    total = p[0] + p[1] + p[5]
    return np.where(h == 0.0, total, total - nested_gamma(h, az, p))

def is_admissible(p, cov_fn=None, n_test=400, seed=0, jitter=1e-9):
    cov_fn = cov_fn or nested_cov
    g = np.random.default_rng(seed)
    px, py = g.uniform(0, SIDE, n_test), g.uniform(0, SIDE, n_test)
    C = cov_fn(px[:, None] - px[None, :],
               py[:, None] - py[None, :], p) + jitter * np.eye(n_test)
    try:
        np.linalg.cholesky(C)
        return True
    except np.linalg.LinAlgError:
        return False

bad = p_hat.copy()
bad[5] = -0.40                      # negative sill on the long structure
print("fitted model admissible :", is_admissible(p_hat))
print("negative-sill admissible:", is_admissible(bad))

fitted model admissible : True
negative-sill admissible: False

8. Cross-validate against the single-structure fit #

Four extra parameters have to earn their place. Fit the best single anisotropic structure to the same empirical variograms, then run leave-one-out ordinary kriging with both models.

def single_gamma(h, az, p):
    nug, c, a, r, t = p
    return nug + sph_reduced(reduced_lag(h, az, a, r, t), c)

def single_cov(dx, dy, p):
    h = np.hypot(dx, dy)
    az = np.degrees(np.arctan2(dy, dx))
    total = p[0] + p[1]
    return np.where(h == 0.0, total, total - single_gamma(h, az, p))

fit_s = least_squares(
    lambda p: np.concatenate([np.sqrt(c) * (g - single_gamma(b, az, p))
                              / single_gamma(b, az, p)
                              for az, b, g, c in emp]),
    [0.3, 1.6, 1200.0, 0.4, 33.0],
    bounds=([0.0, 0.0, 50.0, 0.05, -90.0], [1.5, 4.0, 6000.0, 1.0, 270.0]))
print("single-structure fit:", np.round(fit_s.x, 3))

def loo_kriging(x, y, z, cov_fn, p, total_sill):
    n = len(z)
    C = cov_fn(x[:, None] - x[None, :], y[:, None] - y[None, :], p)
    pred, var = np.empty(n), np.empty(n)
    idx = np.arange(n)
    for k in range(n):
        m = idx != k
        A = np.ones((n, n))
        A[:n - 1, :n - 1] = C[np.ix_(m, m)]
        A[n - 1, n - 1] = 0.0
        b = np.ones(n)
        b[:n - 1] = C[m, k]
        sol = np.linalg.solve(A, b)
        w, mu = sol[:n - 1], sol[n - 1]
        pred[k] = w @ z[m]
        var[k] = total_sill - w @ C[m, k] - mu
    return pred, var

for tag, cov, par in (("single", single_cov, fit_s.x),
                      ("nested", nested_cov, p_hat)):
    tot = par[0] + par[1] + (par[5] if len(par) == 9 else 0.0)
    pr, vr = loo_kriging(x, y, z, cov, par, tot)
    e = z - pr
    print(f"{tag:>7}  RMSE={np.sqrt((e ** 2).mean()):.3f}  "
          f"MAE={np.abs(e).mean():.3f}  MSSE={(e ** 2 / vr).mean():.3f}")

single-structure fit: [   0.31     1.617 1449.6      0.331   34.113]
 single  RMSE=0.612  MAE=0.472  MSSE=1.412
 nested  RMSE=0.548  MAE=0.421  MSSE=1.028

Interpreting the Output #

The nine refined parameters recover the simulation almost exactly: nugget 0.142 against 0.15, spherical sill 0.571 and range 168 m against 0.55 and 180 m, exponential sill 1.284 and length scale 947 m against 1.30 and 900 m, ratio 0.281 against 0.30 and axis 33.4° against 35°. Note how much the joint refinement moved stage 1: the spherical sill fell from 0.612 to 0.571 once the long structure was allowed to claim its share of the near lags.

Two fitted values should be read as absences rather than estimates. ratio1 = 0.940 is close enough to 1 that the short structure is isotropic, and angle1 = 18.7 is then meaningless — the objective is flat in that direction, and re-running with a different seed will return a wildly different angle with no change in fit. That is not a failure; it is the model telling you two parameters are unidentifiable and should be removed.

The cross-validation is where the extra parameters are judged. RMSE falls by 10.5 per cent, which is worthwhile but not dramatic. The mean standardised squared error is the more important number: it should sit near 1 if the kriging variances are calibrated, and the single-structure model’s 1.412 says its predictions were roughly 19 per cent more variable than it claimed. The nested model’s 1.028 is honest. Warning signs run the other way — a nested fit that improves RMSE by less than about 3 per cent, or that pushes the standardised error below about 0.8, is over-parameterised and is smoothing rather than modelling.

Critical Best Practices #

One structure per visible change of slope, and stop at three #

The discipline that stops a nested model becoming a spline is simple: add a structure only when the empirical variogram shows a change of slope you can point at, and only when the new structure’s range lands between the shortest reliable lag and half the study extent. Anything with a range shorter than your minimum sample spacing is indistinguishable from nugget; anything longer than half the extent is indistinguishable from trend, and belongs in the drift, not the variogram. Three structures is the practical ceiling for field data.

Fit the short structure on the lags it can see, and expect bias #

Fitting the near lags in isolation always attributes some of the long structure’s early rise to the short one — 0.612 rather than 0.55 here. That is why stage 1 is a starting point rather than a result, and why the joint refinement in stage 6 is not optional. Skipping the refinement leaves the short structure’s sill overstated by around ten per cent and the total sill correspondingly wrong.

Drop the anisotropy parameters a structure cannot support #

A fitted ratio above roughly 0.85 means the structure is isotropic, and its angle is then pure noise from a flat objective. Refit with ratio1 fixed at 1.0 and angle1 removed: the fit will be statistically indistinguishable, the parameter count falls from nine to seven, and the covariance of the remaining estimates becomes far better conditioned. Reporting an angle you cannot identify invites someone downstream to believe it.

Verify the azimuth convention against a field you built yourself #

Nothing here bites harder. scikit-gstat’s azimuth, gstools’ angles and pykrige’s anisotropy_angle do not all measure from the same zero, and not all measure in the same rotational direction — some take degrees, some radians. A 90° error looks exactly like an anisotropy you have got backwards. The only safe check is the one in step 1: simulate a field with a known axis, recover it, and confirm the number you get back is the number you put in. The same discipline applies when reading fitted spherical, exponential and Gaussian variogram models, where len_scale means the range for one family and a third of the practical range for another.

A single global anisotropy transform cannot carry a nested model #

pykrige.ok.OrdinaryKriging accepts anisotropy_scaling and anisotropy_angle, but it applies that transform to the coordinates before any variogram is evaluated — one ellipse for the whole model. If your two structures have genuinely different ratios or angles, that route cannot represent them, no matter which variogram family you choose. Either build the covariance yourself, as above, or use gstools models, which carry anis and angles per model instance so a sum of instances keeps a sum of transforms.

Troubleshooting #

Next Steps #

Take the fitted nine parameters into cross-validating a variogram model in Python for the full diagnostic battery — standardised error distributions, spatial maps of residuals and $k$-fold alternatives to leave-one-out — before committing the model to a kriging run.

Frequently Asked Questions #

How many structures should a nested variogram have?

One per visible change of slope in the empirical variogram, and in practice almost never more than three. Each structure costs four parameters once it carries its own anisotropy, so a third structure fitted to twenty lag bins per azimuth is fitting noise. The working rule is that a structure must be identifiable from the data you have: its range must fall between the shortest reliable lag and half the study extent, and its sill must exceed about five per cent of the total.

Can each structure really have a different anisotropy angle?

Yes. Each structure is an isotropic covariance evaluated on its own linearly transformed lag vector, and a linear transform of a positive-definite function is still positive-definite. Summing several such terms with non-negative sills therefore yields an admissible model regardless of whether the transforms agree. This is exactly what a single global anisotropy rotation cannot express, which is why libraries that apply one coordinate transform to the whole model cannot represent a genuinely multi-scale anisotropic field.

Why fit in stages instead of optimising all parameters at once?

Because a nine-parameter objective on a binned variogram has many shallow local minima, and a cold start typically collapses the two structures onto each other or drives one sill to zero. Staging gives each structure a starting value in the right order of magnitude: the short structure from the near lags, the long structure from the residual. The joint refinement that follows then moves parameters by tens of per cent rather than orders of magnitude, and it converges reliably.

How do I know the extra parameters are worth it?

Cross-validate. Fit the single-structure model on the same empirical variogram, run leave-one-out ordinary kriging with both, and compare prediction error alongside the mean standardised squared error. On the worked example the nested model cut leave-one-out RMSE from 0.612 to 0.548 and moved the standardised squared error from 1.412 to 1.028. The second number matters more: it says the nested model’s uncertainty estimates are honest where the single-structure model was overconfident.

Related

• Detecting and Modeling Geometric Anisotropy in Python — the single-structure case this page extends
• Fitting Spherical, Exponential & Gaussian Variogram Models — the admissible shapes each nested structure is built from
• Cross-Validating a Variogram Model in Python — the full diagnostic battery behind the comparison above

← Back to Anisotropy & Directional Variograms