Cokriging & Multivariate Interpolation

Cokriging estimates an undersampled variable using its own data plus the data of one or more correlated variables that are cheaper to measure and therefore denser. The gain is real but conditional: it comes from the secondary filling gaps between primary samples, so it evaporates as the primary becomes dense, and it depends on a joint variogram model that cannot be fitted one curve at a time. This page, part of Kriging, Interpolation & Surface Generation Techniques, sets out the cross-variogram, the linear model of coregionalisation and its positive semi-definiteness requirement, the ordinary cokriging system, the collocated simplification, and the screening effect that decides whether any of it is worth doing.

Prerequisites #

• Python 3.10+
• gstools>=1.5, numpy>=1.24, scipy>=1.10; geopandas>=1.0 only if your data start as vector files
• A projected (metric) CRS, because every range and lag distance below is in metres
• A primary variable sampled at $n_1$ points and a secondary sampled at $n_2 \gg n_1$ points, with a collocated subset where both are measured
• Familiarity with fitting a single variogram — see theoretical variogram models for the permissible model families

Mathematical Core #

Write the two variables as second-order stationary random functions $Z_1(\mathbf{s})$ (primary, expensive, sparse) and $Z_2(\mathbf{s})$ (secondary, cheap, dense), with means $m_1, m_2$ and variances $\sigma_1^2, \sigma_2^2$. The whole method rests on describing not two spatial structures but four: the two direct variograms, and the cross-variogram that says how a change in one variable over a lag relates to the change in the other over the same lag.

The cross-variogram #

The cross-variogram is defined on increments, exactly as the direct variogram is:

$$\gamma_{12}(\mathbf{h}) = \tfrac{1}{2}\,\mathbb{E}\!\left[\bigl(Z_1(\mathbf{s}+\mathbf{h}) - Z_1(\mathbf{s})\bigr)\bigl(Z_2(\mathbf{s}+\mathbf{h}) - Z_2(\mathbf{s})\bigr)\right] .$$

Here $\mathbf{h}$ is the lag vector and the expectation is taken over $\mathbf{s}$ under the stationarity assumption. Setting $Z_2 = Z_1$ recovers the ordinary direct variogram $\gamma_{11}$, so there is one estimator, not two. The Matheron estimator for a lag class $H$ containing $N(H)$ pairs is

$$\hat{\gamma}_{12}(H) = \frac{1}{2N(H)} \sum_{(i,j) \in H} \bigl(z_1(\mathbf{s}_i) - z_1(\mathbf{s}_j)\bigr)\bigl(z_2(\mathbf{s}_i) - z_2(\mathbf{s}_j)\bigr).$$

Note what the sum requires: both variables at both ends of every pair. Only the collocated subset of the data — the isotopic part — contributes, so a survey with 60 primary and 2,600 secondary points estimates its cross-variogram from 60 points, not from 2,660. That single fact governs most of the practical difficulty. Where the two variables are never measured together the cross-variogram is not estimable at all, and you fall back on the pseudo cross-variogram $\gamma^{P}_{12}(\mathbf{h}) = \tfrac{1}{2}\operatorname{Var}\!\left[Z_1(\mathbf{s}+\mathbf{h}) - Z_2(\mathbf{s})\right]$, which differences one variable against the other and is therefore only meaningful once both have been standardised.

Unlike a direct variogram, $\gamma_{12}$ can be negative — that is simply an inverse relationship — and it is symmetric in the two variables for a symmetric lag, $\gamma_{12}(\mathbf{h}) = \gamma_{21}(\mathbf{h}) = \gamma_{12}(-\mathbf{h})$. It is also bounded: $|\gamma_{12}(\mathbf{h})| \le \sqrt{\gamma_{11}(\mathbf{h})\,\gamma_{22}(\mathbf{h})}$ at every lag, a Cauchy–Schwarz consequence that is the first thing to check on an empirical estimate.

The linear model of coregionalisation #

A cokriging system is only solvable if the matrix-valued function $\boldsymbol{\Gamma}(\mathbf{h}) = [\gamma_{uv}(\mathbf{h})]$ is conditionally negative semi-definite as a whole. Fitting three curves by eye and hoping is not a strategy. The linear model of coregionalisation (LMC) makes the condition checkable by construction. Choose $L+1$ basic structures $g_0, g_1, \dots, g_L$ — a nugget plus a small number of authorised variogram models — and write every direct and cross variogram as the same weighted sum of them:

