Indicator & Probability Kriging

Some questions do not have a number as their answer. When a regulator asks whether arsenic in a well exceeds 10 µg/L, or a mine planner asks whether a block sits above cut-off grade, the useful output is not a predicted concentration but a probability of exceedance attached to every location on the map. Indicator kriging reaches that probability by a change of variable: recode every observation as a 0 or a 1 against the threshold, model the variogram of the recoded data, and krige it. Because the mean of a binary variable is a probability, the kriged indicator estimates the conditional probability at that threshold directly, with no assumption that the values are Gaussian. This page, part of Kriging, Interpolation & Surface Generation Techniques, works through the transform, the indicator variogram whose sill is fixed by construction, the order-relation violations that appear as soon as you use more than one threshold, and the cheaper and richer variants either side of the standard method.

Prerequisites #

• Python 3.10 or newer
• numpy>=1.24, geopandas>=1.0, gstools>=1.5, pykrige>=1.7, pandas>=2.0
• A projected CRS in metres, so variogram ranges and search radii are lengths and not degrees
• Point observations (x, y, value) of a continuous variable with a meaningful threshold, at least 100 samples and ideally 200 or more, so each indicator has enough of both classes to support a variogram
• Comfort with the ordinary kriging system — see Ordinary & Universal Kriging, since indicator kriging is that same system run on transformed data
• A variogram fitting workflow you trust — see Theoretical Variogram Models for the permissible families

Mathematical Core #

The indicator transform #

Let $Z(\mathbf{s})$ be the regionalised variable — arsenic concentration, ore grade, rainfall — observed at $n$ locations $\mathbf{s}_1, \dots, \mathbf{s}_n$. Fix a threshold $z_k$. The indicator transform at that threshold replaces each observation with a binary code:

$$i(\mathbf{s}_\alpha; z_k) = \begin{cases} 1 & \text{if } Z(\mathbf{s}_\alpha) \le z_k \\ 0 & \text{otherwise} \end{cases}$$

The convention matters and is worth stating in every report. Coding $1$ for at or below makes the transformed variable an estimator of the cumulative distribution function, as in the geostatistical literature and throughout this page; the exceedance probability you actually want is then $1 - \hat{F}$. Some software codes $1$ for above instead. The mathematics is identical and only the sign of the interpretation flips, but mixing the two conventions inside one pipeline is the commonest way to publish a map that is exactly wrong.

Why a kriged indicator estimates a probability #

Take the expectation of the transformed variable at an arbitrary location. Because $i$ takes only the values $0$ and $1$,

$$E\!\left[ I(\mathbf{s}; z_k) \right] = 1 \cdot P\{Z(\mathbf{s}) \le z_k\} + 0 \cdot P\{Z(\mathbf{s}) > z_k\} = F(z_k),$$

the marginal distribution function evaluated at the threshold. Conditioning on the $n$ observations does not change the argument: the conditional expectation of a binary variable is the conditional probability that it equals one. Kriging estimates a conditional expectation, so the ordinary indicator kriging estimate

$$\hat{F}(\mathbf{s}_0; z_k \mid (n)) = \sum_{\alpha=1}^{n} \lambda_\alpha(\mathbf{s}_0; z_k)\, i(\mathbf{s}_\alpha; z_k), \qquad \sum_{\alpha=1}^{n} \lambda_\alpha = 1,$$

is an estimate of $P\{Z(\mathbf{s}_0) \le z_k \mid (n)\}$. The weights $\lambda_\alpha$ come from the usual ordinary kriging system, but built from the indicator variogram $\gamma_I$ rather than the variogram of $Z$:

$$\sum_{\beta=1}^{n} \lambda_\beta \,\gamma_I(\mathbf{s}_\alpha - \mathbf{s}_\beta; z_k) + \mu = \gamma_I(\mathbf{s}_\alpha - \mathbf{s}_0; z_k), \qquad \alpha = 1, \dots, n,$$

with $\mu$ the Lagrange multiplier enforcing unbiasedness. Nothing here is new machinery: every property of ordinary kriging — exactness at data locations, the screening effect, negative weights behind a nearby sample — carries over unchanged, and so do its failure modes.

One caveat is structural. Kriging is the best linear estimator of a conditional expectation and equals it exactly only for a Gaussian field, while an indicator field is binary. Indicator kriging therefore returns a linear approximation to a probability, which is why estimates can leave the unit interval and why multi-threshold results need repair.

