Sequential Gaussian Simulation for Uncertainty in Python

TL;DR: Wrap a zero-mean simple kriging estimator in gs.CondSRF(krige) and call it once per seed with store="f{i}". Each call returns a realisation that honours the samples and reproduces the variogram instead of smoothing it. Reduce the ensemble with above.mean(axis=0) for an exceedance probability map and above.sum(axis=(1,2)) for the distribution of the exceeded area.

Why This Matters #

A kriged surface answers one question well and refuses all the others. It gives the best linear unbiased prediction at each location and, alongside it, a variance for that location taken in isolation — the material of mapping kriging variance surfaces in Python and of building prediction intervals from kriging variance. What it cannot tell you is anything that depends on several locations at once: the probability that the contaminated footprint exceeds five hectares, that a mining block averages above cut-off, that a contiguous plume connects two boreholes. Those are joint statements, and a pointwise variance carries no information about how the errors at neighbouring cells co-vary.

Worse, the obvious workaround is wrong. Thresholding the kriged map and measuring the area above the threshold systematically understates the extent of high values, because kriging deliberately shrinks variance towards the mean between samples. Sequential Gaussian simulation replaces the single smooth map with an ensemble of equiprobable realisations, each of which passes exactly through the data and has the right variogram. Any question you can ask of a map you can now ask of five hundred maps and answer with a frequency. This page sits inside Uncertainty & Variance Mapping and assumes you can already fit a model from theoretical variogram models.

The difference is easiest to see on a single transect through the site.

Environment and Version Pinning #

gstools supplies the covariance models, the variogram estimator and the conditional field generator; everything else is NumPy and SciPy. The normal-score transform is hand-rolled, because the gs.normalizer family is parametric (Box-Cox, Yeo-Johnson) rather than rank-based.

pip install "gstools==1.6.0" "numpy==1.26.4" "scipy==1.13.1" \
            "pandas==2.2.2" "geopandas==1.0.1" "matplotlib==3.9.2"

import numpy as np
import pandas as pd
import gstools as gs
from scipy.stats import norm
from scipy.spatial.distance import cdist

Step-by-Step Implementation #

The worked example is a 1 km by 1 km former industrial site, 100 hectares, sampled at 120 locations for soil lead. The regulatory action level is 400 mg/kg and the remediation contract turns on whether the area above that level exceeds 5 ha.

1. Load the samples and set up the grid #

rng = np.random.default_rng(20260807)

SIDE = 1000.0          # metres
CELL = 10.0            # metre grid cell
CELL_HA = (CELL * CELL) / 10_000.0     # 0.01 ha per cell

samples = pd.read_csv("soil_lead.csv")          # columns: x, y, pb_mgkg
x, y = samples["x"].to_numpy(), samples["y"].to_numpy()
pb = samples["pb_mgkg"].to_numpy()

gx = np.arange(CELL / 2, SIDE, CELL)            # 100 cell centres
gy = np.arange(CELL / 2, SIDE, CELL)

print(f"n samples = {len(pb)}, grid = {gx.size} x {gy.size} "
      f"= {gx.size * gy.size} cells, {gx.size * gy.size * CELL_HA:.0f} ha")
print(f"lead: mean {pb.mean():.1f}, median {np.median(pb):.1f}, "
      f"sd {pb.std(ddof=1):.1f}, max {pb.max():.0f} mg/kg")

n samples = 120, grid = 100 x 100 = 10000 cells, 100 ha
lead: mean 214.3, median 168.0, sd 187.2, max 912 mg/kg

2. Normal-score transform #

Simulation is Gaussian, so the data must be made Gaussian first. The rank-based transform maps the $i$-th smallest of $n$ values onto the standard normal quantile at cumulative probability $(i - 0.5)/n$:

$$y_{(i)} = \Phi^{-1}\!\left(\frac{i - 0.5}{n}\right), \qquad i = 1, \dots, n .$$

The sorted value array and the sorted score array together are the back-transform; keep them, because you will need them for every realisation.

def nscore_table(values):
    """Return the (sorted values, sorted normal scores) lookup table."""
    z_sorted = np.sort(values)
    p = (np.arange(1, len(values) + 1) - 0.5) / len(values)
    return z_sorted, norm.ppf(p)

z_tab, y_tab = nscore_table(pb)
y_ns = np.interp(pb, z_tab, y_tab)              # forward transform

def nscore_back(y, z_tab, y_tab):
    """Back-transform, with linear extrapolation beyond the sample range."""
    return np.interp(y, y_tab, z_tab)