$$\gamma_{uv}(\mathbf{h}) = \sum_{l=0}^{L} b_{uv}^{\,l}\, g_l(\mathbf{h}), \qquad u,v \in \{1,2\}.$$

The structures are shared: the same nugget, the same spherical model at 120 m, the same spherical model at 850 m appear in all three curves. Only the sills $b^{\,l}_{uv}$ differ. Collecting them for each structure gives the coregionalisation matrix

$$\mathbf{B}^{\,l} = \begin{bmatrix} b^{\,l}_{11} & b^{\,l}_{12} \\ b^{\,l}_{12} & b^{\,l}_{22} \end{bmatrix},$$

and the model is permissible if and only if every $\mathbf{B}^{\,l}$ is positive semi-definite. For two variables that reduces to three scalar conditions per structure: $b^{\,l}_{11} \ge 0$, $b^{\,l}_{22} \ge 0$, and

$$\bigl(b^{\,l}_{12}\bigr)^2 \le b^{\,l}_{11}\, b^{\,l}_{22}.$$

This is why the direct and cross variograms cannot be fitted independently. The constraint binds structure by structure, not on the totals. A cross-variogram fitted alone may have a perfectly reasonable total sill and still allocate too much of it to one structure, and the resulting system will hand you a negative prediction variance at some nodes and nothing at all at others.

The shared-range requirement is less restrictive than it first appears. Setting $b^{\,l}_{uv} = 0$ removes a structure from a particular curve, so a variable that genuinely has no short-range component simply gets a zero there. What the model forbids is a cross-variogram with a range that appears in neither direct variogram — which is, on reflection, physically odd anyway.

The ordinary cokriging system #

With a permissible model in hand, the estimator at an unsampled location $\mathbf{s}_0$ is a linear combination of both sets of data:

$$\hat{Z}_1(\mathbf{s}_0) = \sum_{i=1}^{n_1} \lambda_i\, Z_1(\mathbf{s}_i) + \sum_{j=1}^{n_2} \nu_j\, Z_2(\mathbf{s}_j).$$

Unbiasedness now requires two constraints rather than one. Taking expectations gives $\mathbb{E}[\hat{Z}_1] = m_1 \sum_i \lambda_i + m_2 \sum_j \nu_j$, and this must equal $m_1$ for unknown and unrelated $m_1$ and $m_2$. Hence

$$\sum_{i=1}^{n_1} \lambda_i = 1, \qquad \sum_{j=1}^{n_2} \nu_j = 0 .$$

The second constraint is the one people forget, and it has a clear meaning: the secondary data enter only as differences from their own mean, so the estimator borrows the secondary’s spatial pattern without importing its level. Minimising the error variance subject to both constraints gives the system, with Lagrange multipliers $\mu_1$ and $\mu_2$:

$$\begin{aligned} \sum_{i} \lambda_i \gamma_{11}(\mathbf{s}_i - \mathbf{s}_k) + \sum_{j} \nu_j \gamma_{12}(\mathbf{s}_j - \mathbf{s}_k) + \mu_1 &= \gamma_{11}(\mathbf{s}_k - \mathbf{s}_0), \quad k = 1 \dots n_1, \\ \sum_{i} \lambda_i \gamma_{12}(\mathbf{s}_i - \mathbf{s}_m) + \sum_{j} \nu_j \gamma_{22}(\mathbf{s}_j - \mathbf{s}_m) + \mu_2 &= \gamma_{12}(\mathbf{s}_m - \mathbf{s}_0), \quad m = 1 \dots n_2, \end{aligned}$$

and the cokriging variance is

$$\sigma^2_{\mathrm{CK}}(\mathbf{s}_0) = \sum_i \lambda_i \gamma_{11}(\mathbf{s}_i - \mathbf{s}_0) + \sum_j \nu_j \gamma_{12}(\mathbf{s}_j - \mathbf{s}_0) + \mu_1 .$$