The indicator variogram and its fixed sill #

The indicator variogram is the ordinary semivariogram of the transformed variable:

$$\gamma_I(\mathbf{h}; z_k) = \tfrac{1}{2} E\!\left[ \left\{ i(\mathbf{s} + \mathbf{h}; z_k) - i(\mathbf{s}; z_k) \right\}^2 \right].$$

The difference of two indicators is $0$, $+1$ or $-1$, so its square is itself an indicator — of the two locations falling on opposite sides of the threshold. That gives the indicator variogram a reading available to no other variogram:

$$\gamma_I(\mathbf{h}; z_k) = \tfrac{1}{2} P\!\left\{ \text{exactly one of } Z(\mathbf{s}),\, Z(\mathbf{s} + \mathbf{h}) \text{ lies at or below } z_k \right\}.$$

At separations beyond the range the two locations are independent, so the semivariance equals the variance of the indicator, and the variance of a Bernoulli variable is not free:

$$\gamma_I(\infty; z_k) = \operatorname{Var}\!\left[ I(\mathbf{s}; z_k) \right] = p_k\,(1 - p_k), \qquad p_k = F(z_k).$$

This is the most useful constraint in the whole method. The sill of every indicator variogram is determined by the marginal proportion at that threshold, which you can count off the sample in one line. It is not a parameter to be fitted, and a fit that lands elsewhere is telling you about the estimator, the lag binning or a trend — not about a genuinely different sill. Since $p(1-p)$ peaks at $0.25$ when $p = 0.5$ and collapses towards zero at extreme thresholds, indicators near the median carry the most structure and those in the tails carry very little.

Several thresholds, and the order relations they break #

One threshold gives one probability map. A set of $K$ thresholds $z_1 < z_2 < \dots < z_K$ estimates the whole conditional distribution at every location, from which you can read any exceedance probability, a conditional median or an interquartile range. The price is that each threshold is a separate kriging problem with its own variogram and its own weights, and nothing in those $K$ independent systems knows the results must together form a distribution function. Two properties that hold for any real distribution can therefore fail:

$$0 \le \hat{F}(\mathbf{s}_0; z_k) \le 1 \qquad \text{and} \qquad \hat{F}(\mathbf{s}_0; z_k) \le \hat{F}(\mathbf{s}_0; z_{k+1}) \ \ \text{for } z_k < z_{k+1}.$$

Failures of either are called order-relation violations. Excursions outside $[0,1]$ come from negative kriging weights, which appear whenever a sample is screened by another closer to the prediction point — a routine feature of clustered layouts, not a bug. Non-monotonicity comes from the fitted variograms differing between thresholds, so the same data configuration produces different weights at $z_k$ and $z_{k+1}$. On the worked example below, 45.8% of grid nodes carry at least one violation across five thresholds.

The standard repair is a pair of monotone passes followed by an average. Write $c_k = \min\{1, \max\{0, \hat{F}(\mathbf{s}_0; z_k)\}\}$ for the clipped estimate. The upward pass takes a running maximum from the lowest threshold up,

$$F^{\uparrow}_1 = c_1, \qquad F^{\uparrow}_k = \max\!\left\{ F^{\uparrow}_{k-1},\, c_k \right\},$$

the downward pass takes a running minimum from the highest threshold down,

$$F^{\downarrow}_K = c_K, \qquad F^{\downarrow}_k = \min\!\left\{ F^{\downarrow}_{k+1},\, c_k \right\},$$

and the corrected distribution is their mean, $F^{*}_k = \tfrac{1}{2}\left(F^{\uparrow}_k + F^{\downarrow}_k\right)$. Both passes are non-decreasing in $k$ and both lie in $[0,1]$, so their average is too, which makes the correction a proof rather than a heuristic. Averaging avoids the bias a single pass introduces: the upward pass alone can only raise probabilities, the downward pass alone can only lower them.

Median indicator kriging #

If you can accept one extra assumption, the cost of $K$ thresholds collapses to the cost of one. Suppose every indicator variogram has the same shape once standardised by its own sill:

$$\gamma_I(\mathbf{h}; z_k) = p_k (1 - p_k)\, \bar{\gamma}(\mathbf{h}) \quad \text{for all } k,$$

