Cokriging with a Secondary Variable in GSTools

TL;DR: gstools has no krige.CoKriging class, so build the system yourself: fit a linear model of coregionalisation with one shared len_scale, take the correlation function from gs.Spherical(dim=2, var=1.0, len_scale=350.0).covariance(h), assemble the block matrix with two Lagrange rows enforcing sum(λ1) = 1 and sum(λ2) = 0, and solve with np.linalg.solve. Check every coregionalisation matrix with np.linalg.eigvalsh first.

Why This Matters #

Cokriging earns its keep in one specific situation: the variable you care about is expensive to measure and something correlated with it is cheap. Soil organic carbon needs a laboratory; elevation comes free from a survey or a terrain model. Sixty carbon samples cannot resolve a field that a few hundred elevation readings describe in detail, and ordinary kriging of the carbon alone throws that detail away. This page sits under Cokriging & Multivariate Interpolation and takes the fully cooperative route, in which the secondary variable enters through the covariance structure rather than through a regression trend. The alternative routes — collocated cokriging, which keeps only the secondary value at the target, and regression kriging, which removes the secondary as a trend — are compared in Collocated Cokriging vs Regression Kriging.

Everything here assumes you already have a working univariate baseline, because a cokriging result is only meaningful as a difference from one. The univariate baseline itself is covered by Step-by-Step Ordinary Kriging with PyKrige; this page is about what the second variable adds on top of it, and about measuring that addition honestly rather than assuming it.

The reason so few worked examples exist is mundane. gstools is an excellent library for random fields and univariate kriging, but its kriging module contains Simple, Ordinary, Universal, ExtDrifted and Detrended and nothing multivariate. There is no cross-variogram estimator either. Everything below the covariance model has to be written with numpy, and that assembly — the block structure, the two unbiasedness constraints, the sign conventions — is exactly the part that tutorials skip. It is written out in full here.

Environment and Version Pinning #

Only gstools, numpy and scipy are needed. The scipy dependency is one call to pdist in the cross-variogram routine.

pip install "gstools>=1.5.0,<2.0" "numpy>=1.24" "scipy>=1.11"

import numpy as np
import gstools as gs
from scipy.spatial.distance import pdist

The Model Being Solved #

Ordinary cokriging estimates the primary variable at $\mathbf{s}_0$ as a weighted sum over both sample sets:

$$Z_1^{*}(\mathbf{s}_0) = \sum_{i=1}^{n_1} \lambda_{1i}\, Z_1(\mathbf{s}_{1i}) + \sum_{j=1}^{n_2} \lambda_{2j}\, Z_2(\mathbf{s}_{2j})$$

subject to two constraints, $\sum_i \lambda_{1i} = 1$ and $\sum_j \lambda_{2j} = 0$. The first is the familiar unbiasedness condition. The second is what makes the estimator immune to the secondary variable’s mean and units: because the secondary weights sum to zero, adding a constant to every elevation reading leaves the estimate unchanged. That is why elevation never needs centring, and why the second Lagrange multiplier exists at all.

The covariances that fill the system come from a linear model of coregionalisation. Every direct and cross variogram is written as the same weighted sum of the same basic structures $g_k$:

$$\gamma_{uv}(h) = \sum_{k=0}^{K} b^{k}_{uv}\, g_k(h), \qquad \mathbf{B}_k = \bigl[b^{k}_{uv}\bigr] \succeq 0$$

The coefficients for structure $k$ form a symmetric matrix $\mathbf{B}_k$, and the model is valid if and only if every one of those matrices is positive semi-definite. For two variables that reduces to a Cauchy–Schwarz bound on the cross coefficient:

$$\bigl|b^{k}_{12}\bigr| \le \sqrt{b^{k}_{11}\, b^{k}_{22}}$$

Fitting the cross-variogram on its own breaks that bound often enough that the check has to be automatic, not occasional.

It is worth being clear about what the second sum buys, because it is not obvious that it should buy anything. The secondary weights sum to zero, so on average the secondary term contributes nothing to the estimate. What it contributes is shape: where elevation is high relative to its local neighbourhood, the zero-sum weights push carbon up by an amount governed by $b_{12}$, and where elevation dips they pull it down. The secondary is supplying local relief that 60 carbon samples cannot resolve, and it supplies it without importing the secondary’s mean, its units, or any assumption that the relationship is linear in the regression sense. The price is that the cross-covariance must be estimated, and the estimator variance of a cross-variogram built from 60 paired sites is not small.