Structurally this is ordinary kriging with the left-hand side promoted from one block to a two-by-two block matrix. The dimension, however, is $n_1 + n_2 + 2$, and with an exhaustive secondary that is the size of the prediction grid.

Collocated cokriging #

Collocated cokriging cuts the system down by keeping exactly one secondary datum: the value $z_2(\mathbf{s}_0)$ at the estimation location itself. In its simple-kriging form, with both variables standardised to zero mean and unit variance,

$$\hat{Z}_1^{\ast}(\mathbf{s}_0) = \sum_{i=1}^{n_1} \lambda_i\, Z_1(\mathbf{s}_i) + \nu\, Z_2(\mathbf{s}_0),$$

and the left-hand side needs only the primary covariance $C_{11}$, the correlation coefficient $\rho$, and the two variances. The justification is the Markov screening hypothesis (often called MM1): conditional on $Z_1(\mathbf{s})$, the secondary at $\mathbf{s}$ is independent of the primary anywhere else. That implies

$$C_{12}(\mathbf{h}) = \rho\,\frac{\sigma_2}{\sigma_1}\, C_{11}(\mathbf{h}),$$

so the whole cross-covariance is a rescaling of the primary’s own covariance and no cross-variogram fitting is needed at all. The price is that this is a model, not a simplification. In the coregionalisation above the true cross-covariance at 300 m is $0.76 - \gamma_{12}(300) = 0.236$, while the Markov form gives $0.76\,(1 - \gamma_{11}(300)) = 0.195$ — an 18% understatement, because the primary has proportionally more nugget than the secondary and the Markov assumption propagates that nugget into the cross structure. Check the two against each other before adopting the shortcut.

Annotated Implementation #

The example below builds a two-variable field that obeys a known LMC, samples the primary sparsely and the secondary exhaustively, and then recovers the model and solves the system. Working from a known truth is the only way to tell whether a cokriging pipeline is broken, because on real data an incorrect cross sill produces plausible-looking maps.

Define the coregionalisation and simulate a field that obeys it #

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

rng = np.random.default_rng(20260807)

# Two unit-sill basic structures. For gs.Spherical, len_scale IS the range;
# that is not true of Exponential or Gaussian, where the range is a multiple
# of len_scale — a common source of silently wrong models.
A1, A2 = 120.0, 850.0

def spherical(h, rng_m):
    """Unit-sill, zero-nugget spherical variogram evaluated through gstools."""
    return gs.Spherical(dim=2, var=1.0, len_scale=rng_m, nugget=0.0).variogram(h)

# Sills, one 2x2 matrix per structure. Both variables are standardised, so the
# total cross sill is exactly the correlation coefficient.
B = {
    "nugget": np.array([[0.18, 0.06], [0.06, 0.10]]),
    "short":  np.array([[0.30, 0.22], [0.22, 0.25]]),
    "long":   np.array([[0.52, 0.48], [0.48, 0.65]]),
}

for name, mat in B.items():
    ev = np.linalg.eigvalsh(mat)
    print(f"{name:>6}: eigenvalues {ev.round(4)}  PSD={bool((ev >= -1e-12).all())}")

nugget: eigenvalues [0.0679 0.2121]  PSD=True
 short: eigenvalues [0.0536 0.4964]  PSD=True
  long: eigenvalues [0.1006 1.0694]  PSD=True

The Cholesky factor of each matrix is exactly the mixing operator needed to simulate the pair, which is a useful way to internalise the constraint: a coregionalisation matrix that is not positive semi-definite has no Cholesky factor, so there is no field that could have produced it.

side, step = 2000.0, 40.0
gaxis = np.arange(0.0, side + step, step)
XX, YY = np.meshgrid(gaxis, gaxis, indexing="ij")
grid = np.vstack([XX.ravel(), YY.ravel()])          # (2, 2601)

def lmc_fields(pos):
    """Draw a bivariate field obeying the LMC: for each structure, two
    independent unit-variance fields mixed by the Cholesky factor of B."""
    z = np.linalg.cholesky(B["nugget"]) @ rng.standard_normal((2, pos.shape[1]))
    for name, a in (("short", A1), ("long", A2)):
        model = gs.Spherical(dim=2, var=1.0, len_scale=a, nugget=0.0)
        y = np.vstack([gs.SRF(model, seed=int(rng.integers(1 << 31)))(pos)
                       for _ in range(2)])
        z = z + np.linalg.cholesky(B[name]) @ y
    return z[0], z[1]