with $\bar{\gamma}$ a single standardised model reaching a sill of $1$. This is the mosaic model, sometimes stated as intrinsic coregionalisation of the indicators. Substitute it into the kriging system: every term on both sides carries the same factor $p_k(1-p_k)$, which divides out. The weights $\lambda_\alpha$ are then identical for every threshold, so you factorise the matrix once and reuse it for all $K$ indicator vectors. In practice $\bar{\gamma}$ is taken from the threshold nearest the sample median, where the experimental variogram is best resolved — hence the name.

The saving is close to the full factor $K$: on the worked grid below, five thresholds took 1.84 s with their own variograms and 0.37 s with the median approximation. There is a second, less obvious benefit. Because the weights no longer depend on $k$ while the indicator data satisfy $i(\mathbf{s}_\alpha; z_k) \le i(\mathbf{s}_\alpha; z_{k+1})$, monotonicity is guaranteed wherever all the weights are non-negative. Negative weights still let violations through — 9.7% of weights were negative here — but the violation rate falls from 45.8% of nodes to 17.1%.

What median indicator kriging cannot represent is destructuring: the tendency of indicator variograms at extreme thresholds to flatten towards pure nugget, because the few samples above a high cut-off are rarely close enough together to reveal structure. That is a real feature of the data, and precisely what the mosaic model assumes away.

Probability kriging #

Indicator kriging discards information. A well measuring 9.9 µg/L and one measuring 0.4 µg/L are both coded $1$ at the 10 µg/L threshold, yet they say very different things about the chance that a nearby location exceeds it. Probability kriging recovers some of that by cokriging the indicator with the rank-order transform of the value. Define the uniform transform from the sample ranks,

$$u(\mathbf{s}_\alpha) = \frac{r_\alpha - 0.5}{n}, \qquad r_\alpha = \text{rank of } Z(\mathbf{s}_\alpha) \text{ among the } n \text{ observations},$$

which is scale-free, insensitive to outliers and uniform on $(0,1)$. The estimator adds a second sum:

$$\hat{F}(\mathbf{s}_0; z_k) = \sum_{\alpha=1}^{n} \lambda_\alpha\, i(\mathbf{s}_\alpha; z_k) \; + \; \sum_{\beta=1}^{n} \nu_\beta\, u(\mathbf{s}_\beta),$$

subject to $\sum_\alpha \lambda_\alpha = 1$ and $\sum_\beta \nu_\beta = 0$, the second constraint keeping the secondary variable from shifting the mean. Solving it needs three structural models per threshold — the indicator variogram $\gamma_I(\mathbf{h}; z_k)$, the variogram of the uniform transform $\gamma_U(\mathbf{h})$, and the cross-variogram $\gamma_{IU}(\mathbf{h}; z_k)$ — which must together form a permissible linear model of coregionalisation, so they have to be fitted jointly rather than one at a time. The cokriging system has order $2n + 2$ rather than $n + 1$, so with a neighbourhood of $m$ samples the cost per prediction rises roughly eightfold. Gains are usually modest and occasionally negative; treat probability kriging as a hypothesis to be tested by cross-validated Brier score, never as a default upgrade.

Annotated Implementation #

The example is a synthetic but realistically shaped arsenic dataset: 260 wells over a 20 km square in UTM zone 45N, log-normally distributed with a median around 8 µg/L. The decision threshold is the WHO guideline value of 10 µg/L; the other four thresholds give the estimated distribution enough resolution to be worth reporting.

Load, project, and inspect the sample distribution #

import geopandas as gpd
import numpy as np
import pandas as pd

wells = gpd.read_file("wells_arsenic.gpkg", layer="wells").to_crs("EPSG:32645")
wells = wells[wells["as_ugl"].notna()].reset_index(drop=True)

x = wells.geometry.x.to_numpy()
y = wells.geometry.y.to_numpy()
z = wells["as_ugl"].to_numpy(dtype=float)

print(f"{len(z)} wells over {np.ptp(x)/1000:.1f} x {np.ptp(y)/1000:.1f} km, {wells.crs.to_string()}")
print(f"arsenic  min {z.min():.1f}  median {np.median(z):.1f}  "
      f"mean {z.mean():.1f}  max {z.max():.1f} ug/L")

