Indicator Kriging for Threshold Exceedance Probability

TL;DR: Recode the measurements as ind = (pb > 200).astype(float), fit a spherical variogram to that 0/1 field with skgstat.Variogram, and pass its range, partial sill and nugget to pykrige.ok.OrdinaryKriging(..., variogram_parameters=[psill, vrange, nugget]). Clip the result with np.clip(z, 0, 1). The surface is a probability only once a held-out reliability check says so.

Why This Matters #

Most kriging questions are asked in the wrong units. A regulator does not want to know the estimated lead concentration at a node; they want to know how likely it is that the soil there breaches 200 mg/kg. Those are different questions, and the second cannot be answered by thresholding a kriged concentration surface. Kriging smooths, so a kriged concentration map understates the frequency of high values, and a hard cut applied to it will systematically under-flag the site. Indicator kriging sidesteps the problem by moving the threshold to the front of the pipeline: transform first, then interpolate the transformed variable.

This page works one cut-off end to end. Choosing the cut-off in the first place, and deciding how many of them to model, belongs to choosing thresholds and indicator variograms, and the wider family of methods sits in indicator and probability kriging. The mechanics of the kriging system itself are unchanged from step-by-step ordinary kriging with PyKrige; only the variable being kriged is different.

Define the indicator at cut-off $z_c$ in the exceedance direction:

$$i(\mathbf{s}; z_c) = \begin{cases} 1 & \text{if } Z(\mathbf{s}) > z_c \\ 0 & \text{otherwise.} \end{cases}$$

Because $i$ takes only the values 0 and 1, its expectation collapses to a probability:

$$E[\,i(\mathbf{s}; z_c)\,] = 1 \cdot P(Z(\mathbf{s}) > z_c) + 0 \cdot P(Z(\mathbf{s}) \le z_c) = P(Z(\mathbf{s}) > z_c).$$

Kriging estimates a conditional expectation, so the kriged indicator $\hat{\imath}(\mathbf{s}_0) = \sum_\alpha \lambda_\alpha\, i(\mathbf{s}_\alpha)$ is an estimate of $P(Z(\mathbf{s}_0) > z_c \mid \text{data})$. That single line is the whole justification for reading the output map as a probability surface. Note also that the classical geostatistical convention runs the other way, with $i = 1$ below the cut-off so that the result is a conditional cumulative distribution; the two conventions are complements, and mixing them silently inverts every map you produce.

Environment and Version Pinning #

Three libraries carry the workflow: gstools to simulate a realistic field so the example runs without proprietary data, scikit-gstat to fit the indicator variogram, and pykrige to solve the kriging system.

pip install "numpy>=1.26,<2.2" "scipy>=1.11" "pandas>=2.1" \
            "scikit-gstat==1.0.18" "pykrige==1.7.2" "gstools==1.6.0" \
            "matplotlib>=3.8"

import numpy as np
import pandas as pd
import skgstat as skg
import gstools as gs
from pykrige.ok import OrdinaryKriging

Step-by-Step Implementation #

1. Build the sample dataset and split off a validation set #

The field is simulated as lognormal with an exponential covariance, which is a fair caricature of soil contamination: strongly skewed, spatially continuous at a few hundred metres. Sixty of the 220 samples are held back before anything is fitted, so the reliability check later is honest.

rng = np.random.default_rng(20260807)

SIDE = 2000.0          # metres, square site
N = 220
Z_CUT = 200.0          # mg/kg regulatory cut-off

xs = rng.uniform(0.0, SIDE, N)
ys = rng.uniform(0.0, SIDE, N)

# Lognormal soil lead with an exponential spatial structure.
model = gs.Exponential(dim=2, var=0.55, len_scale=260.0)
srf = gs.SRF(model, mean=np.log(170.0), seed=20260807)
pb = np.exp(srf((xs, ys)))

perm = rng.permutation(N)
train, test = perm[:160], perm[160:]

ind = (pb > Z_CUT).astype(float)      # the indicator transform

print(f"n samples             : {N}")
print(f"median Pb (mg/kg)     : {np.median(pb):.1f}")
print(f"above {Z_CUT:.0f} mg/kg       : {int(ind.sum())}  ({ind.mean():.3f})")
print(f"train / test split    : {len(train)} / {len(test)}")
print(f"train exceedance rate : {ind[train].mean():.3f}")
print(f"test  exceedance rate : {ind[test].mean():.3f}")

n samples             : 220
median Pb (mg/kg)     : 168.4
above 200 mg/kg       : 93  (0.423)
train / test split    : 160 / 60
train exceedance rate : 0.425
test  exceedance rate : 0.417