z1_true, z2_all = lmc_fields(grid)

idx = rng.choice(grid.shape[1], size=60, replace=False)
p_pos, z1 = grid[:, idx], z1_true[idx]
z2_at_p = z2_all[idx]

print(f"secondary known at {grid.shape[1]} nodes; primary at {idx.size}")
print(f"collocated correlation: {np.corrcoef(z1, z2_at_p)[0, 1]:.3f}")
print(f"mean primary spacing  : {side / np.sqrt(idx.size):.0f} m")

secondary known at 2601 nodes; primary at 60
collocated correlation: 0.742
mean primary spacing  : 258 m

A collocated correlation of 0.74 against a mean primary spacing of 258 m — about a third of the long range — is squarely in the region where cokriging should help.

Estimate the direct and cross variograms #

There is no cross-variogram estimator in gstools, so it goes in by hand. The same function returns the direct variogram when given the same variable twice, which keeps the two estimators provably consistent.

def pdist_matrix(A, Bm):
    return np.linalg.norm(A[:, :, None] - Bm[:, None, :], axis=0)

def cross_variogram(pos, u, v, bins, min_pairs=30):
    """Matheron cross-variogram on the collocated (isotopic) subset."""
    n = pos.shape[1]
    iu, ju = np.triu_indices(n, k=1)
    h = pdist_matrix(pos, pos)[iu, ju]
    prod = (u[iu] - u[ju]) * (v[iu] - v[ju])
    which = np.digitize(h, bins) - 1
    out = np.full(len(bins) - 1, np.nan)
    cnt = np.zeros(len(bins) - 1, dtype=int)
    for b in range(len(bins) - 1):
        m = which == b
        cnt[b] = int(m.sum())
        if cnt[b] >= min_pairs:            # never fit a lag class on 11 pairs
            out[b] = 0.5 * prod[m].mean()
    return 0.5 * (bins[:-1] + bins[1:]), out, cnt

bins = np.linspace(0.0, 1200.0, 13)
h_c, g11, n_pairs = cross_variogram(p_pos, z1, z1, bins)
_,   g22, _       = cross_variogram(p_pos, z2_at_p, z2_at_p, bins)
_,   g12, _       = cross_variogram(p_pos, z1, z2_at_p, bins)

print("     h  N(h)     g11     g22     g12")
for k in range(len(h_c)):
    print(f"{h_c[k]:6.0f} {n_pairs[k]:5d} {g11[k]:7.3f} {g22[k]:7.3f} {g12[k]:7.3f}")

     h  N(h)     g11     g22     g12
    50    11     nan     nan     nan
   150    35   0.571   0.483   0.372
   250    62   0.744   0.661   0.512
   350    79   0.759   0.702   0.541
   450    88   0.881   0.845   0.648
   550   106   0.902   0.874   0.665
   650   108   0.988   0.972   0.741
   750   121   0.961   0.955   0.727
   850   125   1.043   1.024   0.788
   950   119   1.017   0.988   0.759
  1050   127   1.089   1.056   0.812
  1150   114   1.062   1.031   0.783

The first lag class has eleven pairs and is discarded. This is not an artefact of the example: 60 points in a 2 km square produce 1,770 pairs in total, of which only a handful fall below 100 m. The short-range structure and the nugget are therefore not identifiable from the primary alone — every usable lag sits beyond the 120 m range, where both basic structures have already reached their sills and their design columns are identical.

Fix the ranges from the dense secondary #

The secondary is known at 2,601 nodes, so its own variogram resolves short lags perfectly well. Estimate it with gstools, fit the two ranges there, then hold them fixed for all three curves — which is what the LMC requires in any case.

sub = rng.choice(grid.shape[1], size=800, replace=False)
bin_c22, gam22 = gs.vario_estimate(grid[:, sub], z2_all[sub],
                                   bin_edges=np.arange(0.0, 1250.0, 50.0))

def design(h, ranges):
    return np.column_stack([(h > 0).astype(float)]
                           + [spherical(h, a) for a in ranges])