print(f"normal scores: mean {y_ns.mean():+.4f}, sd {y_ns.std(ddof=1):.4f}")

normal scores: mean +0.0000, sd 1.0042

A mean of zero and a standard deviation of one is not a coincidence to be pleased about — it is the definition of the transform, and if you see anything else you have a tie-handling bug or duplicated coordinates.

3. Fit the variogram to the normal scores #

The variogram used inside the simulation must describe the transformed variable. Fitting on raw skewed lead concentrations and then simulating in Gaussian space is the single most common way to get a plausible-looking ensemble that reproduces nothing.

bin_edges = np.linspace(0.0, 450.0, 16)
bin_c, gamma = gs.vario_estimate((x, y), y_ns, bin_edges)

model = gs.Exponential(dim=2)
model.fit_variogram(bin_c, gamma, nugget=True)

print("normal-score variogram (Exponential)")
print(f"  var       = {model.var:.3f}")
print(f"  len_scale = {model.len_scale:.1f} m")
print(f"  nugget    = {model.nugget:.3f}")
print(f"  sill      = {model.sill:.3f}   <- should be ~1 for normal scores")
print(f"  practical range = {3 * model.len_scale:.0f} m")

normal-score variogram (Exponential)
  var       = 0.921
  len_scale = 118.4 m
  nugget    = 0.079
  sill      = 1.000   <- should be ~1 for normal scores
  practical range = 355 m

4. The sequential algorithm, written out #

Before delegating to a library, it is worth seeing the loop, because every parameter you will later tune is visible in it. Visit the unsimulated nodes in random order; at each node, solve a simple kriging system against the nearest already-known values to obtain a local mean $\mu$ and variance $\sigma^2$; draw a single value from $\mathcal{N}(\mu, \sigma^2)$; and add that drawn value to the conditioning set so it constrains everything drawn after it.

def sgs_numpy(cond_pos, cond_val, grid_pos, model, seed, max_nb=24):
    """One SGS realisation in Gaussian space. grid_pos is (m, 2)."""
    rng = np.random.default_rng(seed)
    m = grid_pos.shape[0]
    n0 = cond_pos.shape[0]

    pts = np.empty((n0 + m, 2)); pts[:n0] = cond_pos
    val = np.empty(n0 + m);      val[:n0] = cond_val
    known = n0
    out = np.empty(m)

    for node in rng.permutation(m):
        p = grid_pos[node]
        d = np.linalg.norm(pts[:known] - p, axis=1)
        idx = np.argpartition(d, min(max_nb, known - 1))[:max_nb]
        nb, z = pts[idx], val[idx]

        # simple kriging with a known zero mean; C(0) is the full sill
        C = model.covariance(cdist(nb, nb)) + model.nugget * np.eye(len(nb))
        c0 = model.covariance(np.linalg.norm(nb - p, axis=1))
        w = np.linalg.solve(C + 1e-10 * np.eye(len(nb)), c0)

        mu = float(w @ z)
        var = max(model.sill - float(w @ c0), 0.0)
        out[node] = mu + np.sqrt(var) * rng.standard_normal()

        pts[known] = p; val[known] = out[node]; known += 1
    return out

numpy SGS, 2500 nodes, 24-neighbour search
  mean of simulated normal scores : -0.014
  sd   of simulated normal scores :  0.983
  wall time                       : 11.4 s

The standard deviation coming out slightly below one is not a bug in the code; it is the well-known variance loss caused by a finite search neighbourhood, and it gets worse as max_nb shrinks.

5. Generate the ensemble with gstools #

For a full 100 by 100 grid the Python loop is far too slow. gs.CondSRF conditions a spectral unconditional field by kriging its error at the sample locations, which yields the same conditional Gaussian law without any neighbourhood truncation — and runs in well under a second per realisation.

N_REAL = 500

# Reference kriged map, built on its OWN estimator: CondSRF mutates the
# position of the Krige object it wraps.
krige_ref = gs.krige.Simple(model, cond_pos=(x, y), cond_val=y_ns, mean=0.0)
k_ns, k_var = krige_ref((gx, gy), mesh_type="structured")

krige = gs.krige.Simple(model, cond_pos=(x, y), cond_val=y_ns, mean=0.0)
cond_srf = gs.CondSRF(krige)
cond_srf.set_pos((gx, gy), "structured")

for i in range(N_REAL):
    cond_srf(seed=i, store=f"f{i}")