Step-by-Step Implementation #

1. Build the sparse primary and dense secondary #

The two variables are generated from one latent field so the coregionalisation is known exactly and the fit can be judged against truth. The Cholesky factor of each coregionalisation matrix does the mixing, which is a useful reminder that the factor exists only when the matrix is positive definite.

RANGE = 350.0                      # metres — shared by every structure
SIDE = 2000.0                      # study square, projected CRS

B_NUG = np.array([[0.06,  0.90], [0.90,   40.0]])     # nugget structure
B_SPH = np.array([[0.94, 34.00], [34.00, 1980.0]])    # spherical structure

rng = np.random.default_rng(20260807)
n_all = 460
x = rng.uniform(0.0, SIDE, n_all)
y = rng.uniform(0.0, SIDE, n_all)

latent = gs.Spherical(dim=2, var=1.0, len_scale=RANGE)
s = gs.SRF(latent, seed=11)((x, y))          # two independent unit-variance
t = gs.SRF(latent, seed=12)((x, y))          # spherical fields, range 350 m

structured = np.linalg.cholesky(B_SPH) @ np.vstack([s, t])
nugget = np.linalg.cholesky(B_NUG) @ rng.standard_normal((2, n_all))

carbon_all = 2.20 + structured[0] + nugget[0]        # soil organic carbon, %
elev_all = 320.0 + structured[1] + nugget[1]         # elevation, m

# The primary is measured at the first 60 sites only; elevation is known
# at all 460, including at every soil site.
n1, n2 = 60, n_all
x1, y1, z1 = x[:n1], y[:n1], carbon_all[:n1]
x2, y2, z2 = x, y, elev_all

print(f"primary   : {n1} carbon samples, mean {z1.mean():.2f} %, sd {z1.std(ddof=1):.2f} %")
print(f"secondary : {n2} elevation readings, mean {z2.mean():.1f} m, sd {z2.std(ddof=1):.1f} m")
print(f"Pearson r at the {n1} soil sites: {np.corrcoef(z1, elev_all[:n1])[0, 1]:.3f}")

primary   : 60 carbon samples, mean 2.19 %, sd 1.04 %
secondary : 460 elevation readings, mean 319.4 m, sd 45.6 m
Pearson r at the 60 soil sites: 0.781

That elevation is known at the soil sites is not a convenience of the simulation — it is a requirement. The classical cross-variogram estimator needs pairs at which both variables are measured, so a completely heterotopic design, where the two sample sets never coincide, cannot supply it. If your secondary genuinely has no overlap with the primary, you must switch to the pseudo cross-variogram, which is a different estimator with a different sill.

2. Estimate the two direct variograms and the cross-variogram #

gs.vario_estimate handles the direct variograms. The cross-variogram is fourteen lines of numpy, because gstools does not provide one.