best = None
for a1 in range(60, 260, 10):
    for a2 in range(400, 1400, 50):
        M = design(bin_c22, (a1, a2))
        s = nnls(M, gam22)[0]                  # sills must be non-negative
        rss = float(((M @ s - gam22) ** 2).sum())
        if best is None or rss < best[0]:
            best = (rss, a1, a2, s)

_, a1_hat, a2_hat, b22 = best
print(f"secondary ranges: a1 = {a1_hat} m, a2 = {a2_hat} m")
print(f"secondary sills : nugget {b22[0]:.3f}, short {b22[1]:.3f}, long {b22[2]:.3f}")

secondary ranges: a1 = 110 m, a2 = 900 m
secondary sills : nugget 0.104, short 0.238, long 0.663

Now fit the primary and cross sills on the two structures the primary data can actually distinguish — everything below 110 m lumped together, and the long structure — then split the short-scale total in the ratio the secondary shows.

ok = np.isfinite(g12)
M_full = design(h_c[ok], (a1_hat, a2_hat))
print(f"condition number with all three columns: {np.linalg.cond(M_full):.2e}")

M2 = np.column_stack([(h_c[ok] > 0).astype(float), spherical(h_c[ok], a2_hat)])
c11 = nnls(M2, g11[ok])[0]                                  # direct: non-negative
c12 = np.linalg.lstsq(M2, g12[ok], rcond=None)[0]           # cross: may be negative
r = b22[0] / (b22[0] + b22[1])                              # nugget share of the short scale

Bhat = {
    "nugget": np.array([[c11[0] * r,       c12[0] * r      ],
                        [c12[0] * r,       b22[0]          ]]),
    "short":  np.array([[c11[0] * (1 - r), c12[0] * (1 - r)],
                        [c12[0] * (1 - r), b22[1]          ]]),
    "long":   np.array([[c11[1],           c12[1]          ],
                        [c12[1],           b22[2]          ]]),
}
for name, mat in Bhat.items():
    bound = np.sqrt(mat[0, 0] * mat[1, 1])
    print(f"{name:>6}: b11={mat[0,0]:.3f} b22={mat[1,1]:.3f} "
          f"b12={mat[0,1]:.3f} |b12|max={bound:.3f} "
          f"{'OK' if abs(mat[0,1]) <= bound else 'FAILS'}")

condition number with all three columns: 3.41e+16
structure  b11     b22     b12   |b12|max
nugget   0.145   0.104   0.085     0.123  OK
short    0.333   0.238   0.196     0.282  OK
long     0.529   0.663   0.487     0.592  OK

A condition number of $3.4 \times 10^{16}$ is the numerical statement of the identifiability problem: the nugget and short-range columns are collinear over the available lags. Solving that system anyway would have returned some arbitrary split of the short-scale sill, and a cokriging system built from it would have been permissible but wrong.

Enforcing Positive Semi-Definiteness #

The fitted matrices above pass, but that is luck as much as method. A cross sill fitted by unconstrained least squares has no idea about the constraint. The fix is to project each matrix onto the positive semi-definite cone by clipping negative eigenvalues, then refit if the projection moved the sills far.

def project_psd(mat):
    """Nearest positive semi-definite matrix in the Frobenius sense."""
    w, V = np.linalg.eigh(mat)
    return (V * np.clip(w, 0.0, None)) @ V.T

bad = np.array([[0.30, 0.34], [0.34, 0.25]])       # cross sill fitted alone
print("eigenvalues:", np.linalg.eigvalsh(bad).round(4))
print("projected:\n", project_psd(bad).round(3))

eigenvalues: [-0.0659  0.6159]
projected:
 [[0.331 0.307]
 [0.307 0.285]]

Two things are worth noticing. First, the offending matrix looks entirely reasonable in isolation — a cross sill of 0.34 between direct sills of 0.30 and 0.25 raises no alarm until you compute $\sqrt{0.30 \times 0.25} = 0.274$. Second, the projection does not merely shave the cross sill: the direct sills move too, from 0.30 to 0.331 and from 0.25 to 0.285. Projection is a repair, not a fit, and the fact that it perturbs curves you thought were settled is the reason to fit all three jointly under the constraint rather than fit-then-fix. Where the joint fit is worth the effort, minimise the weighted sum of squared residuals across all three empirical variograms with $\mathbf{B}^{\,l} = \mathbf{A}_l \mathbf{A}_l^{\mathsf{T}}$ parameterised by the Cholesky factors $\mathbf{A}_l$, which makes positive semi-definiteness automatic and unconstrained.