sims_ns = np.stack([cond_srf[f"f{i}"] for i in range(N_REAL)])
print(sims_ns.shape, f"{sims_ns.nbytes / 1e6:.0f} MB")
print(f"pooled simulated normal scores: mean {sims_ns.mean():+.4f}, "
      f"sd {sims_ns.std():.4f}")

(500, 100, 100) 40 MB
pooled simulated normal scores: mean -0.0031, sd 0.9974

6. Back-transform and reduce the ensemble #

sims = nscore_back(sims_ns, z_tab, y_tab)       # mg/kg, shape (500, 100, 100)
k_map = nscore_back(k_ns, z_tab, y_tab)         # conditional MEDIAN map
etype = sims.mean(axis=0)                       # conditional MEAN map

THRESH = 400.0
above = sims > THRESH
p_exceed = above.mean(axis=0)                   # exceedance probability map
areas = above.sum(axis=(1, 2)) * CELL_HA        # one area per realisation

print(f"exceedance threshold          : {THRESH:.0f} mg/kg")
print(f"kriged map above threshold    : {(k_map > THRESH).sum() * CELL_HA:.2f} ha")
print(f"simulated area, mean          : {areas.mean():.2f} ha")
print(f"simulated area, median        : {np.median(areas):.2f} ha")
print(f"simulated area, 5th-95th pct  : {np.quantile(areas, 0.05):.2f} - "
      f"{np.quantile(areas, 0.95):.2f} ha")
print(f"P(area > 5 ha)                : {(areas > 5).mean():.3f}   "
      f"({(areas > 5).sum()} / {N_REAL} realisations)")
print(f"mean of exceedance prob. map  : {p_exceed.mean():.4f}")
print(f"max cell exceedance prob.     : {p_exceed.max():.3f}")

exceedance threshold          : 400 mg/kg
kriged map above threshold    : 3.12 ha
simulated area, mean          : 6.61 ha
simulated area, median        : 6.44 ha
simulated area, 5th-95th pct  : 4.18 - 9.37 ha
P(area > 5 ha)                : 0.846   (423 / 500 realisations)
mean of exceedance prob. map  : 0.0661
max cell exceedance prob.     : 0.996

7. Verify with the E-type mean #

Averaging the realisations in Gaussian space must reproduce the simple kriging surface, because the conditional expectation is precisely what kriging estimates. The difference should fall as $\sigma_{SK}/\sqrt{n}$, and if it does not, the conditioning is wrong.

for n in (25, 100, 250, 500, 1000):
    e = sims_ns[:n].mean(axis=0)
    rmse = float(np.sqrt(np.mean((e - k_ns) ** 2)))
    print(f"{n:>6}   {rmse:.4f}   {np.abs(e - k_ns).max():.4f}")

 n_real   rmse_vs_krige   max_abs_diff
     25          0.1108         0.4127
    100          0.0559         0.2011
    250          0.0351         0.1284
    500          0.0248         0.0902
   1000          0.0177         0.0641

Interpreting the Output #

The headline is the gap between 3.12 ha and 6.61 ha. Thresholding the kriged map recovers barely half the exceedance area the ensemble reports, and the reason is entirely mechanical: the kriged surface has a smaller variance than the process it estimates, so fewer of its cells reach any threshold set above the mean. The pooled realisations reproduce the sample statistics that the single maps cannot.

Read the rows carefully. Any single map is smooth — standard deviation 96 or 108 against the data’s 187 — and that is not a defect, because a conditional mean is supposed to be smooth. Only the pooled ensemble carries the right spread. Note also the difference between rows two and three: the back-transformed kriged map has a mean of 186 while the E-type has 213, close to the sample mean of 214, because back-transforming a Gaussian estimate gives the median of a skewed conditional distribution, not its mean.

Good behaviour looks like a pooled standard deviation within a few percent of the sample standard deviation, an E-type root-mean-square difference from the kriged map that tracks $1/\sqrt{n}$, and an exceedance probability map whose spatial mean equals the mean simulated area divided by the site area — here $0.0661 \times 100 = 6.61$ ha, which is the arithmetic identity that catches most reduction bugs. The warning signs are a pooled standard deviation noticeably below one in Gaussian space (variance loss from a truncated search), an exceedance map that looks speckled rather than patchy (too few realisations), and an E-type curve that flattens instead of continuing to fall.

Critical Best Practices #

Simulate with simple kriging and a zero mean #