2. Fit the indicator variogram #

The variogram is fitted to the 0/1 field, not to the concentrations. Its interpretation is unusually concrete: the squared difference of two indicators is 1 exactly when one sample is above the cut-off and the other below, so

$$\gamma_I(\mathbf{h}) = \tfrac{1}{2}E\!\left[\big(i(\mathbf{s}+\mathbf{h}) - i(\mathbf{s})\big)^2\right] = \tfrac{1}{2}\,P\!\left(\text{the pair straddles } z_c\right).$$

At large separations the pair straddles independently, giving a sill of $p(1-p)$ where $p$ is the exceedance rate. That is a free, non-negotiable check on the fit — see fitting empirical variograms with SciKit-GStat for the estimator and binning options behind it.

coords_tr = np.column_stack([xs[train], ys[train]])

V = skg.Variogram(
    coords_tr,
    ind[train],
    model="spherical",
    estimator="matheron",
    bin_func="even",
    n_lags=14,
    maxlag=1000.0,
    normalize=False,
)

vrange, psill, nugget = V.parameters   # scikit-gstat order: range, sill, nugget
p_hat = ind[train].mean()

print(f"range  (m)   : {vrange:.1f}")
print(f"sill         : {psill:.4f}")
print(f"nugget       : {nugget:.4f}")
print(f"sill+nugget  : {psill + nugget:.4f}   vs  p(1-p) = {p_hat*(1-p_hat):.4f}")

range  (m)   : 482.6
sill         : 0.2076
nugget       : 0.0313
sill+nugget  : 0.2389   vs  p(1-p) = 0.2444

The total, 0.2389, sits just under the theoretical 0.2444 — normal, because a finite sample under-estimates the sill of a spatially correlated field. If instead the reported sill alone had landed on 0.2444, your version of scikit-gstat is giving the full sill rather than the partial sill, and you must subtract the nugget before step 3.

3. Krige the indicator onto a grid and clip #

PyKrige takes the variogram parameters as [psill, range, nugget] — a different order from the one scikit-gstat reports, which is the single most common source of a nonsense surface here.

gridx = np.arange(10.0, SIDE, 20.0)    # 100 nodes
gridy = np.arange(10.0, SIDE, 20.0)    # 100 nodes

ok_ind = OrdinaryKriging(
    xs[train], ys[train], ind[train],
    variogram_model="spherical",
    variogram_parameters=[psill, vrange, nugget],   # PyKrige order
    exact_values=True,
    pseudo_inv=True,
    verbose=False,
    enable_plotting=False,
)

z_raw, kvar = ok_ind.execute("grid", gridx, gridy)
z_raw = np.asarray(z_raw)

n_out = int(((z_raw < 0.0) | (z_raw > 1.0)).sum())
prob = np.clip(z_raw, 0.0, 1.0)

print(f"grid                 : {len(gridx)} x {len(gridy)} nodes at 20 m")
print(f"raw kriged indicator : min {z_raw.min():+.4f}  max {z_raw.max():+.4f}")
print(f"nodes outside [0, 1] : {n_out}  ({100 * n_out / z_raw.size:.2f}%)")
print(f"after clipping       : min {prob.min():.4f}  max {prob.max():.4f}")
print(f"mean kriging variance: {np.asarray(kvar).mean():.4f}")

grid                 : 100 x 100 nodes at 20 m
raw kriged indicator : min -0.0613  max +1.0428
nodes outside [0, 1] : 183  (1.83%)
after clipping       : min 0.0000  max 1.0000
mean kriging variance: 0.1176

Ordinary kriging is a linear estimator whose weights may go negative, so nothing in the system constrains the answer to $[0,1]$. Here 1.83% of nodes strayed, by at most 0.06 — the ordinary, benign case, and clipping resolves it.

4. Check reliability on the held-out points #

A probability surface is only a probability surface if it is calibrated. The test is direct: predict at the 60 withheld locations, sort the predictions into bands, and ask whether roughly 30% of the points predicted at around 0.3 actually exceeded the cut-off.

p_test, _ = ok_ind.execute("points", xs[test], ys[test])
p_test = np.clip(np.asarray(p_test), 0.0, 1.0)
obs = ind[test]

edges = np.array([0.0, 0.2, 0.4, 0.6, 0.8, 1.0])
which = np.clip(np.digitize(p_test, edges[1:-1]), 0, 4)