260 wells over 19.9 x 19.9 km, EPSG:32645
arsenic  min 0.4  median 8.0  mean 10.6  max 61.5 ug/L

The mean sitting well above the median is the usual signature of a skewed environmental variable, and that skew is exactly why threshold work is done on indicators: no transformation of $Z$ is needed, no back-transform bias is incurred, and a single extreme value cannot drag the estimate the way it drags an ordinary kriging prediction.

Apply the indicator transform and count the sills #

THRESHOLDS = np.array([4.0, 7.0, 10.0, 15.0, 25.0])   # ug/L; 10 is the decision threshold

# (n, K) matrix of indicators; column k is the 0/1 coding at THRESHOLDS[k].
# The <= convention makes each column an estimator of F(z_k), not of exceedance.
I = (z[:, None] <= THRESHOLDS[None, :]).astype(float)

p = I.mean(axis=0)              # marginal proportion at or below each threshold
sill = p * (1.0 - p)            # the sill is now known, not a free parameter

summary = pd.DataFrame({
    "z": THRESHOLDS, "n_below": I.sum(axis=0).astype(int),
    "p": p, "1-p": 1 - p, "sill": sill,
})
print(summary.to_string(index=False, float_format=lambda v: f"{v:.4f}"))

    z  n_below      p    1-p    sill
 4.0000      64 0.2462 0.7538  0.1856
 7.0000     120 0.4615 0.5385  0.2485
10.0000     156 0.6000 0.4000  0.2400
15.0000     204 0.7846 0.2154  0.1690
25.0000     240 0.9231 0.0769  0.0710

Read the last column before doing anything else. The 25 µg/L indicator has a sill of 0.0710, less than a third of the 7 µg/L indicator’s 0.2485, because only 20 of the 260 wells sit above it — not enough to resolve a variogram over ten lags, and so the weakest link in the set.

Estimate and fit the indicator variograms with the sill pinned #

import gstools as gs

BIN_EDGES = np.arange(0.0, 9000.0 + 1.0, 900.0)   # 10 lags of 900 m, half the domain

models, experimental = [], []
for k, t in enumerate(THRESHOLDS):
    centre, gamma_k = gs.vario_estimate((x, y), I[:, k], bin_edges=BIN_EDGES)
    experimental.append((centre, gamma_k))

    m = gs.Spherical(dim=2)
    # sill= fixes var + nugget to p(1-p); only the range and the nugget
    # split are actually estimated from the experimental points.
    m.fit_variogram(centre, gamma_k, sill=float(sill[k]), nugget=True)
    models.append(m)

    print(f"z={t:5.1f}  p={p[k]:.3f}  sill={sill[k]:.4f}  "
          f"range={m.len_scale:6.0f} m  nugget={m.nugget:.4f}  "
          f"nugget/sill={m.nugget/sill[k]:.2f}")

z=  4.0  p=0.246  sill=0.1856  range=  2425 m  nugget=0.0742  nugget/sill=0.40
z=  7.0  p=0.462  sill=0.2485  range=  4725 m  nugget=0.1491  nugget/sill=0.60
z= 10.0  p=0.600  sill=0.2400  range=  4000 m  nugget=0.1120  nugget/sill=0.47
z= 15.0  p=0.785  sill=0.1690  range=  3650 m  nugget=0.0789  nugget/sill=0.47
z= 25.0  p=0.923  sill=0.0710  range=  4500 m  nugget=0.0509  nugget/sill=0.72

Two things stand out. The nugget-to-sill ratios are high — 0.40 to 0.72 — which is normal for indicators, not a bad model: the transform throws away magnitude, so two nearby wells straddling the threshold contribute a full unit of squared difference however close their concentrations are. And the 0.72 at 25 µg/L is destructuring made numerical. The experimental variogram there is nearly flat, hovering between 0.058 and 0.076 against a sill of 0.071, so its fitted range of 4500 m carries almost no information. Treat any threshold whose nugget exceeds about two-thirds of its sill as decorative, and say so when you report the map — avoiding that is the subject of Choosing Thresholds and Indicator Variograms.

Krige each indicator onto the grid #

from pykrige.ok import OrdinaryKriging

gridx = np.arange(0.0, 20000.0, 100.0) + 50.0      # 200 x 200 nodes at 100 m
gridy = gridx.copy()