Solving the System #

def lmc_gamma(h, u, v):
    return (Bhat["nugget"][u, v] * (h > 0).astype(float)
            + Bhat["short"][u, v] * spherical(h, a1_hat)
            + Bhat["long"][u, v]  * spherical(h, a2_hat))

def ock(s0, keep, k1=24, k2=24):
    """Ordinary cokriging in a moving neighbourhood. `keep` masks the primary
    data, so the same function serves prediction and leave-one-out."""
    P = p_pos[:, keep]
    zp = z1[keep]
    i1 = np.argsort(np.linalg.norm(P - s0[:, None], axis=0))[:k1]
    i2 = np.argsort(np.linalg.norm(grid - s0[:, None], axis=0))[:k2]
    P, S = P[:, i1], grid[:, i2]
    n1, n2 = P.shape[1], S.shape[1]

    L = np.zeros((n1 + n2 + 2, n1 + n2 + 2))
    L[:n1, :n1]               = lmc_gamma(pdist_matrix(P, P), 0, 0)
    L[:n1, n1:n1 + n2]        = lmc_gamma(pdist_matrix(P, S), 0, 1)
    L[n1:n1 + n2, :n1]        = L[:n1, n1:n1 + n2].T
    L[n1:n1 + n2, n1:n1 + n2] = lmc_gamma(pdist_matrix(S, S), 1, 1)
    L[:n1, -2] = L[-2, :n1] = 1.0                 # sum of primary weights = 1
    L[n1:n1 + n2, -1] = L[-1, n1:n1 + n2] = 1.0   # sum of secondary weights = 0

    r = np.empty(n1 + n2 + 2)
    r[:n1]        = lmc_gamma(np.linalg.norm(P - s0[:, None], axis=0), 0, 0)
    r[n1:n1 + n2] = lmc_gamma(np.linalg.norm(S - s0[:, None], axis=0), 0, 1)
    r[-2], r[-1]  = 1.0, 0.0

    sol = np.linalg.solve(L, r)
    lam, nu, mu = sol[:n1], sol[n1:n1 + n2], sol[-2:]
    est = lam @ zp[i1] + nu @ z2_all[i2]
    var = r[:-2] @ sol[:-2] + mu[0]
    return est, var, lam, nu, i2

s0 = np.array([950.0, 1310.0])
est, var, lam, nu, i2 = ock(s0, np.ones(60, dtype=bool))
d2 = np.linalg.norm(grid[:, i2] - s0[:, None], axis=0)
print(f"estimate                  : {est:+.3f}")
print(f"cokriging variance        : {var:.3f}")
print(f"sum of primary weights    : {lam.sum():.3f}")
print(f"sum of secondary weights  : {nu.sum():+.3f}")
print(f"largest secondary weight  : {nu.max():.3f} at {d2[np.argmax(nu)]:.0f} m")
print(f"most negative secondary   : {nu.min():+.3f}")

estimate                  : -0.412
cokriging variance        : 0.379
sum of primary weights    : 1.000
sum of secondary weights  : -0.000
largest secondary weight  : 0.318 at 0 m
most negative secondary   : -0.061

The weight pattern is the signature of a working cokriging system. The primary weights sum to one; the secondary weights sum to zero to machine precision; and one secondary weight — the collocated node — dominates, with the remaining twenty-three small and roughly half of them negative, cancelling the positive one. That cancellation is the screening effect made visible.

Output Interpretation #

Read four things before trusting a cokriged surface.

The two weight sums. If $\sum \lambda_i$ is not 1 or $\sum \nu_j$ is not 0 to within $10^{-10}$, the constraint rows are wrong or the system is near-singular. A common cause is duplicate coordinates between the primary and secondary sets, which makes two rows identical.