rows = []
for k in range(5):
    m = which == k
    if not m.any():
        continue
    rows.append({
        "bin": f"{edges[k]:.1f}-{edges[k+1]:.1f}",
        "n": int(m.sum()),
        "mean_p": p_test[m].mean(),
        "observed": obs[m].mean(),
    })

rel = pd.DataFrame(rows)
rel["gap"] = rel["observed"] - rel["mean_p"]
print(rel.to_string(index=False, float_format=lambda v: f"{v:.3f}"))

brier = float(np.mean((p_test - obs) ** 2))
base = float(ind[train].mean())                 # climatological forecast
brier_ref = float(np.mean((base - obs) ** 2))
print(f"\nBrier score          : {brier:.3f}")
print(f"climatological Brier : {brier_ref:.3f}")
print(f"Brier skill score    : {1 - brier / brier_ref:.3f}")
print(f"mean predicted       : {p_test.mean():.3f}   observed rate: {obs.mean():.3f}")

    bin  n  mean_p  observed    gap
0.0-0.2 19   0.082     0.105  0.023
0.2-0.4 12   0.291     0.250 -0.041
0.4-0.6  9   0.503     0.444 -0.059
0.6-0.8 11   0.688     0.727  0.039
0.8-1.0  9   0.884     0.889  0.005

Brier score          : 0.158
climatological Brier : 0.243
Brier skill score    : 0.350
mean predicted       : 0.418   observed rate: 0.417

Interpreting the Output #

Read the reliability table before the map. Every band sits within 0.06 of its observed frequency, and the gaps alternate in sign rather than all leaning one way, which is what unbiased looks like at this sample size. The overall mean prediction of 0.418 matching the observed base rate of 0.417 confirms there is no global drift. A one-sided pattern — every gap negative, say — would mean the surface is systematically over-stating exceedance, usually because the training set happened to be sited on hot spots or because the nugget was set too low and the estimator has become an interpolator of local ones.

The Brier score $BS = \frac{1}{m}\sum_j (\hat{p}_j - i_j)^2$ is 0.158 against a climatological reference of 0.243, giving a skill score of 0.350. That reference is the score you would get by predicting the training exceedance rate of 0.425 everywhere, so a skill score of zero means the map is worth nothing beyond a single site-wide number, and a negative one means it is worse than that. Values in the 0.2 to 0.5 band are typical of real contamination surveys; anything above 0.7 on a first pass usually means the held-out points were not really held out.

The kriging variance is not the uncertainty of the probability. The ordinary kriging variance

$$\sigma^2_{OK}(\mathbf{s}_0) = \sum_\alpha \lambda_\alpha\, \gamma_I(\mathbf{s}_\alpha - \mathbf{s}_0) + \mu$$

is built from the variogram and the sample geometry alone. No data value enters it, so two nodes with the same neighbour configuration receive the same variance whether their estimated probabilities are 0.05 or 0.50. It is a map of data density, useful for planning where to sample next but not an error bar. The variability of the actual outcome at a node is the Bernoulli variance $\hat{p}(1-\hat{p})$, which is maximal at $\hat{p}=0.5$ and vanishes at the extremes — behaviour the kriging variance does not have. The general treatment is in uncertainty and variance mapping; the indicator case is simply the sharpest example of the distinction.

Turning the probability surface into a decision map #

A probability surface is not yet a decision. Turning it into one requires a stated risk tolerance $\alpha$: flag every node where $\hat{p}(\mathbf{s}) > \alpha$. The tolerance is not arbitrary. If a false alarm costs $C_{FP}$ (needless excavation) and a miss costs $C_{FN}$ (contaminated soil left in place), acting is worth it when $(1-p)\,C_{FP} < p\,C_{FN}$, which rearranges to

$$\alpha^{*} = \frac{C_{FP}}{C_{FP} + C_{FN}}.$$

A miss that costs nine times a false alarm gives $\alpha^{*} = 0.1$. Report the flagged area at two or three tolerances so the reader sees how sensitive the answer is to a judgement they, not you, should be making.

cell_area_km2 = (20.0 * 20.0) / 1e6

print(f"{'tolerance':>10}  {'nodes':>7}  {'area (km2)':>10}  {'share':>6}")
for alpha in (0.30, 0.50, 0.70):
    flagged = int((prob > alpha).sum())
    print(f"{alpha:>10.2f}  {flagged:>7d}  {flagged * cell_area_km2:>10.3f}  "
          f"{100 * flagged / prob.size:>5.1f}%")

 tolerance    nodes  area (km2)   share
      0.30     3142       1.257   31.4%
      0.50     2015       0.806   20.2%
      0.70     1188       0.475   11.9%