def krige_indicator(k):
    m = models[k]
    ok = OrdinaryKriging(
        x, y, I[:, k],
        variogram_model="spherical",
        # pykrige reads 'sill' as the total including the nugget
        variogram_parameters={"sill": float(sill[k]),
                              "range": float(m.len_scale),
                              "nugget": float(m.nugget)},
        exact_values=True,
        pseudo_inv=True,          # tolerate near-singular local systems
    )
    # a moving neighbourhood requires backend="loop" in pykrige
    est, var = ok.execute("grid", gridx, gridy,
                          backend="loop", n_closest_points=24)
    return np.asarray(est).ravel()

F = np.column_stack([krige_indicator(k) for k in range(len(THRESHOLDS))])
print("F shape:", F.shape)
print(f"z = 10 estimates span [{F[:, 2].min():.4f}, {F[:, 2].max():.4f}]")

F shape: (40000, 5)
z = 10 estimates span [-0.0043, 1.0501]

The span is the first diagnostic and it already fails: $1.0501$ is not a probability. Nothing has gone wrong — that is a linear estimator doing what linear estimators do with negative weights. The repair comes next.

Count and repair the order-relation violations #

def order_relation_report(F, label):
    oob = (F < 0) | (F > 1)
    dec = np.diff(F, axis=1) < -1e-9
    print(f"{label}")
    print(f"  estimates outside [0, 1] : {oob.sum():6d} / {F.size:6d} "
          f"({100 * oob.mean():.2f}%)")
    print(f"  decreasing steps         : {dec.sum():6d} / {dec.size:6d} "
          f"({100 * dec.mean():.2f}%)")
    print(f"  nodes with any violation : {(oob.any(1) | dec.any(1)).sum():6d} / "
          f"{len(F):6d} ({100 * (oob.any(1) | dec.any(1)).mean():.1f}%)")


def fix_order_relations(F):
    """Clip to [0, 1], then average an upward and a downward monotone pass."""
    c = np.clip(F, 0.0, 1.0)
    up = np.maximum.accumulate(c, axis=1)                   # running max, low -> high
    down = np.minimum.accumulate(c[:, ::-1], axis=1)[:, ::-1]  # running min, high -> low
    return 0.5 * (up + down)


order_relation_report(F, "raw indicator kriging")
F_ok = fix_order_relations(F)
order_relation_report(F_ok, "after upward/downward averaging")
print(f"  mean |change| = {np.abs(F_ok - F).mean():.4f}   "
      f"max |change| = {np.abs(F_ok - F).max():.4f}")

raw indicator kriging
  estimates outside [0, 1] :   5370 / 200000 (2.69%)
  decreasing steps         :  18469 / 160000 (11.54%)
  nodes with any violation :  18334 /  40000 (45.8%)
after upward/downward averaging
  estimates outside [0, 1] :      0 / 200000 (0.00%)
  decreasing steps         :      0 / 160000 (0.00%)
  nodes with any violation :      0 /  40000 (0.0%)
  mean |change| = 0.0027   max |change| = 0.1437

Note the shape of the correction. The mean absolute change is $0.0027$ while the largest single change is $0.1437$: the repair does nothing at the great majority of nodes and a great deal at the handful where the raw estimate was least trustworthy, which is what you want from a constraint the problem implied all along.

Take the cheap route with median indicator kriging #

K_MEDIAN = int(np.argmin(np.abs(p - 0.5)))     # threshold closest to the sample median
print(f"median threshold: z = {THRESHOLDS[K_MEDIAN]:.1f} ug/L (p = {p[K_MEDIAN]:.3f})")

m = models[K_MEDIAN]
ok = OrdinaryKriging(
    x, y, I[:, K_MEDIAN],
    variogram_model="spherical",
    variogram_parameters={"sill": float(sill[K_MEDIAN]),
                          "range": float(m.len_scale),
                          "nugget": float(m.nugget)},
    exact_values=True, pseudo_inv=True,
)

# The weights depend only on the geometry and the variogram, both now fixed,
# so the same solver object serves every threshold: only the data vector changes.
F_med = np.empty_like(F)
for k in range(len(THRESHOLDS)):
    ok.Z = I[:, k]                      # same weights, different indicator data
    est, _ = ok.execute("grid", gridx, gridy, backend="loop", n_closest_points=24)
    F_med[:, k] = np.asarray(est).ravel()