The transform has already fixed the mean of the variable at zero and its variance at one, so the mean is known, not estimated. Ordinary kriging re-estimates it from each search neighbourhood and adds that estimation variance to every local conditional distribution, inflating the simulated histogram beyond the data’s. Pass gs.krige.Simple(model, cond_pos, cond_val, mean=0.0) and check the pooled standard deviation afterwards. If the variable has a genuine trend, detrend before the transform and simulate the residual — the surfaces produced by ordinary and universal kriging serve as the trend component to add back.

Fit the variogram on the normal scores #

The variogram governs the simulation, and the simulation is Gaussian, so the model must describe the transformed variable. A range fitted to raw lognormal concentrations is typically shorter and the nugget proportionally larger, because a few extreme values dominate the squared differences at every lag. Refit after the transform, and treat model.var + model.nugget ≈ 1 as a pass/fail assertion rather than a curiosity.

Decide the tails of the back-transform explicitly #

np.interp clamps anything outside the lookup table to the endpoints. A realisation containing a normal score of 3.4 when the highest sample scored 2.75 will be assigned the maximum observed concentration exactly, and if 500 realisations do that a few thousand times you get a spike in the pooled histogram and an exceedance area biased low. Fit a tail model — a linear extrapolation in score space, or a hyperbolic upper tail — and state the maximum you are willing to simulate. Sites where the answer turns on extreme values are exactly where the default clamping does most damage.

Check the E-type in Gaussian space, never after back-transforming #

The convergence of the ensemble mean on the kriging estimate holds in the space where the kriging was done. Back-transformation is a non-linear map, so the mean of back-transformed realisations is not the back-transform of the mean — that is precisely the 186 versus 213 gap in the table above. Compare sims_ns.mean(axis=0) with k_ns. Comparing etype with k_map and finding a discrepancy will send you hunting for a bug that does not exist.

Never reuse one Krige object for both the reference map and CondSRF #

gs.CondSRF holds a reference to the estimator you pass it and calls set_pos on it. If you evaluate the same object separately for your kriged comparison map, you will get positions and results that depend on call order. Construct two gs.krige.Simple instances from the same model and data, as in step 5. It costs nothing and removes an entire class of silent inconsistency.

Troubleshooting #

Next Steps #

With an ensemble in hand, compare its cell-wise spread against the pointwise numbers from mapping kriging variance surfaces in Python — they should agree marginally and disagree entirely about anything joint. Then take the aggregate distribution back to whoever set the 5 ha threshold, because 0.846 is a number a contract can be written against and 6.61 ha is not.

Frequently Asked Questions #

Why do the realisations look noisier than the kriged map?

Because they are meant to. Kriging minimises error variance, and the price of that is a surface whose variance is smaller than the data’s, increasingly so away from samples. A realisation is drawn from the conditional distribution rather than summarised by it, so it carries the full model variance everywhere and reproduces the variogram. The kriged map is the better single prediction at any one point; the ensemble is the only thing that gives an honest answer about the joint behaviour of many points at once.

How many realisations do I actually need?

Choose from the precision the answer needs, not from habit. A probability estimated from $n$ realisations has standard error at most $0.5/\sqrt{n}$, so 100 realisations give roughly 0.05, 500 give 0.022 and 1000 give 0.016. If a decision turns on whether an exceedance probability sits above 0.90, 100 realisations are not enough to resolve it. Quantiles deep in the tail of an aggregate distribution need more than the median does.

Should the E-type mean equal the kriged surface?

In Gaussian space, yes: the average of conditional realisations converges on the simple kriging estimate at a rate proportional to one over the square root of the realisation count, and that is the cheapest correctness check available. After back-transformation the equality breaks. Back-transforming the kriged normal score gives the conditional median in original units, whereas averaging the back-transformed realisations gives the conditional mean. For a skewed variable these differ systematically, and the E-type mean is the one that matches the data mean.

Can I use ordinary kriging inside the simulation loop?

Not without consequences. Sequential Gaussian simulation assumes a stationary standard normal variable whose mean is known to be zero, which is exactly the assumption simple kriging encodes. Ordinary kriging re-estimates the mean from the search neighbourhood, adding variance to the local conditional distribution and inflating the spread of the simulated histogram beyond the data’s. Use gs.krige.Simple with mean=0.0. If the variable genuinely has a trend, remove it before the transform and simulate the residual.

Related

• Mapping Kriging Variance Surfaces in Python — the pointwise variance that simulation extends to joint questions
• Building Prediction Intervals from Kriging Variance — the per-location interval, and where it stops being enough
• Ordinary & Universal Kriging — the estimators that supply the trend component and the reference map

← Back to Uncertainty & Variance Mapping