Critical Best Practices #

Anchor the indicator sill to $p(1-p)$ #

The sill of an indicator variogram is not free. It must approach $p(1-p)$, where $p$ is the sample exceedance rate, because that is the variance of a Bernoulli variable with that mean. A fit whose total sill is 0.40 when $p = 0.42$ is wrong regardless of how good the least-squares residual looks; the usual cause is a maxlag extending past the point where the empirical variogram becomes noisy, dragging the model up. Set maxlag to roughly half the site diagonal and check the arithmetic every time.

Watch the parameter order between the two libraries #

scikit-gstat reports [range, sill, nugget]; pykrige expects [psill, range, nugget]. Swapping the first two is silent — kriging runs happily with a partial sill of 483 and a range of 0.21 — and produces a surface that is nearly constant at the global mean, because with a range far shorter than the sample spacing every weight collapses to $1/n$. If your probability map looks suspiciously flat and close to the base rate, check the ordering first.

Do not read the kriging variance as probability uncertainty #

Covered above but worth repeating, because it is the mistake most often shipped in a report. Publish the kriging variance under its own heading, labelled as a measure of sample support, and never as a confidence band around $\hat{p}$. If you need a genuine interval on the probability, you need sequential indicator simulation, not the kriging variance.

Clip, but log what you clipped #

np.clip is the right response to excursions outside the unit interval, and it is also an excellent diagnostic that costs one line. Record the count and the maximum excursion. Under about 5% of nodes and excursions under 0.1 is routine. Beyond that, the variogram is usually the culprit — a nugget too small relative to the true short-scale variability, or a range far longer than the data support — and clipping is papering over it.

Fix the exceedance convention once and state it #

Half the literature defines the indicator as $1$ below the cut-off, so the kriged field is a conditional CDF and the exceedance probability is its complement. Both are correct; mixing them is not. Put the convention in the docstring, and add an assertion that the mean of the kriged surface is within a few hundredths of the sample exceedance rate. That single check catches an inverted map instantly, and it is the kind of assertion that will actually fire one day.

Troubleshooting #

Next Steps #

One cut-off answers one question. To build a full conditional distribution at every node, repeat this workflow across several cut-offs and enforce order relations between them, as set out in choosing thresholds and indicator variograms, then pair the result with uncertainty and variance mapping so the sample-support map is published beside the probability map.

Frequently Asked Questions #

Why is the kriged indicator a probability rather than just an interpolated number?

Because the indicator only takes the values 0 and 1, its expectation is exactly the probability that it equals 1. Kriging produces a best linear unbiased estimate of the conditional expectation, so the estimate at an unsampled location is an estimate of the probability that the underlying value exceeds the cut-off there. Nothing about the arithmetic changes from ordinary kriging; what changes is that the quantity being averaged is a Bernoulli variable, which is what makes the result readable on a probability scale.

Does the kriging variance tell me how uncertain the exceedance probability is?

No. The ordinary kriging variance is a function of the sample geometry and the variogram alone; it never sees the data values, so two locations with identical neighbour configurations get the same variance whether the estimated probability is 0.05 or 0.50. What it measures is how well supported the estimate is by nearby samples. The uncertainty of the underlying outcome at a node is the Bernoulli variance of the estimated probability itself. Use the kriging variance as a data-density map, not as an error bar on the probability.

What should I do about kriged indicator values below 0 or above 1?

Clip them, and record how many nodes were affected and by how much. Ordinary kriging is a linear estimator with weights that may be negative, so nothing constrains it to the unit interval. Small excursions of a few hundredths on one or two per cent of nodes are ordinary and clipping costs nothing. Excursions beyond about 0.1, or affecting more than roughly five per cent of nodes, are a symptom: usually a nugget set too low, a range that is far too long, or a search neighbourhood that is too small.

How many held-out points do I need for the reliability check to mean anything?

Enough that each probability band holds roughly ten observations, so about fifty to sixty points for five bands, and considerably more if the exceedance rate is far from one half. With fewer, the observed frequency in a band swings wildly on one or two outcomes and the diagram tells you nothing. If the survey is too small to hold points back, use leave-one-out cross-validation across the whole dataset instead, and report the number of points behind every band alongside the frequency.

Related

• Choosing Thresholds and Indicator Variograms — how many cut-offs to model and where to put them
• Step-by-Step Ordinary Kriging with PyKrige — the kriging system that indicator kriging reuses unchanged
• Fitting Empirical Variograms with SciKit-GStat — binning, estimators and fit diagnostics for the indicator variogram

← Back to Indicator & Probability Kriging