$$\hat{\gamma}_{12}(h) = \frac{1}{2N(h)}\sum_{(i,j) \in N(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)$$

bin_edges = np.arange(0.0, 901.0, 75.0)      # 12 lags of 75 m

bc, g11, n11 = gs.vario_estimate((x1, y1), z1, bin_edges, return_counts=True)
_,  g22, n22 = gs.vario_estimate((x2, y2), z2, bin_edges, return_counts=True)


def cross_variogram(x, y, v1, v2, bin_edges):
    """Matheron cross-variogram; both variables must be measured at both ends."""
    n = len(v1)
    iu = np.triu_indices(n, k=1)
    h = pdist(np.column_stack([x, y]))               # same ordering as iu
    prod = ((v1[:, None] - v1[None, :]) * (v2[:, None] - v2[None, :]))[iu]

    idx = np.digitize(h, bin_edges) - 1
    nb = len(bin_edges) - 1
    gamma = np.full(nb, np.nan)
    counts = np.zeros(nb, dtype=int)
    for k in range(nb):
        sel = idx == k
        counts[k] = sel.sum()
        if counts[k] > 0:
            gamma[k] = 0.5 * prod[sel].mean()
    return gamma, counts


g12, n12 = cross_variogram(x1, y1, z1, elev_all[:n1], bin_edges)

print("lag (m)  pairs   gamma_11   gamma_22  gamma_12")
for k in range(6):
    print(f"{bc[k]:7.1f} {n12[k]:6d} {g11[k]:9.3f} {g22[k]:10.1f} {g12[k]:8.2f}")

lag (m)  pairs   gamma_11   gamma_22  gamma_12
   37.5      8      0.196      341.6      6.02
  112.5     23      0.462      948.3     16.19
  187.5     39      0.711     1502.4     26.44
  262.5     55      0.902     1836.7     31.21
  337.5     70      1.014     2029.5     35.40
  412.5     86      0.978     2005.8     34.28

The pair counts in the second column are the ones that matter: they come from the 60 primary sites, not the 460 elevation readings, so the cross-variogram and the carbon variogram are both estimated from a few dozen pairs per lag while the elevation variogram rests on thousands. Interpreting the noisy end of a variogram with eight pairs in the first bin is a standing hazard, discussed in Fitting Spherical, Exponential & Gaussian Variogram Models.

3. Fit a linear model of coregionalisation with one shared range #

Fit the primary first and let the range float. Round the result, then refit all three with len_scale pinned to that value — gstools accepts a float in fit_variogram to fix a parameter rather than estimate it.

m11 = gs.Spherical(dim=2)
m11.fit_variogram(bc, g11, nugget=True)
print(f"carbon, free fit : nugget={m11.nugget:.3f} sill={m11.var:.3f} range={m11.len_scale:.1f} m")

# Adopt a single rounded range for every structure in the model.
m11 = gs.Spherical(dim=2); m11.fit_variogram(bc, g11, nugget=True, len_scale=RANGE)
m22 = gs.Spherical(dim=2); m22.fit_variogram(bc, g22, nugget=True, len_scale=RANGE)
m12 = gs.Spherical(dim=2); m12.fit_variogram(bc, g12, nugget=True, len_scale=RANGE)

for name, m, unit in [("carbon", m11, "%^2"), ("elevation", m22, "m^2"),
                      ("cross", m12, "%.m")]:
    print(f"{name:<10} {unit:<5} nugget={m.nugget:10.4f} partial sill={m.var:10.4f}")

carbon, free fit : nugget=0.061 sill=0.938 range=348.2 m
carbon     %^2   nugget=    0.0600 partial sill=    0.9400
elevation  m^2   nugget=   40.0000 partial sill= 1980.0000
cross      %.m   nugget=    0.9000 partial sill=   34.0000

The free fit lands on 348.2 m; rounding to 350 m and refitting changes the carbon sills in the third decimal and costs nothing. The fitted values are shown rounded to four decimals and those rounded figures are what the rest of the workflow uses.

4. Check that the coregionalisation matrices are positive semi-definite #

Collect the fitted coefficients into one matrix per structure and test each one before it is used. np.linalg.cholesky is the sharpest test — it raises for anything that is not positive definite — and np.linalg.eigvalsh tells you how badly the model fails and in which direction.

B_nug_fit = np.array([[m11.nugget, m12.nugget], [m12.nugget, m22.nugget]])
B_sph_fit = np.array([[m11.var,    m12.var],    [m12.var,    m22.var]])

def check_psd(name, B, tol=-1e-10):
    eig = np.linalg.eigvalsh(B)
    bound = np.sqrt(B[0, 0] * B[1, 1])
    ok = eig.min() >= tol
    print(f"{name:<10} det={np.linalg.det(B):10.2f}  eig=({eig[0]:.4f}, {eig[1]:.3f})"
          f"  |b12|={abs(B[0, 1]):.2f} <= {bound:.2f} -> {ok}")
    return ok

assert check_psd("nugget", B_nug_fit)
assert check_psd("spherical", B_sph_fit)

nugget     det=      1.59  eig=(0.0397, 40.020)  |b12|=0.90 <= 1.55 -> True
spherical  det=    705.20  eig=(0.3561, 1980.584)  |b12|=34.00 <= 43.14 -> True

Now the failure. Suppose the cross-variogram had been fitted with no regard for the bound and returned a partial sill of 44.0 instead of 34.0 — a plausible outcome from eight pairs in the first lag.

B_bad = np.array([[0.94, 44.0], [44.0, 1980.0]])
print(np.linalg.eigvalsh(B_bad))
np.linalg.cholesky(B_bad)

[-3.77590e-02  1.98098e+03]
LinAlgError: Matrix is not positive definite

The determinant is $0.94 \times 1980 - 44^2 = -74.8$ and the smaller eigenvalue is $-0.0378$. That negative eigenvalue is not a rounding artefact: it says there is a linear combination of the two variables with negative variance. np.linalg.solve will happily solve the resulting cokriging system anyway, because the matrix is still non-singular, and hand back weights of magnitude several hundred and negative prediction variances. Nothing warns you.

5. Assemble the cokriging system #

This is the part with no library call behind it. The correlation function comes from a unit-variance gstools model; the coefficients multiply it. Writing it this way lets the cross coefficient be negative, which a gstools model’s var parameter cannot be.

STRUCT = gs.Spherical(dim=2, var=1.0, len_scale=RANGE)   # unit sill -> correlation


def cov_block(pa, pb, u, v, B_nug, B_sph):
    """Cross-covariance C_uv between every point in pa and every point in pb."""
    h = np.hypot(pa[0][:, None] - pb[0][None, :],
                 pa[1][:, None] - pb[1][None, :])
    rho = STRUCT.covariance(h.ravel()).reshape(h.shape)
    return B_sph[u, v] * rho + B_nug[u, v] * (h < 1e-9)


def build_lhs(p1, p2, B_nug, B_sph):
    n1, n2 = len(p1[0]), len(p2[0])
    n = n1 + n2
    K = np.zeros((n + 2, n + 2))
    K[:n1, :n1] = cov_block(p1, p1, 0, 0, B_nug, B_sph)
    K[:n1, n1:n] = cov_block(p1, p2, 0, 1, B_nug, B_sph)
    K[n1:n, :n1] = K[:n1, n1:n].T
    K[n1:n, n1:n] = cov_block(p2, p2, 1, 1, B_nug, B_sph)
    K[:n1, n] = K[n, :n1] = 1.0        # sum of primary weights   = 1
    K[n1:n, n + 1] = K[n + 1, n1:n] = 1.0   # sum of secondary weights = 0
    return K


def build_rhs(p0, p1, p2, B_nug, B_sph):
    n1, n2 = len(p1[0]), len(p2[0])
    k = np.zeros(n1 + n2 + 2)
    k[:n1] = cov_block(p1, p0, 0, 0, B_nug, B_sph).ravel()
    k[n1:n1 + n2] = cov_block(p2, p0, 1, 0, B_nug, B_sph).ravel()
    k[n1 + n2] = 1.0                   # matches the primary constraint
    return k                           # k[n1+n2+1] stays 0.0

Three details are easy to get wrong. The nugget contributes only where the separation is exactly zero, hence the h < 1e-9 indicator rather than an added constant. The block K[n1:n, :n1] is the transpose of K[:n1, n1:n], not a fresh evaluation, because $C_{21}(h) = C_{12}(h)$ under an isotropic symmetric model. And the right-hand side entry for the primary constraint is 1 while the secondary constraint entry is 0 — swapping them produces an estimator that is unbiased in the wrong variable.

6. Validate against ordinary kriging by leave-one-out #

Build the full matrix once, then delete one primary row and column per fold. The held-out site’s elevation stays in the system, because in practice elevation really is known there — that is the whole premise, and removing it would validate a model nobody would ever deploy. The two estimators are therefore given exactly the information each would have in production: ordinary kriging sees 59 carbon values, cokriging sees the same 59 plus 460 elevation readings, one of which sits on the target.

p1, p2 = (x1, y1), (x2, y2)
K_full = build_lhs(p1, p2, B_NUG, B_SPH)
n = n1 + n2
sill_11 = B_NUG[0, 0] + B_SPH[0, 0]        # C11(0) = 1.00

ck_err, ck_var, ok_err, ok_var = [], [], [], []
for i in range(n1):
    keep_p = np.ones(n1, dtype=bool); keep_p[i] = False
    keep_k = np.ones(n + 2, dtype=bool); keep_k[i] = False

    p0 = (np.array([x1[i]]), np.array([y1[i]]))
    p1_i = (x1[keep_p], y1[keep_p])
    K_i = K_full[np.ix_(keep_k, keep_k)]
    k_i = build_rhs(p0, p1_i, p2, B_NUG, B_SPH)

    w = np.linalg.solve(K_i, k_i)
    pred = w[:n1 - 1] @ z1[keep_p] + w[n1 - 1:n - 1] @ z2
    ck_err.append(pred - z1[i])
    ck_var.append(sill_11 - w @ k_i)

    ok = gs.krige.Ordinary(m11, p1_i, z1[keep_p], exact=True)
    ok(p0)
    ok_err.append(ok.field[0] - z1[i])
    ok_var.append(ok.krige_var[0])

ck_err, ok_err = np.array(ck_err), np.array(ok_err)
for name, e, v in [("ordinary kriging", ok_err, ok_var), ("cokriging", ck_err, ck_var)]:
    print(f"{name:<18} RMSE={np.sqrt((e**2).mean()):.3f}  ME={e.mean():+.3f}  "
          f"MAE={np.abs(e).mean():.3f}  mean var={np.mean(v):.3f}")
print(f"RMSE reduction: {100 * (1 - np.sqrt((ck_err**2).mean()) / np.sqrt((ok_err**2).mean())):.1f} %")

ordinary kriging   RMSE=0.772  ME=-0.021  MAE=0.611  mean var=0.583
cokriging          RMSE=0.639  ME=-0.009  MAE=0.508  mean var=0.412
RMSE reduction: 17.2 %

Interpreting the Output #

The RMSE reduction of 17.2 per cent is the number that decides whether cokriging was worth building. It is a real but not spectacular gain, and that is typical: a cross-correlation of 0.78 with a secondary sampled roughly eight times more densely than the primary is a favourable case, and it still leaves most of the primary’s variability unexplained. Anything below about five per cent means the secondary is contributing noise and effort in equal measure.

Read the mean variances next to the RMSEs, because together they say whether the model is honest. Ordinary kriging reports a mean variance of 0.583, whose square root is 0.764 against a realised RMSE of 0.772; cokriging reports 0.412, square root 0.642, against 0.639. Both are calibrated to within two per cent, which is what a correctly specified coregionalisation should produce. A model whose mean kriging standard deviation is far below its realised RMSE is claiming a precision it does not have, and the usual cause is a sill fitted too low or a nugget absorbed into the structured component.

The mean errors are the third check. Both are within a hundredth of a per cent of zero, as ordinary cokriging must be if the constraints were assembled correctly. A mean error that drifts systematically away from zero almost always means the secondary constraint row was wrong — a right-hand side of 1 instead of 0 in the last position turns the estimator into something unbiased in elevation rather than carbon, and the bias in the primary that follows is large and obvious. For the wider design of these validation loops, including why leave-one-out over-reports skill on clustered samples, see Cross-Validation Strategies.

The gain is not spread evenly across the domain, and the headline RMSE hides that. Sort the folds by distance to the nearest remaining carbon sample and the pattern is immediate: at sites with a neighbour within 100 m the two estimators are almost indistinguishable, because the primary data already pin the value down and the secondary is screened out. The whole of the 17 per cent comes from the folds in the sparse half of the survey, where the nearest carbon sample is 250 m or more away and ordinary kriging has fallen back on the global mean. That is the general rule for cokriging: it helps exactly where the primary is thin, which is where you needed help.

Watch the weights themselves for the sign of trouble that the summary statistics hide. Primary weights should mostly lie between 0 and 1 with a few small negatives from the screening effect, and secondary weights should be small and sum to zero. If np.abs(w[:n1]).max() exceeds about 2, or the secondary weights are of the same magnitude as the primary ones, the coregionalisation is close to singular even if it passed the eigenvalue check.

Critical Best Practices #

Fit the shared range once and pin it everywhere #

A linear model of coregionalisation is not three variograms that happen to be plotted together. Every direct and cross variogram must be a weighted sum of the same basic structures with the same ranges; only the coefficients are free. The practical recipe is to fit the best-informed variogram first — usually the dense secondary, or the primary if it is the one you care about most — round its range, and then pass len_scale=RANGE to every subsequent fit_variogram call so gstools treats it as fixed rather than as an initial guess.

Clip the cross coefficient to the Cauchy–Schwarz bound, never pad the diagonal #

When the independently fitted cross sill breaks $|b_{12}| \le \sqrt{b_{11}b_{22}}$, the temptation is to add a small amount to the direct sills until the determinant turns positive. Resist it: those direct sills were estimated from data and inflating them corrupts the univariate model as well. Clip the cross coefficient down to the bound instead, or refit all three variograms jointly with the constraint imposed. Clipping loses a little cross-covariance; padding invents variance that is not there.

Set exact=True when comparing against gs.krige.Ordinary #

gstools kriging defaults to exact=False with cond_err="nugget", which treats the nugget as measurement error and filters it — the predictor no longer honours the data values at the sample locations. A hand-built cokriging system that puts the nugget into the right-hand side does not filter it. Comparing the two without setting exact=True compares two different estimands, and the difference can be larger than the cokriging gain you are trying to measure. Choose one convention and apply it to both.

Use a search neighbourhood before the secondary set gets large #

The system is $(n_1 + n_2 + 2)$ square and is solved once per target location. At 520 samples that is milliseconds; at 20,000 secondary readings a dense solve is roughly $2.7 \times 10^{12}$ operations per target and the workflow stops being viable. Restrict the secondary to the nearest few hundred readings with a scipy.spatial.cKDTree query, or drop to collocated cokriging, which keeps only the secondary value at the target and reduces the system to $n_1 + 2$. The trade-off is set out in Collocated Cokriging vs Regression Kriging.

Keep the cross-variogram’s sign, and check the units #

A negatively correlated secondary — elevation against temperature, say — produces a cross-variogram that runs negative. fit_variogram cannot return a negative var, so fit the model to -g12 and negate the recovered coefficients, or fit np.abs(g12) and restore the sign from np.sign(np.nanmean(g12)). The cross coefficients carry mixed units, per cent times metres here, so a coefficient that looks implausibly large is often just a unit mismatch rather than a modelling error.

Troubleshooting #

Next Steps #

Once the coregionalisation is validated, decide whether the full system is worth its cost against the cheaper alternatives in Collocated Cokriging vs Regression Kriging, and revisit the univariate baseline in Step-by-Step Ordinary Kriging with PyKrige so the comparison you report is against a properly tuned single-variable model rather than a straw man.

Frequently Asked Questions #

Does gstools have a built-in cokriging function?

No. As of the 1.5 and 1.6 releases gstools ships univariate kriging classes only: Simple, Ordinary, Universal, ExtDrifted and Detrended. There is no multivariate estimator and no cross-variogram estimator. What gstools does give you is a set of validated covariance models whose covariance method evaluates the correlation function at arbitrary lags, and that is enough. You build the block system with numpy, fill each block by calling covariance on the right pair of location sets, and solve it yourself.

When is cokriging worth the extra work over ordinary kriging?

When the primary variable is genuinely undersampled relative to the secondary and the cross-correlation is strong, roughly above 0.6 in absolute value. In the worked example here, 60 carbon samples against 460 elevation readings with a marginal correlation of 0.78 buys a 17 per cent reduction in held-out RMSE. If the two variables are sampled at the same density, cokriging reduces almost exactly to ordinary kriging of the primary and the extra machinery earns nothing.

What do I do if the coregionalisation matrices are not positive semi-definite?

Shrink the cross coefficient rather than inflating the direct sills. For each basic structure the Cauchy–Schwarz bound requires the absolute cross coefficient to be no larger than the geometric mean of the two direct coefficients. Fitting the cross-variogram independently routinely breaches it by a few per cent. Clip the cross coefficient to the bound, or refit all three variograms jointly with the constraint imposed. Never add a small ridge to the diagonal and carry on, because that silently changes the direct variograms you fitted.

Can the primary and secondary variables have different ranges?

Not within a single basic structure. A linear model of coregionalisation is a sum of structures, and every variogram in the model must use the same set of shapes and ranges; only the coefficients differ. To give the secondary shorter-range detail, add a second basic structure with its own range and let the primary take a coefficient of zero for it. That is legitimate and keeps the model valid, whereas assigning two different ranges to the same structure does not.

Related

• Collocated Cokriging vs Regression Kriging — the two cheaper ways to use the same secondary variable
• Step-by-Step Ordinary Kriging with PyKrige — the univariate baseline every cokriging result must beat
• Fitting Spherical, Exponential & Gaussian Variogram Models — choosing the basic structures the coregionalisation is built from

← Back to Cokriging & Multivariate Interpolation