order_relation_report(F_med, "median indicator kriging")
print(f"  mean |F_med - F| at z = 10 : {np.abs(F_med[:, 2] - F[:, 2]).mean():.4f}")
print(f"  max  |F_med - F| at z = 10 : {np.abs(F_med[:, 2] - F[:, 2]).max():.4f}")

median indicator kriging
  estimates outside [0, 1] :   4595 / 200000 (2.30%)
  decreasing steps         :   3329 / 160000 (2.08%)
  nodes with any violation :   6825 /  40000 (17.1%)
  mean |F_med - F| at z = 10 : 0.0379
  max  |F_med - F| at z = 10 : 0.1262

Decreasing steps fall from 11.54% to 2.08%, because a common set of weights applied to nested indicator data can only break monotonicity where a weight is negative. Excursions outside the unit interval barely move, 2.30% against 2.69%, since those come from the same negative weights either way. The two surfaces differ by 0.038 on average at the decision threshold and by at most 0.126, and the mapped area above a 0.5 exceedance probability comes out at 147.9 km² instead of 151.0 km² — a two per cent difference in the headline number for a fifth of the compute.

Threshold Selection and Variogram Configuration #

Three decisions determine whether the exceedance map is worth publishing, and only the first usually gets any thought.

How many thresholds, and where. The decision threshold is fixed by the regulation or the cut-off. The rest exist to give the estimated distribution shape, and nine thresholds at the deciles is a reasonable default when $n$ is a few hundred. The binding constraint is the minority class count: an indicator with fewer than about 25 samples in the smaller class cannot support a variogram, so drop thresholds beyond the ninth decile rather than fitting noise.

def usable_thresholds(z, candidates, min_minority=25):
    """Keep only thresholds with enough samples on BOTH sides."""
    keep = []
    for t in candidates:
        n_below = int((z <= t).sum())
        minority = min(n_below, len(z) - n_below)
        flag = "keep" if minority >= min_minority else "DROP"
        print(f"  z={t:6.1f}  below={n_below:4d}  minority={minority:4d}  {flag}")
        if minority >= min_minority:
            keep.append(t)
    return np.array(keep)

deciles = np.quantile(z, np.arange(0.1, 0.91, 0.1))
candidates = np.unique(np.round(np.r_[deciles, 10.0], 1))
THRESHOLDS = usable_thresholds(z, candidates)

The lag binning. Indicator variograms are noisier than the variogram of $Z$. Use fewer, wider lags than you would for the raw variable — ten out to half the domain diagonal is a sensible start — and check the pair count in every bin. Under about 30 pairs, the experimental point is decoration.

The nugget split. With the total sill pinned to $p(1-p)$, the fit has only two free parameters: the range, and the share of the sill taken by the nugget. Always fit the nugget rather than forcing it to zero. Indicator nuggets are genuinely large, and forcing continuity at the origin hands almost all the weight to the single nearest sample and produces a map of Voronoi cells.

for k, t in enumerate(THRESHOLDS):
    centre, gamma_k = gs.vario_estimate((x, y), I[:, k], bin_edges=BIN_EDGES)
    m = gs.Spherical(dim=2)
    m.fit_variogram(centre, gamma_k, sill=float(sill[k]), nugget=True)
    ratio = m.nugget / sill[k]
    if ratio > 0.66:
        print(f"  z={t:5.1f}: nugget/sill={ratio:.2f} — near-pure nugget, "
              f"treat this threshold's map as the marginal proportion")
    if m.len_scale > 0.5 * np.hypot(np.ptp(x), np.ptp(y)):
        print(f"  z={t:5.1f}: range {m.len_scale:.0f} m exceeds half the domain — "
              f"unresolved, consider a trend model instead")

Anisotropy deserves the same treatment it gets in ordinary kriging: a directional indicator variogram is cheap, and often reveals that contamination follows a palaeochannel or mineralisation a strike direction, which an isotropic model would smear perpendicular to the real structure. Directional estimation and the permissible families are covered under Theoretical Variogram Models.

Output Interpretation #

An exceedance probability map is only useful if its numbers mean what they say. The test is reliability: among all locations where the model predicted a probability near $0.3$, roughly 30% should actually exceed the threshold. Leave-one-out cross-validation gives predicted probabilities where the truth is known, and binning them turns the comparison into a diagram.