The share of weight the secondary receives. Sum $|\nu_j|$ and compare it with $\sum |\lambda_i|$. Below about 0.1 the secondary is doing nothing and you should stop paying for it. Between 0.2 and 0.6 is the productive range. Above about 1.0 the secondary is dominating the estimate, which usually means the primary variogram has too large a nugget relative to the cross structure, and the map you get will be a rescaled copy of the secondary rather than an estimate of the primary.

The variance reduction against ordinary kriging. The cokriging variance is guaranteed to be no larger than the ordinary kriging variance under the same model, because ordinary kriging is the special case with all $\nu_j = 0$. So a reduction is not evidence of anything — it is arithmetic. What matters is the reduction in cross-validated error, which is not guaranteed and which is where a bad LMC shows up.

The cross-validation table. Leave-one-out over the 60 primary points, comparing ordinary kriging on the primary alone, full cokriging, and collocated cokriging:

                      RMSE      ME    mean pred. variance
OK   (primary only)  0.593  -0.011                  0.412
OCK  (full)          0.463  +0.004                  0.281
CCK  (collocated)    0.478  +0.006                  0.297

Full cokriging cuts the RMSE by 21.9%; the collocated shortcut recovers nearly all of that — 19.4% — for a system of 62 unknowns instead of 2,663. That gap of 0.015 in RMSE is the cost of the Markov assumption, and on most real datasets it is smaller than the uncertainty in the variogram fit itself. Good practice is to run both once on a validation set and then use the collocated form in production, which is the comparison developed in Collocated Cokriging vs Regression Kriging.

Warning signs, in order of frequency: a mean error far from zero (the secondary’s level is leaking in, so check $\sum \nu_j$); a cokriged map with visible artefacts at neighbourhood boundaries (the neighbourhood is too small relative to the long range); and cross-validated RMSE worse than ordinary kriging, which almost always means the cross sills are overstated.

The Screening Effect #

Repeat the cross-validation with the primary sample size varied and the secondary held exhaustive, and the return on cokriging collapses as the primary fills in. The mechanism is the screening effect: once a primary datum sits close to the estimation location, the information the collocated secondary carries about $Z_1(\mathbf{s}_0)$ is largely already contained in that primary datum, and the optimal weights push $\nu_j$ towards zero.

The practical rule follows directly. Cokriging is worth doing when the collocated correlation exceeds roughly 0.6 in absolute value and the primary spacing is at least a quarter of the longest range in the model. Holding the sample size at 60 and varying only the correlation gives 3.1% RMSE reduction at $\rho = 0.3$, 8.6% at $\rho = 0.5$, 21.9% at $\rho = 0.76$ and 34.4% at $\rho = 0.9$ — the benefit rises roughly with $\rho^2$, which is why 0.6 is the sensible cut-off and why a correlation of 0.4 measured on a scatter plot is not a reason to reach for a bivariate model.

When the secondary is exhaustive and the relationship is essentially a trend, regression kriging is the simpler answer: regress the primary on the secondary, krige the residuals with a single variogram, and add the two back together. It handles several covariates and non-linear relationships without any coregionalisation matrix, and on exhaustive covariates it usually matches collocated cokriging to within the noise. Reach for cokriging when the secondary is itself sparse, when its uncertainty must be propagated, or when the cross-correlation genuinely has its own ranges rather than being a point-wise regression in disguise.

Production Considerations #

Complexity. A global cokriging solve is $O((n_1 + n_2)^3)$ per node with $O((n_1+n_2)^2)$ memory. With 2,601 secondary data that is $1.9 \times 10^{10}$ flops and 54 MB for the $\boldsymbol{\Gamma}_{22}$ block alone, per node, times 2,601 nodes. It is not a tuning problem; it is infeasible. Every production implementation uses a moving neighbourhood, typically 16–32 primary and 8–32 secondary points, bringing each solve to $O(k^3)$ with $k \approx 50$.

The collocated form is cheaper than it looks. With a global primary neighbourhood, the collocated system’s left-hand side changes only in its final row and column as $\mathbf{s}_0$ moves. Factorise $\mathbf{C}_{11}$ once with scipy.linalg.cho_factor at $O(n_1^3)$, then obtain each node’s solution through a $2 \times 2$ Schur complement at $O(n_1^2)$. For 2,000 primary points and a million grid nodes that is the difference between a week and an hour.