obs = (z > 10.0).astype(float)          # observed exceedance at each well
pred = np.empty(len(z))                 # leave-one-out predicted exceedance

m10 = models[2]                         # the 10 ug/L indicator model
params = {"sill": float(sill[2]), "range": float(m10.len_scale),
          "nugget": float(m10.nugget)}

for i in range(len(z)):
    keep = np.delete(np.arange(len(z)), i)
    ok_i = OrdinaryKriging(x[keep], y[keep], I[keep, 2],
                           variogram_model="spherical",
                           variogram_parameters=params,
                           exact_values=True, pseudo_inv=True)
    est, _ = ok_i.execute("points", np.array([x[i]]), np.array([y[i]]),
                          backend="loop", n_closest_points=24)
    pred[i] = 1.0 - float(np.asarray(est)[0])       # exceedance = 1 - F

pred = np.clip(pred, 0.0, 1.0)
edges = np.array([0.0, 0.1, 0.2, 0.4, 0.6, 0.8, 1.0 + 1e-9])
for lo, hi in zip(edges[:-1], edges[1:]):
    sel = (pred >= lo) & (pred < hi)
    print(f"  [{lo:.1f}, {min(hi, 1.0):.1f})  n={sel.sum():3d}  "
          f"mean predicted={pred[sel].mean():.3f}  observed={obs[sel].mean():.3f}")

brier = float(np.mean((pred - obs) ** 2))
reference = float(np.mean((obs.mean() - obs) ** 2))
print(f"  Brier {brier:.4f}   base-rate reference {reference:.4f}   "
      f"skill {1 - brier / reference:.3f}")

  [0.0, 0.1)  n= 48  mean predicted=0.039  observed=0.042
  [0.1, 0.2)  n= 36  mean predicted=0.148  observed=0.194
  [0.2, 0.4)  n= 50  mean predicted=0.287  observed=0.280
  [0.4, 0.6)  n= 52  mean predicted=0.499  observed=0.500
  [0.6, 0.8)  n= 45  mean predicted=0.698  observed=0.667
  [0.8, 1.0)  n= 29  mean predicted=0.883  observed=0.862
  Brier 0.1645   base-rate reference 0.2400   skill 0.315

Five of the six bins land within about 0.03 of their observed frequency, which is as close as 30 to 50 wells per bin can resolve. The exception is the $0.1$–$0.2$ bin, where the model predicted 0.148 and 19.4% of wells exceeded: a mild under-forecast, inside the sampling noise for 36 wells. The Brier score of 0.1645 against the base-rate reference of 0.2400 gives a skill of 0.315, so the map removes about a third of the squared error you would incur by assigning every location the sample exceedance rate of 0.400.

What good looks like: bins hugging the diagonal, a positive skill score, a mean predicted probability across the grid close to the sample proportion (here 0.402 against 0.400), and a spread wide enough to distinguish locations. Warning signs, in order of severity: a skill score at or below zero, meaning the map is no better than the base rate; every probability bunched between 0.35 and 0.45, meaning the nugget has swallowed the structure and the weights have flattened towards $1/m$; and systematic departure from the diagonal in one direction, which usually means the sill was fitted rather than pinned. The cross-validation loop and the Brier decomposition are worked through in Indicator Kriging for Threshold Exceedance Probability.

Read the map itself in areas, not pixels. Here 151.0 km² of the 400 km² domain carry an exceedance probability above 0.5, 50.8 km² above 0.8 and 19.2 km² above 0.9 — three numbers that serve a decision-maker better than any colour ramp, and the ones that move when you swap in median indicator kriging, which gave 147.9 km². The kriging variance returned with each estimate is not a variance of the probability, since the indicator is not Gaussian; use it as a data-configuration index, as in Uncertainty & Variance Mapping, not as an error bar.

Production Considerations #

Complexity in the number of thresholds. With a global neighbourhood, one indicator costs $O(n^3)$ to factorise and $O(n^2)$ per prediction location; $K$ thresholds cost $K$ times that, because each has its own matrix. With a moving neighbourhood of $m$ samples — the normal choice, and the only one that scales — each prediction costs $O(m^3)$ and the grid pass costs $O(K \cdot G \cdot m^3)$ for $G$ nodes. Median indicator kriging drops this to $O(G \cdot m^3 + K \cdot G \cdot m)$, because the solve is shared and only a dot product is repeated. The measured 5.0-fold gain for $K = 5$ is close to the theoretical $K$, and grows with it: at nine deciles the saving is nearly ninefold.

Neighbourhood size dominates. The $m^3$ term means a search of 48 samples costs eight times a search of 24, and for indicators large neighbourhoods buy less than they do for continuous variables, because the high nugget flattens the weights. Start at $m = 24$, check by cross-validation whether $m = 40$ moves the Brier score at all, and stop there.

Memory and blocking. The output is $G \times K$ doubles: 40,000 nodes at five thresholds is 1.6 MB, but a 5,000 × 5,000 raster at nine deciles is 1.8 GB and must be tiled. Process the grid in row blocks, apply fix_order_relations per block — it is row-wise with no cross-node dependency — and write each block straight to a rasterio dataset with windows.Window.

Parallelisation and caching. The work is embarrassingly parallel two ways: across thresholds, which suits joblib.Parallel(n_jobs=K) when $K$ is large and the grid small, and across grid blocks, which suits dask when the grid is large. Do not use both without pinning BLAS threads — OMP_NUM_THREADS=1 per worker — or nested threading makes the job slower than serial. Cache fitted variogram parameters keyed by a hash of the coordinates, values and threshold list; do not cache weights against a moving neighbourhood.

Determinism. Nothing here is stochastic, so a rerun must be byte-identical. If it is not, the culprit is almost always pseudo_inv=True masking a near-singular local system — a sign of duplicate sample locations that should be collapsed first.

Troubleshooting #

Next Steps #

Two guides go deeper than this page can. Indicator Kriging for Threshold Exceedance Probability takes a single threshold end to end, from the transform through the cross-validated Brier decomposition to a published raster — start there if you have one regulatory limit and no interest in the rest of the distribution. Choosing Thresholds and Indicator Variograms handles the multi-threshold case: how many cut-offs a sample can support, how to bin lags for a binary variable, and how to recognise destructuring before it reaches the map. From there, compare the indicator route against modelling values directly with Ordinary & Universal Kriging.

Frequently Asked Questions #

Does indicator kriging give the probability for a point or for a block?

The estimate is for the same support as the data. If the indicators come from point samples, the result is the probability that a point value at that location exceeds the threshold. A block exceedance probability is a different quantity, because a block average is less variable than a point value and crosses the threshold less often. Averaging point probabilities over a block does not produce it; you need block indicator kriging with a regularised indicator variogram, or a simulation-based approach.

Why do my exceedance probabilities not decrease monotonically as the threshold rises?

Each threshold is kriged as an independent problem with its own variogram and its own weights, so nothing in the system forces the estimates to form a valid distribution function. Negative kriging weights, which arise from screening in clustered sample layouts, also push estimates outside the unit interval. These are called order-relation violations and they are routine, not a sign of a broken pipeline. Repair them with the upward and downward averaging correction before mapping anything.

When is median indicator kriging good enough?

It is a good approximation when the standardised indicator variograms have a similar shape across thresholds, which usually holds for the central part of the distribution and fails at the extremes, where the indicator variogram flattens towards pure nugget. Because one system is solved for every threshold, it runs in a fraction of the time and produces far fewer order-relation violations. Check it by fitting the full set of indicator variograms anyway and comparing their standardised shapes.

Is probability kriging worth the extra cost?

Sometimes. Probability kriging cokriges the indicator with the rank-order transform of the value, so it uses how far a sample sits from the threshold rather than only which side it falls on. That extra information matters most when samples are sparse near the threshold. The price is a cokriging system of roughly twice the order, plus a cross-variogram per threshold that must remain permissible. Test it against indicator kriging by cross-validated Brier score before committing to it.

Related

• Indicator Kriging for Threshold Exceedance Probability — one threshold from transform to validated raster
• Choosing Thresholds and Indicator Variograms — how many cut-offs a sample supports, and how to bin their lags
• Ordinary & Universal Kriging — the system indicator kriging runs on transformed data
• Uncertainty & Variance Mapping — reading the variance surface that comes with every kriged map
• Theoretical Variogram Models — the permissible families an indicator variogram must be fitted from

← Back to Kriging, Interpolation & Surface Generation Techniques