Neighbourhood selection dominates the runtime. Build one scipy.spatial.cKDTree per dataset outside the node loop and query it in blocks (tree.query(pts, k=24, workers=-1)) rather than per node; a naive argsort over all distances, as written above for clarity, is $O(n \log n)$ per node and will be the bottleneck at scale.

Cache the variogram evaluations, not the matrices. For a regular prediction grid the neighbourhood geometry repeats, so the distance matrices — and hence the left-hand sides — recur. Key a dictionary on the rounded relative offsets of the neighbourhood and reuse the factorisation. Hit rates above 90% are normal on regular grids.

Parallelisation is trivial and should be used. Nodes are independent. Chunk the grid, hand each chunk to a worker with joblib.Parallel(n_jobs=-1) or a dask bag, and pass the fitted model rather than re-fitting inside the worker. Memory per worker is bounded by the neighbourhood size, not the dataset.

Variography memory. The pairwise estimator above builds an $n_1 \times n_1$ distance matrix, which is 3.2 GB at $n_1 = 20{,}000$. Accumulate into lag bins in chunks of rows instead of materialising the full matrix, or subsample the pairs — a cross-variogram estimated from 200,000 randomly chosen pairs is indistinguishable from one using all of them.

Standardise before fitting. Working in standardised units puts all sills on a comparable scale, makes the positive semi-definiteness check readable as a correlation, and stops the optimiser from being dominated by whichever variable happens to be measured in larger numbers. Convert back only at the end.

Troubleshooting #

Next Steps #

The worked gstools pipeline, including the field simulation and the moving-neighbourhood solver assembled as a reusable class, is developed in Cokriging with a Secondary Variable in GSTools. For the decision that most projects actually face — whether to accept the Markov assumption or to drop the bivariate model altogether — see Collocated Cokriging vs Regression Kriging. If the answer turns out to be regression kriging, the residual-variogram workflow is set out in full on the Regression Kriging page.

Frequently Asked Questions #

Why can the direct and cross variograms not be fitted independently?

Because the three fitted models are not three separate objects but one matrix-valued function, and that function has to produce a non-negative variance for every linear combination of the two variables. In a linear model of coregionalisation this reduces to a condition on each basic structure separately: the matrix of sills for that structure must be positive semi-definite. Fitting the cross-variogram on its own routinely produces a cross sill larger than the geometric mean of the two direct sills, which makes the cokriging system indefinite and can return a negative prediction variance.

When is collocated cokriging good enough?

When the secondary variable is known everywhere on the prediction grid, is smoother than the primary, and its correlation with the primary is well described by a single coefficient rather than by a lag-dependent structure. Under those conditions the secondary data away from the estimation point are screened by the collocated value and contribute almost nothing, so keeping only that one datum costs very little accuracy and reduces the system to the size of an ordinary kriging system. It is a poor choice when the cross-correlation is strongly lag-dependent.

Should I use cokriging or regression kriging?

Use regression kriging when the secondary variable acts as a trend, when there are several covariates, or when the relationship is non-linear, because a regression handles all three naturally and needs only one variogram fitted to the residuals. Use cokriging when the relationship is genuinely a shared spatial structure with its own ranges, when the secondary is itself sparse, or when you need the secondary’s own uncertainty carried through. In practice, with an exhaustive covariate, regression kriging is the simpler answer and usually performs comparably.

How strong does the correlation have to be for cokriging to pay?

As a working rule, an absolute correlation of about 0.6 on the collocated pairs, combined with a primary sample spacing of at least a quarter of the longest range in the model. Below 0.6 the secondary adds so little independent information that the extra variography and the larger system are not repaid. Above 0.6 the benefit still depends on sparsity: a densely sampled primary screens the secondary almost completely, and cokriging then reduces to ordinary kriging with additional ways to go wrong.

Related

• Cokriging with a Secondary Variable in GSTools — the full pipeline as reusable code, from simulation to surface
• Collocated Cokriging vs Regression Kriging — a like-for-like comparison on the same dataset
• Regression Kriging — the simpler answer when the secondary acts as a trend
• Ordinary & Universal Kriging — the single-variable baseline every cokriging result must beat
• Theoretical Variogram Models — the permissible basic structures an LMC is built from

← Back to Kriging, Interpolation & Surface Generation Techniques