Geostatistics

The problem, stated plainly

You have 155 soil samples from a 4 km × 4 km block in Kiambu. You need a map of soil organic carbon for the whole block. Classical statistics — the statistics you learned in second year — will give you a mean and a standard deviation, and both are useless for the question, because they assume your 155 observations are independent draws from a population. They are not. Two samples 20 m apart are telling you nearly the same thing. Two samples 3 km apart are telling you different things. That structure is not a nuisance to be corrected for; it is the entire signal you need in order to predict at the 3,999,845 locations you did not sample.

This is the founding observation of geostatistics, made by Georges Matheron in the early 1960s while working on the estimation of ore grades in French mines, building on the empirical work of Danie Krige in the South African gold fields. Their problem was commercial and unforgiving: mining a block of ore that turns out to be waste costs money, and so does leaving payable ore in the ground. They needed estimates with a stated error, not just estimates.

The random function, and the one-realisation problem

Formally we treat the value at every location x as a random variable Z(x). The collection of all of them across the region is a random function. Our data — the 155 measurements — are one realisation of that random function, denoted z(x₁), …, z(x₁₅₅).

There is an immediate and serious problem with this. In ordinary statistics you estimate the properties of a random variable by taking many independent draws from it. Here we have one draw from Z(x₁), one from Z(x₂), and so on: one realisation of the whole field. There is only one Kiambu, and we sampled it once. How can you estimate a variance from a single observation?

The answer is to trade an assumption for the data we lack. If we are willing to assume that the statistical properties are the same everywhere in the region, then observations at different locations can substitute for repeated observations at the same location. That assumption is stationarity, and it is the price of admission to the entire subject.

Three grades of stationarity

Stationarity comes in strengths. We want the weakest one that still lets us do the work, because every assumption we make is a claim about reality that could be wrong.

Assumption	What it claims	Verdict
Strict	The entire joint distribution is invariant under translation. Everything about the statistics is the same everywhere.	Far too strong; unverifiable
Second-order	The mean is constant, and the covariance between two points depends only on the vector separating them, not on where they are.	Usable, but requires a finite variance
Intrinsic	The mean of the difference between two points is zero, and the variance of that difference depends only on the separation.	Weakest; what we actually assume

Intrinsic stationarity is weaker in a way that matters. It says nothing about the variance of Z itself — which may not even be finite — only about the variance of differences. Many real phenomena drift upward indefinitely across a region and have no meaningful global variance, but still have well-behaved local differences. This is precisely why the variogram, which is built from differences, became the central tool of the subject rather than the covariance function.

Stated formally, the intrinsic hypothesis is two conditions:

Support, and why it changes everything

Support is the size, shape and orientation of the volume on which a measurement is made. A soil sample taken with a 5 cm auger has a different support from a bulked sample composited across a 10 m plot, which has a different support again from a 30 m Landsat pixel.

Support is not a technicality. As support increases, the variance of the measurements decreases — averaging over a larger volume smooths out small-scale variation. A dataset of 5 cm cores and a dataset of 10 m composites from the same field will have genuinely different histograms, different variances, and different variograms, and neither is wrong. This is why mining engineers estimate block grades rather than point grades: the mine extracts blocks, so the block is the relevant support, and the variance of block grades is much lower than the variance of drill-core grades. We return to this in week 6 under change of support.

What sampling design does to what follows

You will usually inherit a dataset rather than design the survey, but you must be able to read what the design has done to your data. Four patterns recur:

• Simple random. Unbiased, but leaves large gaps and clusters by chance, and yields few short-lag pairs — which is exactly where the variogram is most informative.
• Systematic (grid). Even coverage and easy to execute, but risks aliasing if the phenomenon has periodicity matching the grid spacing, and again supplies no short-lag pairs.
• Stratified random. One random point per grid cell. Usually the best general-purpose compromise.
• Nested or hierarchical. Clusters of points at deliberately varied separations — 1 m, 10 m, 100 m, 1 km. Ugly on a map, superb for variogram estimation, because it populates every lag class including the short ones.

The last point is worth holding on to. A survey designed to estimate a mean and a survey designed to estimate a variogram are different surveys. If you are ever asked to design one, ask which it is.

import geopandas as gpd
import pandas as pd
import matplotlib.pyplot as plt

# CSV with columns: lon, lat, soc  (soil organic carbon, g/kg)
df = pd.read_csv("data/kiambu_soc.csv")

gdf = gpd.GeoDataFrame(
    df,
    geometry=gpd.points_from_xy(df.lon, df.lat),
    crs="EPSG:4326",          # degrees — NOT a distance unit
)
gdf = gdf.to_crs("EPSG:32737")   # UTM 37S — metres

gdf["X"] = gdf.geometry.x
gdf["Y"] = gdf.geometry.y

ax = gdf.plot(column="soc", markersize=gdf.soc * 3,
              legend=True, cmap="viridis", figsize=(7, 7))
ax.set_title("Soil organic carbon, Kiambu (g/kg)")
ax.set_axis_off()
plt.show()

Reading. Webster & Oliver, Geostatistics for Environmental Scientists, ch. 2–3. Isaaks & Srivastava, ch. 1. For the historical account, Cressie (1990), "The origins of kriging", Mathematical Geology 22(3).

Do not skip this week

There is a strong temptation, once you know how to fit a variogram, to load the data and fit one. Resist it. Almost every failed geostatistical analysis fails here, in the hour before any modelling, because someone did not look at their data. The variogram is a summary statistic and, like every summary statistic, it can be dominated by a handful of values that should never have entered the calculation.

Distributions and transforms

Kriging is a linear estimator and, strictly, it does not require normality. But two things make the distribution matter anyway. First, the variogram estimator is a squared difference, so it is highly sensitive to extreme values — a single erroneous outlier can inflate the entire variogram. Second, the simulation methods in week 8, and the confidence intervals you will want to attach to your predictions, do assume a Gaussian distribution.

So: plot the histogram, plot the Q–Q plot, compute the skewness. As a working rule, a coefficient of skewness beyond about ±1 is worth acting on. Two transforms are standard:

• Logarithmic. Appropriate for strictly positive, right-skewed variables — most concentrations, most trace metals, rainfall. Simple, interpretable, but back-transforming a kriged estimate is not as simple as exponentiating: that gives you the median, not the mean, and the bias correction is a real subject (lognormal kriging).
• Normal score. Rank the data, map each rank onto the corresponding quantile of a standard normal. Guarantees an exactly normal marginal distribution, which is what sequential Gaussian simulation requires. Back-transformation is by the same lookup table, run in reverse.

Outliers, and the spatial version of the question

Classical outlier detection asks whether a value is extreme relative to the distribution. Spatial data lets you ask a better question: is a value extreme relative to its neighbours? A cadmium concentration of 18 mg/kg might be perfectly ordinary in the distribution and yet be a clear error if every sample within 200 m reads below 2.

The practical tool is the local deviation: for each point, compare its value to the mean of its k nearest neighbours, and standardise. Points with a large standardised deviation are candidates for investigation. Note the word: investigation, not deletion. Delete a value only when you have a reason external to the statistics — a transcription error, a known contaminated sample, a GPS fix that puts the point in a lake. An extreme value that is real is often the most important observation in the dataset, and a hotspot you deleted for tidiness is a hotspot you failed to map.

Trend, and the decision it forces

A trend (or drift) is a systematic change in the mean across the region — soil depth increasing downslope, temperature decreasing with altitude, contamination decreasing with distance from a smelter. Trend violates the constant-mean assumption, and if you ignore it the variogram will not reach a sill: it climbs indefinitely, because points far apart differ systematically rather than randomly.

Detect it by regressing the variable on the coordinates, and by looking at the shape of the empirical variogram at long lags. Then choose:

1. Restrict the region. If the trend is a large-scale feature and your area of interest is small, work locally where the mean is plausibly constant. The cheapest solution when it is available.
2. Model the trend and krige the residuals. Fit the trend (as a function of coordinates, or better of a covariate such as elevation), work with what is left over, then add the trend back. This is universal kriging and regression kriging — week 6.
3. Ignore it. Occasionally defensible if the trend is weak relative to the local variation and your predictions are all short-range interpolations. Say so explicitly if you do.

Preferential sampling and declustering

Real surveys are rarely evenly spread. Geologists drill densely where the ore is; environmental agencies sample intensively where they suspect contamination; field teams sample near roads. The result is that the naive sample mean is biased toward whatever the clustered areas contain — usually the high values, because that is why people clustered there.

Cell declustering is the standard correction. Lay a grid over the region; each occupied cell receives a total weight of one, shared equally among the points inside it. A point alone in its cell gets weight 1; each of ten points crowded into one cell gets weight 0.1. The declustered mean is then the weighted mean.

import numpy as np
from scipy import stats
from sklearn.neighbors import NearestNeighbors

def local_deviation(coords, values, k=8):
    """Standardised difference between each value and its k nearest neighbours."""
    nn = NearestNeighbors(n_neighbors=k + 1).fit(coords)
    _, idx = nn.kneighbors(coords)
    nbr_mean = values[idx[:, 1:]].mean(axis=1)   # column 0 is the point itself
    resid = values - nbr_mean
    return resid / resid.std()

coords = gdf[["X", "Y"]].values
gdf["ld"] = local_deviation(coords, gdf.soc.values)

print("skewness:", stats.skew(gdf.soc))
print(gdf.reindex(gdf.ld.abs().sort_values(ascending=False).index).head(5)
         [["X", "Y", "soc", "ld"]])

# --- cell declustering -------------------------------------------------
def declustered_mean(coords, values, cell):
    keys = np.floor(coords / cell).astype(int)
    _, inv = np.unique(keys, axis=0, return_inverse=True)
    counts = np.bincount(inv)
    w = 1.0 / counts[inv]           # each occupied cell contributes weight 1
    return np.average(values, weights=w)

for cell in [100, 200, 400, 800, 1600]:
    print(cell, round(declustered_mean(coords, gdf.soc.values, cell), 3))
print("naive  ", round(gdf.soc.mean(), 3))

Reading. Isaaks & Srivastava, ch. 2–4 — still the best treatment of exploratory work in the subject. Deutsch & Journel, GSLIB, on declustering. Diggle, Menezes & Su (2010), "Geostatistical inference under preferential sampling", JRSS-C 59(2), for the serious version of the problem.

From "nearby things are similar" to a number

Last week's claim was that neighbouring values are related. This week we measure it. Take every pair of points separated by roughly the same distance h, square the difference in their values, and average. If nearby things really are similar, that average will be small for small h and grow as h increases.

That is the whole idea. The estimator, due to Matheron, is:

Under second-order stationarity the variogram and the covariance function carry the same information:

Binning: the three choices you must make

In two dimensions, no two pairs have exactly the same separation, so we bin. Three decisions follow, and all three change the answer.

• Lag width. Narrow bins resolve fine structure but each holds few pairs and the result is noisy. Wide bins are stable but smear the short-lag behaviour, which is precisely the part that matters most for kriging. Start with the mean nearest-neighbour distance as your first lag width.
• Maximum lag. Half the diagonal of the sampled region, at most. Beyond that the estimate is built from a shrinking number of pairs drawn from the edges of the survey and is not to be trusted.
• Lag tolerance. Usually half the lag width, so bins tile without gaps. For directional variograms (week 4) you also need an angular tolerance — typically ±22.5°, which divides the plane into four non-overlapping sectors.

Two diagnostics worth the time

The variogram cloud plots every individual pair — separation on the x-axis, half the squared difference on the y-axis — before any averaging. Points sitting far above the general mass at short lags mark pairs of nearby locations with wildly different values, and tracing them back to specific samples is often how a data error is caught.

The h-scatterplot plots z(x) against z(x+h) for one lag class. At short lags the cloud should hug the 1:1 line; as the lag grows it should fatten toward a circular blob. Watching that sequence is the most intuitive picture of spatial correlation decaying that the subject offers, and it makes the variogram feel less like an incantation.

import numpy as np
import skgstat as skg
import matplotlib.pyplot as plt

coords = gdf[["X", "Y"]].values
vals   = np.log(gdf.zinc.values)          # week 2 said transform; we did

# half the diagonal of the bounding box — the sensible maximum lag
xmin, ymin, xmax, ymax = gdf.total_bounds
maxlag = 0.5 * np.hypot(xmax - xmin, ymax - ymin)

V = skg.Variogram(coords, vals,
                  n_lags=15,
                  maxlag=maxlag,
                  estimator="matheron",
                  normalize=False)

print(V)                                   # bins, pair counts, fitted parameters
print("pairs per bin:", V.bin_count)      # check none is below ~30

# sensitivity to binning — never report just one
fig, ax = plt.subplots(figsize=(7, 4.5))
for n in (8, 15, 30):
    Vn = skg.Variogram(coords, vals, n_lags=n, maxlag=maxlag)
    ax.plot(Vn.bins, Vn.experimental, "o-", label=f"{n} bins", alpha=0.8)
ax.set_xlabel("lag h (m)"); ax.set_ylabel("semivariance γ(h)"); ax.legend()
plt.show()

Reading. Webster & Oliver, ch. 4. Isaaks & Srivastava, ch. 7. Cressie & Hawkins (1980) on the robust estimator, Mathematical Geology 12(2). The scikit-gstat documentation is unusually good and should be read as a text, not consulted as a reference.

Why we cannot use the empirical variogram directly

Kriging needs a semivariance value for every separation that arises in the calculation — including separations that fall between your bins, and separations shorter than anything you sampled. The empirical variogram gives you fifteen points. So we fit a continuous function through them, and from that point on the model, not the data, is what kriging uses.

There is a second and less obvious reason. The kriging system requires that the matrix of semivariances be invertible for any configuration of points. Only functions with a particular mathematical property — conditional negative-definiteness — guarantee this. Fit an arbitrary curve, and you can obtain negative kriging variances, which are meaningless. This is why we choose from a small catalogue of permissible models rather than fitting whatever passes closest through the points.

The anatomy of a variogram

Each parameter has a physical reading, and you should always give one.

• Nugget — named for the gold nuggets whose presence or absence made two adjacent assays differ wildly. It has two sources you cannot separate from the variogram alone: genuine variation at scales below your minimum sample spacing, and measurement error. A nugget that is 60% of the sill is a warning: most of your variance is unstructured, and kriging will not do much for you. If you need to know which source dominates, take duplicate samples at the same location — the duplicate variance is the measurement error.
• Sill — under second-order stationarity this equals the variance of the data, which gives you a useful sanity check: fit a model whose sill is wildly different from the sample variance and something is wrong.
• Range — the distance over which the phenomenon "remembers" itself. A range of 800 m for soil carbon says that samples more than 800 m apart tell you nothing about each other, which immediately tells you what sampling density a future survey needs.

The permissible models

Nested structures

A sum of permissible models is itself permissible, and real data often needs one. Soil properties frequently show variation at two scales — a short-range structure from local soil formation processes at 50–200 m, and a long-range structure from geology or topography at several kilometres. Fitting a single model to such a variogram fits neither structure. Fit a nested model: c₀ plus a short-range spherical plus a long-range spherical. Each component then has a physical interpretation you can name, which is a useful discipline — if you cannot name the process behind a component, ask whether you are fitting noise.

Anisotropy

So far we have assumed the variogram depends only on the length of the lag vector, not its direction. That is often false. Soil properties along a valley behave differently from soil properties across it; a plume of contamination is elongated in the direction of groundwater flow; rainfall correlates further along a mountain range than across it.

Detect it by computing directional variograms: restrict pairs to those whose separation vector falls within an angular tolerance of a chosen direction. Compute at 0°, 45°, 90° and 135° with ±22.5° tolerance, which covers the plane without overlap. Then read the result:

• Geometric anisotropy — same sill, different ranges by direction. Common, and easy to handle: the tool rescales the coordinate system by the anisotropy ratio and direction, turning ellipses back into circles, and everything downstream proceeds as before.
• Zonal anisotropy — different sills by direction. Usually a sign of a trend in one direction, and better handled by removing that trend (week 6) than by modelling the anisotropy.

Report anisotropy as a ratio and an angle: "geometric anisotropy, ratio 2.4, major axis 35°" tells a reader everything they need. And beware the trap — with four directions and typical sample sizes, each directional variogram has a quarter of the pairs of the omnidirectional one, so apparent anisotropy in a small dataset is often noise. Check whether it survives a change of binning before you believe it.

Fitting

Two approaches, and you should know both:

• Weighted least squares on the empirical variogram points, weighting each bin by its pair count (and often inversely by γ²) so that well-populated short lags dominate the fit. This is what skgstat does by default. Fast, robust, and standard in practice.
• Maximum likelihood or REML on the raw data, skipping the empirical variogram entirely. Statistically more efficient and provides standard errors on the parameters; more expensive and more sensitive to the assumption of normality. Standard in the statistical literature, less so in applied geoscience.

Whichever you use, do not treat fitting as automatic. Plot the fit over the empirical points and look at it. A fit that is excellent at long lags and poor at short ones is a bad fit for kriging purposes, no matter what the objective function says, because short lags are what determine the weights.

# --- compare models on the same empirical variogram --------------------
V = skg.Variogram(coords, vals, n_lags=15, maxlag=maxlag)

for m in ["spherical", "exponential", "gaussian", "matern"]:
    V.model = m
    c0 = V.parameters[-1]                 # nugget is reported last
    print(f"{m:12s} range={V.parameters[0]:8.1f}  "
          f"sill={V.parameters[1]:6.3f}  nugget={c0:6.3f}  rmse={V.rmse:6.4f}")

# --- directional variograms -------------------------------------------
fig, ax = plt.subplots(figsize=(7, 4.5))
for az in (0, 45, 90, 135):
    DV = skg.DirectionalVariogram(coords, vals,
                                  azimuth=az,
                                  tolerance=22.5,     # four non-overlapping sectors
                                  n_lags=12,
                                  maxlag=maxlag)
    ax.plot(DV.bins, DV.experimental, "o-", label=f"{az}°", alpha=0.85)
ax.set_xlabel("lag h (m)"); ax.set_ylabel("γ(h)")
ax.legend(title="azimuth"); plt.show()

Reading. Webster & Oliver, ch. 5–6. Chilès & Delfiner, Geostatistics: Modeling Spatial Uncertainty, ch. 2, for the mathematics of permissible models. Stein, Interpolation of Spatial Data (1999), ch. 2, on why the Matérn deserves its status.

What we are asking for

We want an estimate at an unsampled location x₀, formed as a weighted average of the samples:

Every interpolation method you have met is of this form. Nearest neighbour sets one weight to 1 and the rest to 0. Inverse distance weighting sets λ_i ∝ 1/d_i^p. What distinguishes kriging is that it does not assert the weights; it derives them, from two requirements.

Unbiasedness. On average the estimate should equal the truth, E[ẑ(x₀) − Z(x₀)] = 0. With an unknown but constant mean, this forces the weights to sum to one.

Minimum variance. Among all weight sets satisfying that constraint, choose the one minimising the estimation variance Var[ẑ(x₀) − Z(x₀)].

This is a constrained optimisation, handled with a Lagrange multiplier ψ. Differentiating and setting to zero gives a linear system of n+1 equations:

Having solved for the weights, the minimised variance itself is available:

Look carefully at the left-hand matrix in (5.2). It contains γ between samples — how redundant the samples are with each other — and the right-hand side contains γ between samples and the target — how relevant each sample is. Kriging balances relevance against redundancy. That is the entire idea, and the worked example below makes it concrete.

What drives the weights

• Distance to the target. The obvious one: nearer samples get more weight, through the right-hand side of the system.
• Clustering among the samples. As above. Redundant samples share the weight that one of them would have received alone.
• Screening. A sample directly behind another, along the same line from the target, receives very little weight — the nearer one shields it. In extreme configurations screened samples receive small negative weights. These are legitimate and are what allows kriging to extrapolate slightly beyond the data range, but they can also produce impossible values such as a negative concentration. If that matters for your variable, either constrain the solution or note the occurrences.
• The nugget. As the nugget grows relative to the sill, the weights flatten toward equal. In the limit of a pure nugget model, kriging returns the mean of the neighbourhood at every location — there is no structure left to exploit, and it tells you so.

The kriging variance, and the trap in it

Equation (5.3) yields a variance at every location, and mapping it alongside the prediction is one of the most useful things geostatistics gives you. But read (5.3) again and notice what is not in it: the data values. The kriging variance is built entirely from the variogram model and the geometry of the sampling configuration. Move all your sample values around at random, keeping the locations fixed, and the kriging variance map does not change by a single unit.

Search neighbourhoods

Using all n samples for every prediction means inverting an n×n matrix at every grid node, which is impractical beyond a few thousand points and unnecessary anyway: samples beyond the range contribute essentially nothing. So we use a moving neighbourhood — typically the nearest 16 to 32 samples, sometimes with quadrant or octant constraints forcing the search to draw from all directions rather than taking 20 points from one dense cluster.

The choice is a real one. Too few samples gives noisy, discontinuous maps with visible artefacts at the boundaries where the neighbourhood changes composition. Too many is slow for no benefit. Sixteen to twenty-four with quadrant balancing is a reasonable default; report what you used.

import numpy as np

def spherical(h, c0=2.0, c1=8.0, a=250.0):
    h = np.asarray(h, dtype=float)
    g = np.where(h <= a,
                 c0 + c1 * (1.5 * h / a - 0.5 * (h / a) ** 3),
                 c0 + c1)
    return np.where(h == 0, 0.0, g)      # γ(0) = 0 exactly

pts  = np.array([[0., 100.], [30., 95.], [0., -100.]])
z    = np.array([10., 11., 6.])
x0   = np.array([0., 0.])
n    = len(pts)

# Γ: semivariances between samples, bordered with ones for the constraint
D = np.linalg.norm(pts[:, None, :] - pts[None, :, :], axis=-1)
G = np.ones((n + 1, n + 1))
G[:n, :n] = spherical(D)
G[n, n]   = 0.0

# γ₀: semivariances between each sample and the target
g0 = np.ones(n + 1)
g0[:n] = spherical(np.linalg.norm(pts - x0, axis=1))

sol    = np.linalg.solve(G, g0)
w, psi = sol[:n], sol[n]

print("weights   ", np.round(w, 4))          # [0.2755 0.2749 0.4496]
print("estimate  ", round(float(w @ z), 3))    # 8.475
print("variance  ", round(float(w @ g0[:n] + psi), 3))
print("naive mean", round(z.mean(), 3))       # 9.0 — what IDW would give

from pykrige.ok import OrdinaryKriging

gridx = np.arange(xmin, xmax, 40.0)
gridy = np.arange(ymin, ymax, 40.0)

OK = OrdinaryKriging(
    coords[:, 0], coords[:, 1], vals,
    variogram_model="spherical",
    variogram_parameters={"sill": 0.62, "range": 900.0, "nugget": 0.05},
    nlags=15, enable_plotting=False,
)

zhat, var = OK.execute("grid", gridx, gridy)
se = np.sqrt(var)

fig, (a1, a2) = plt.subplots(1, 2, figsize=(12, 5))
a1.imshow(zhat, origin="lower", cmap="viridis"); a1.set_title("OK prediction, log-zinc")
a2.imshow(se,   origin="lower", cmap="magma");   a2.set_title("kriging standard error")
a2.scatter(*((coords.T - [[xmin], [ymin]]) / 40.0), s=5, c="white")
plt.show()

Reading. Isaaks & Srivastava, ch. 12 — the clearest derivation in print, worked with actual numbers. Webster & Oliver, ch. 8. Goovaerts, ch. 5.

The assumption we have been making

Ordinary kriging assumes the mean is constant but unknown, at least within the search neighbourhood. Often it is not. Rainfall increases with altitude. Soil depth increases downslope. Contamination decreases away from a source. When the mean varies systematically, ordinary kriging is inefficient at best and biased at worst, and the variogram will not reach a sill.

The family of methods below all decompose the variable into a deterministic part and a stochastic part:

Method	Assumption about m(x)	Use when
Simple kriging	Constant and known	Rarely — but it is the building block for simulation in week 8
Ordinary kriging	Constant, unknown	The default; no covariate available
Universal kriging	A polynomial function of the coordinates	A smooth regional gradient with no physical covariate
KED / external drift	A linear function of a covariate known everywhere	The best case: a covariate that is physically related and mapped densely
Regression kriging	Any regression model, fitted separately	Several covariates, or a non-linear relationship

Kriging with an external drift

This is the method to reach for, because it exploits the situation you are usually in: a sparse set of expensive point measurements, and a cheap covariate available everywhere. The pairings are familiar — rainfall with elevation, temperature with elevation, soil organic carbon with NDVI or with terrain wetness index, groundwater level with topography.

KED does everything at once: it estimates the trend coefficients and the spatial prediction within a single system, using the residual variogram. Two conditions must hold, and both are frequently violated in student work.

• The covariate must be known at every prediction location, not just at the sample points. A DEM qualifies. A second sparse point dataset does not — that is a co-kriging problem.
• The relationship should be linear and reasonably strong. If the scatterplot of z against the covariate is a shapeless cloud, KED will add complexity and nothing else.

The critical technical point: the variogram used in KED is the variogram of the residuals, not of the raw variable. Fit the trend, take the residuals, compute and model their variogram, and use that. The residual variogram will typically have a shorter range and a lower sill than the raw one — that difference is precisely the structure the covariate has explained.

Regression kriging

Regression kriging is the same idea run as three explicit steps, which makes it easier to understand, easier to debug, and more flexible:

1. Fit a regression of z on the covariates at the sample locations, by any method — OLS, GLM, or in week 10, random forest.
2. Compute the residuals, model their variogram, and krige the residuals across the region.
3. Add the kriged residual surface to the regression prediction surface.

With OLS in step 1 and consistent variogram estimation, regression kriging and KED are theoretically equivalent. In practice RK is preferred whenever you want several covariates, a non-linear model, or simply to see the two components separately — and being able to map the trend surface and the residual surface individually is genuinely informative about where your covariate is working.

Change of support: block kriging

Everything so far has predicted a value at a point. Often the decision is about an area: the mean grade of a 10 m mining block, the mean soil carbon of a one-hectare field, the mean rainfall over a catchment. Block kriging estimates the average over a block V directly, by replacing the point-to-point semivariances on the right-hand side of the kriging system with point-to-block average semivariances.

The consequence to remember: block estimates have substantially lower kriging variance than point estimates, because averaging over an area removes short-scale variability. This is not a trick — the block mean genuinely is more predictable than any individual point within it. It also means that comparing a point-kriging standard error with a block-kriging standard error is meaningless unless you say which is which.

Co-kriging, briefly

When the secondary variable is also sparsely sampled — not available everywhere — KED is unavailable and co-kriging is the formally correct method. It requires modelling the variogram of each variable plus the cross-variogram between them, all under the linear model of coregionalisation, which constrains the models to share ranges and produces a positive-definite system.

It is a lot of machinery. Be honest about when it pays: co-kriging helps materially only when the secondary variable is markedly more densely sampled than the primary and strongly correlated with it. If both are sampled at the same locations, the gain is usually negligible and not worth the modelling burden. Know that it exists, know the condition under which it is worth doing, and reach for KED or RK first.

from sklearn.linear_model import LinearRegression
from pykrige.uk import UniversalKriging

# ---- regression kriging, done as three explicit steps -----------------
X = stations[["elev"]].values
y = stations["rain"].values

lm = LinearRegression().fit(X, y)
print(f"R² = {lm.score(X, y):.3f}   slope = {lm.coef_[0]*100:.1f} mm per 100 m")

resid = y - lm.predict(X)

# the variogram that matters is the RESIDUAL variogram, not the raw one
Vr = skg.Variogram(stations[["X", "Y"]].values, resid,
                   n_lags=12, model="exponential")
print("residual variogram:", Vr.parameters)

OKr = OrdinaryKriging(stations.X, stations.Y, resid,
                      variogram_model="exponential",
                      variogram_parameters={"sill": Vr.parameters[1],
                                            "range": Vr.parameters[0],
                                            "nugget": Vr.parameters[2]})
resid_grid, _ = OKr.execute("grid", gridx, gridy)

trend_grid = lm.predict(dem_grid.reshape(-1, 1)).reshape(dem_grid.shape)
rk_grid    = trend_grid + resid_grid      # the two components, added back

# ---- the same thing in one step, as KED -------------------------------
UK = UniversalKriging(stations.X, stations.Y, y,
                      variogram_model="exponential",
                      drift_terms=["specified"],
                      specified_drift=[stations.elev.values])
ked_grid, ked_var = UK.execute("grid", gridx, gridy,
                               specified_drift_arrays=[dem_grid])

Reading. Hengl, Heuvelink & Stein (2004), "A generic framework for spatial prediction of soil variables based on regression-kriging", Geoderma 120 — the standard reference for RK. Goovaerts, ch. 6. Webster & Oliver, ch. 9.

The question this week answers

You have a map. Someone is going to spend money on the basis of it. How good is it, and how do you know?

The kriging variance will not answer this — week 5 established that it never saw your data values. The only honest route is to withhold observations, predict them, and compare. That is cross-validation, and it is the closest thing this subject has to a truth serum.

The diagnostics, and what each one catches

Leave-one-out cross-validation removes each observation in turn, predicts it from the remainder, and collects the errors e_i = z(x_i) − ẑ(x_i). Four statistics follow, and you should report all four, because each catches a different failure.

MSDR is the only one of the four that tests the variance model rather than the prediction. A model can have an excellent RMSE and an MSDR of 3, which means the map is reasonably accurate but the error bars you would quote from it are wrong by a factor of nearly two. If you are going to attach uncertainty to a map — and this unit insists you do — this is the number that says whether that uncertainty is honest.

The smoothing effect

Kriging minimises estimation variance, and one consequence is that the resulting map is smoother than reality: high values are under-estimated, low values over-estimated, and the variance of the kriged surface is lower than the variance of the data. Plotting predicted against observed makes it obvious — the regression line has a slope below one, fanning away from the 1:1 line at both ends.

This is not a bug. It is the correct behaviour for minimising squared error, and it is unavoidable in any optimal linear estimator. But it matters enormously for a particular class of question. If you are asked "what fraction of this district exceeds the WHO fluoride guideline?", counting cells above the threshold on a kriged map will give you a systematically wrong answer, because the map has flattened exactly the extremes you are counting. Kriging is the right tool for "what is the value here"; it is the wrong tool for "how much of the area exceeds a threshold". That is a simulation question, and it is next week.

Benchmark honestly

Kriging is not automatically the best method. Compare it, on identical folds, against inverse distance weighting, nearest neighbour and a thin-plate spline. There are real situations where the simpler methods win: a very high nugget-to-sill ratio leaves little structure for kriging to exploit; very dense, evenly spaced sampling makes almost any interpolator adequate; small datasets make the variogram too unstable to be worth the trouble.

Report the comparison whichever way it comes out. A study that concludes "IDW performed as well as kriging here, and is simpler, so we used it" is a good study. A study that reports only the kriging results because they were expected is not.

import numpy as np
from sklearn.model_selection import KFold

def loo_cv(coords, values, model, params):
    """Leave-one-out CV returning the four diagnostics that matter."""
    n = len(values)
    pred = np.empty(n); var = np.empty(n)
    for i in range(n):
        keep = np.arange(n) != i
        ok = OrdinaryKriging(coords[keep, 0], coords[keep, 1], values[keep],
                             variogram_model=model, variogram_parameters=params)
        z, s = ok.execute("points", coords[i:i+1, 0], coords[i:i+1, 1])
        pred[i], var[i] = z[0], s[0]

    e = values - pred
    return {
        "ME":   e.mean(),                    # → 0 if unbiased
        "RMSE": np.sqrt((e ** 2).mean()),     # accuracy, in data units
        "MAE":  np.abs(e).mean(),             # less outlier-sensitive
        "MSDR": ((e ** 2) / var).mean(),      # → 1 if the variance model is honest
    }

def spatial_blocks(coords, block_size):
    """Fold labels from a spatial grid — hold out whole blocks, not points."""
    keys = np.floor(coords / block_size).astype(int)
    _, fold = np.unique(keys, axis=0, return_inverse=True)
    return fold

# the comparison that matters: same data, two ways of splitting it
fold_random  = KFold(10, shuffle=True, random_state=0)
fold_spatial = spatial_blocks(coords, block_size=2 * fitted_range)
# ... score both; the gap between them is your optimism

Reading. Webster & Oliver, ch. 8 §8.4. Roberts et al. (2017), "Cross-validation strategies for data with temporal, spatial, hierarchical or phylogenetic structure", Ecography 40(8) — the definitive treatment. Meyer & Pebesma (2022) on area of applicability, Methods in Ecology and Evolution 12(9).

Two different questions

Compare these:

• "What is the fluoride concentration at this proposed borehole site?" — an estimation question. Kriging answers it, and answers it optimally.
• "What fraction of this sub-county has groundwater fluoride above the 1.5 mg/L guideline, and where should we prioritise testing?" — a risk question. Kriging answers it badly.

The reason is last week's smoothing effect. A kriged map has less variance than reality: it has flattened the peaks and filled in the troughs. Counting cells above a threshold on such a map systematically undercounts the area at risk, and the undercount is worst exactly where the threshold is far from the mean — which is usually where the policy question sits.

Simulation gives up optimality at each individual location in exchange for getting the statistical properties of the whole field right. A simulated realisation is not the best estimate anywhere; it is a plausible version of reality that reproduces the correct histogram, the correct variogram, and honours the data at the sample points. Generate many of them and you have a picture of what the field might be, not just a single flattened best guess.

Sequential Gaussian simulation

SGS is the standard algorithm and it is simpler than it sounds. The key insight is that under a multivariate Gaussian model, the conditional distribution at any location — given the data and the previously simulated values — is normal, with mean equal to the simple kriging estimate and variance equal to the simple kriging variance. So you can build a realisation one node at a time:

1. Normal-score transform the data, so the marginal distribution is exactly standard normal. (This is why week 2 mattered.)
2. Model the variogram of the transformed data.
3. Define a random path visiting every grid node exactly once.
4. At each node in turn: simple-krige from the original data and all previously simulated nodes to obtain a mean and a variance; draw one random value from the resulting normal distribution; add that drawn value to the conditioning dataset.
5. When the path is complete, back-transform the whole realisation through the inverse normal-score table.
6. Repeat from step 3 with a different random path and different random draws, for as many realisations as you need — 100 for exploratory work, 500 or more for stable tail probabilities.

Two details are worth pausing on. Step 4's phrase "and all previously simulated nodes" is what makes the realisation reproduce the variogram — each drawn value becomes data for its neighbours, propagating spatial structure across the field. And the random path in step 3 matters: a systematic path would introduce directional artefacts, which is why it must be random and different for each realisation.

Reading a set of realisations

You now have, say, 200 values at every grid node — one from each realisation. That is an empirical distribution, and everything useful follows from it:

• The E-type estimate — the mean across realisations at each node. It converges toward the kriged map as the number of realisations grows, which is a useful check that your simulation is behaving.
• The conditional variance — the spread across realisations at each node. Unlike the kriging variance, this does depend on the data values: it is wider where the local values are variable and narrower where they agree.
• Probability maps — the proportion of realisations exceeding a threshold at each node. This is the deliverable for a risk question, and it is what you cannot obtain from kriging at all.
• Uncertainty in an aggregate — compute any downstream quantity separately on each realisation, and you get a distribution for that quantity. The total tonnage above cut-off; the total area at risk; the number of households affected. This is uncertainty propagation, and it is the most valuable thing in the chapter.

import numpy as np
import gstools as gs        # pip install gstools — the practical SGS route in Python

model = gs.Exponential(dim=2, var=0.85, len_scale=420.0, nugget=0.15)

krige = gs.krige.Ordinary(model, cond_pos=(x_obs, y_obs), cond_val=ns_values)
cond  = gs.CondSRF(krige)

N = 100
sims = np.empty((N, gridy.size, gridx.size))
for i in range(N):
    cond.set_pos((gridx, gridy), "structured")
    sims[i] = cond(seed=i)                       # one realisation

sims = back_transform(sims)                      # inverse normal score

etype  = sims.mean(axis=0)                      # → converges toward the kriged map
cstd   = sims.std(axis=0)                       # depends on DATA, unlike kriging variance
p_exc  = (sims > 1.5).mean(axis=0)             # P(exceeding the guideline)

# uncertainty in an aggregate — the number a decision-maker can use
cell_km2  = 0.04
area_risk = (sims > 1.5).sum(axis=(1, 2)) * cell_km2
print(f"area above guideline: {np.median(area_risk):.1f} km² "
      f"(90% interval {np.percentile(area_risk, 5):.1f}–"
      f"{np.percentile(area_risk, 95):.1f})")

Reading. Goovaerts, ch. 8 — the standard treatment. Deutsch & Journel, GSLIB, on sgsim. Chilès & Delfiner, ch. 7. The gstools documentation for the Python implementation.

When the Gaussian assumption will not hold

Simulation last week required a normal-score transform, which requires a continuous variable with a well-defined distribution. Plenty of real problems do not fit: heavily skewed data with a long tail no transform tames; data with a spike of values at a detection limit; and variables that are not numeric at all — soil class, land cover, rock type, presence or absence.

The indicator approach handles all of these by giving up on estimating the value and estimating a probability instead.

The indicator transform

Pick a threshold z_k and replace every observation with a 0 or a 1:

This is a remarkable trick, because the expected value of the indicator is the probability of being below the threshold. So if you krige the indicator variable, the estimate you get back — a number between 0 and 1 — is a direct estimate of P(Z(x₀) ≤ z_k), with no distributional assumption at all.

Two further benefits come free. Extreme values lose their leverage: a wildly high outlier becomes a 0, exactly like a moderately high value, so the indicator variogram is far more robust than the raw one. And you can use a different variogram at each threshold, which lets the spatial structure of the high values differ from that of the low ones — often physically true, and something ordinary kriging cannot express.

Building a full local distribution

Repeat the procedure at a series of thresholds — typically the nine deciles of the sample distribution — and krige each indicator separately. At every prediction location you now have nine estimated probabilities, which together trace out the local cumulative distribution function. From that you can read the median, any quantile, an interquartile range as an uncertainty measure, and the probability of exceeding any threshold you care about.

Categorical variables

For a categorical variable with K classes, define one indicator per class and krige each. At every location you obtain K probabilities, which after correction sum to one. Two products follow: assign each location its most probable class for a map, and map the probability itself for a picture of where the classification is confident and where it is not. That second map is usually the more useful of the two, and it is the sort of thing that would improve a great many published soil and land-cover maps.

Sequential indicator simulation follows the same logic as SGS but draws from the indicator-derived local CDF instead of a Gaussian one — which is how you produce realisations of a categorical field, such as alternative plausible soil maps, each internally consistent.

From probability to decision

A probability map is not yet a decision. Turning it into one requires a threshold on the probability, and choosing that threshold is a value judgement about the relative cost of the two errors, not a statistical one.

Consider a fluoride exceedance map used to prioritise borehole testing. Flagging every location with P > 0.1 casts a wide net: few contaminated boreholes missed, many clean ones tested unnecessarily. Flagging at P > 0.7 tests far fewer, but misses more contaminated sites. The right answer depends on the cost of a test against the cost of a child drinking fluoride-contaminated water for a year, and that is a decision for the public health authority, not the analyst.

Your job is to make the trade-off visible: present the probability map, state the consequences of several candidate thresholds, and let the decision-maker choose with their eyes open. This is what it means to communicate uncertainty rather than merely compute it, and it is a skill this unit is assessing.

# --- indicator kriging at a regulatory threshold -----------------------
threshold = 1.5                                 # mg/L, WHO guideline
ind = (gdf.fluoride.values > threshold).astype(float)
print(f"proportion above guideline: {ind.mean():.3f}")

Vi = skg.Variogram(coords, ind, n_lags=12, model="spherical")

OKi = OrdinaryKriging(coords[:, 0], coords[:, 1], ind,
                      variogram_model="spherical",
                      variogram_parameters={"sill": Vi.parameters[1],
                                            "range": Vi.parameters[0],
                                            "nugget": Vi.parameters[2]})
prob, _ = OKi.execute("grid", gridx, gridy)
prob = np.clip(prob, 0, 1)                       # kriging does not respect [0,1]

# --- make the decision trade-off explicit ------------------------------
for p in (0.1, 0.3, 0.7):
    flagged = (prob > p).sum()
    print(f"P > {p:.1f}: flag {flagged:5d} cells "
          f"({100*flagged/prob.size:.1f}% of the area)")

Reading. Goovaerts, ch. 7. Journel (1983), "Nonparametric estimation of spatial distributions", Mathematical Geology 15(3). Deutsch & Journel, GSLIB, on ik3d and sisim.

The honest framing

Most operational mapping in soils, ecology and Earth observation is now done with machine learning rather than kriging. SoilGrids, which produces global soil property maps at 250 m, is built on gradient-boosted trees and random forests, not on variograms. It would be dishonest to teach this unit as though that had not happened.

It would be equally dishonest to present it as a replacement. The two traditions answer overlapping but different questions, and the sharpest way to see the difference is this: a random forest gives you a prediction; kriging gives you a prediction and a model of the error. Where a stated uncertainty is required — a resource statement, a regulatory exceedance assessment, anything that will be challenged — that difference decides the matter.

Why a random forest on coordinates fails

Give a random forest X and Y as features and it will fit the training data well and generalise badly. Trees partition space with axis-parallel splits, producing blocky, rectangular artefacts; the model has no notion of distance and cannot interpolate smoothly; and outside the convex hull of the training data it predicts a constant, because trees cannot extrapolate at all.

Two established fixes:

• Regression kriging with an ML trend (week 6, generalised). Fit the random forest on genuine covariates — elevation, slope, NDVI, rainfall, distance to river — then model and krige the residuals and add them back. The forest captures non-linear covariate relationships; kriging captures the spatial structure the covariates missed. This is usually the best-performing and most defensible option.
• RFsp — distance features. Instead of coordinates, give the forest the distance from each prediction location to every training point as features. This gives the trees the spatial information they lack in a form they can split on. It works, though it produces a very wide feature matrix and scales poorly.

What each side gives up

	Kriging	Random forest
Uncertainty	Explicit variance from the model; simulation gives full distributions	None natively; quantile forests and conformal methods approximate it
Covariates	Few, and linear unless you work at it	Many, non-linear, with interactions, essentially free
Extrapolation	Degrades gracefully toward the mean with growing variance	Fails silently — predicts a constant with no warning
Interpretability	Three parameters, each physically meaningful	Variable importance and partial dependence, both easy to over-read
Small samples	Workable from about 50–100 points	Needs several hundred at minimum

The practical conclusion is unglamorous and correct: use both, compare them under a fair protocol, and report what you find. A candidate who can say "random forest with kriged residuals beat plain kriging by 8% RMSE under spatial block cross-validation, but I used kriging for the exceedance assessment because the client needed confidence intervals" is demonstrating exactly the judgement this unit exists to build.

from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import GroupKFold

# Spatial blocks become GROUPS — GroupKFold then never splits a block.
groups = spatial_blocks(coords, block_size=2 * fitted_range)
gkf    = GroupKFold(n_splits=5)

results = []
for train, test in gkf.split(covars, y, groups=groups):
    rf = RandomForestRegressor(n_estimators=500, min_samples_leaf=3,
                               random_state=0, n_jobs=-1)
    rf.fit(covars[train], y[train])                # fit INSIDE the fold

    # --- RF + regression kriging of the residuals ----------------------
    res = y[train] - rf.predict(covars[train])
    ok  = OrdinaryKriging(coords[train, 0], coords[train, 1], res,
                          variogram_model="exponential")
    res_hat, _ = ok.execute("points", coords[test, 0], coords[test, 1])

    pred_rf   = rf.predict(covars[test])
    pred_rfrk = pred_rf + res_hat

    results.append({
        "rf":   np.sqrt(((y[test] - pred_rf)   ** 2).mean()),
        "rfrk": np.sqrt(((y[test] - pred_rfrk) ** 2).mean()),
    })

import pandas as pd
print(pd.DataFrame(results).mean().round(3))

Reading. Hengl et al. (2018), "Random forest as a generic framework for predictive modeling of spatial and spatio-temporal variables", PeerJ 6:e5518. Meyer & Pebesma (2021) on the area of applicability. Ploton et al. (2020), Nature Communications 11, on spatial validation and overoptimism in the ecological literature.

Four fields, four different pressures

Digital soil mapping. The largest current application. Sparse expensive samples, abundant cheap covariates from remote sensing and terrain analysis — the ideal setting for regression kriging and its machine-learning descendants. AfSIS and SoilGrids are the reference programmes; both publish their methods, and reading one properly is worth a week of lectures.

Mineral resource estimation. Where the subject began, and the one field where the reporting is legally constrained. Codes such as JORC and SAMREC govern how a resource may be stated publicly, and a Competent Person signs off personally. Block kriging on the mining block support, careful treatment of the high-grade tail, and conditional simulation for grade-tonnage uncertainty are all standard. This is geostatistics at its most rigorous, because being wrong is actionable.

Environmental and public health. Contamination, air quality, disease incidence. Almost always a threshold question rather than a value question, which means indicator methods and simulation rather than plain kriging. The audience is a regulator or a health officer, so communication of uncertainty is half the work.

Climate and agriculture. Rainfall and temperature surfaces from station networks, almost always with elevation as a covariate, feeding drought indices and yield models. Kenya's station network is sparse and unevenly distributed, which makes declustering and honest error reporting matter more here than in the textbooks.

What a defensible report contains

Before it is a map, an analysis is a set of decisions someone else must be able to check. Every one of these belongs in the report:

1. The question, and the support and units of the answer.
2. The data: source, count, sampling design, date, and known quality problems.
3. Exploratory findings: distribution, transform applied and its justification, outliers investigated and the decision on each, trend, declustering.
4. The variogram: estimator, binning choices, model, parameters, a plot of the fit, and directional analysis with the anisotropy conclusion.
5. The prediction method, with the reason it was chosen over the alternatives.
6. Validation: the full diagnostic set, under spatial blocking, with a comparison against at least one simpler method.
7. Uncertainty: how it was quantified and what it does and does not cover.
8. The limitations, stated plainly — including where the map should not be used.

Point 8 separates a professional from a technician. Every map has regions where the data do not support the prediction: outside the sampling envelope, in areas of a covariate range never observed, in unsampled corners. Say where they are. A report that identifies its own limits is far more likely to be trusted than one that does not.

The one thing to carry out of this unit

You will forget the spherical model's formula; it is in every textbook and in Appendix B. What should stay with you is the habit the whole unit has been arranged around: every prediction carries an uncertainty, that uncertainty can be quantified, and quoting a number without one is not a professional act.

Software will keep making the computation easier. It has never made the judgement easier — deciding whether the data support the question, whether the model is defensible, whether the map should be trusted where someone is about to act on it. That judgement is what you were doing all semester, and it is what you are being hired for.

Reading. Hengl & MacMillan, Predictive Soil Mapping with R (open access). Rossiter (2018), "Past, present and future of information technology in pedometrics", Geoderma 324. The JORC Code (2012) for the mineral reporting context.

You met Python in GIS programming. This appendix tells you — privately, before anyone is watching — whether that is enough for what follows. Work through the twelve tasks below without help of any kind: no notes, no search, no AI assistant. The point is to find out what you actually have, and using help defeats the entire purpose.

The twelve tasks

1. Read samples.csv into a pandas DataFrame and print the number of rows and the column names.
2. Select all rows where soc is greater than 20 and depth is less than 30.
3. Add a column soc_pct equal to soc divided by 10.
4. Write a function euclid(x1, y1, x2, y2) that returns the distance between two points. Test it on a case whose answer you know.
5. Given lists xs and ys, use a loop to build a list of distances from every point to the point (0, 0).
6. Do the same with NumPy and no loop at all.
7. Create a 5×5 NumPy array of random numbers. Extract the third column. Extract the bottom-right 2×2 corner.
8. Given a 2D array D of pairwise distances, return the mean of all values strictly greater than 0 and less than 100.
9. Plot a scatter of x against y, coloured by soc, with a colourbar and axis labels.
10. This code raises IndexError: index 10 is out of bounds for axis 0 with size 10. Explain why, and fix it:
    for i in range(len(v) + 1): print(v[i])
11. Write a function that takes a NumPy array and returns a new array with values above the 95th percentile replaced by the 95th percentile. Do not modify the input.
12. Explain in one sentence what the following does, and why it is preferable to a loop:
    np.linalg.norm(pts[:, None, :] - pts[None, :, :], axis=-1)

Tasks 1–5 are the floor. Tasks 6–9 are the working level for this unit. Tasks 10–12 are what the Extension work assumes — task 12 in particular is the broadcasting idiom that appears in week 5's kriging code, and understanding it is worth more than memorising any formula in this reader.

Five things that will actually bite you

• Views versus copies. b = a[2:5] gives you a view — writing to b changes a. Use .copy() when you mean a copy. This causes silent, baffling bugs.
• Broadcasting. The rules by which NumPy combines arrays of different shapes. Learning them turns most loops into one line, and pairwise-distance computation is unusable without them.
• Chained assignment in pandas. df[df.a > 5]["b"] = 0 silently fails. Use df.loc[df.a > 5, "b"] = 0.
• Index alignment. Pandas aligns on the index when you combine Series. After filtering, indices have gaps, and combining a filtered Series with a fresh array produces NaNs. Use .values or .reset_index(drop=True) deliberately.
• Floating-point equality. Never test a == b on floats. Use np.isclose. This matters in week 5, where γ(0) must be exactly zero.

Where to go if you scored low

• NumPy — the absolute basics, in the official NumPy documentation. Two hours, and it covers tasks 6–8.
• Pandas — 10 minutes to pandas, official documentation. Longer than ten minutes, worth it.
• Python Data Science Handbook by Jake VanderPlas, chapters 2 and 3, free online. The best single explanation of broadcasting in print.
• Automate the Boring Stuff with Python, chapters 1–6, free online, if the gap is in core Python rather than the libraries.

This sheet is supplied in the CATs and in the final examination. You are not expected to memorise these; you are expected to know which one to reach for, and what each symbol means.

Structure

Permissible models

Ordinary kriging

Validation

Indicators and declustering

Rules of thumb

Maximum lag	≤ ½ the survey diagonal
Pairs per lag class	≥ 30
First lag width	≈ mean nearest-neighbour distance
Angular tolerance, directional	±22.5°
Kriging neighbourhood	16–24 points, quadrant-balanced
Nugget / sill > 0.5	weak structure — reconsider
CV block size	≥ 2 × variogram range
Realisations for tail probabilities	≥ 500
Kenya UTM zones	EPSG:32736 (W of 36°E) · EPSG:32737 (E of 36°E)

All datasets are distributed through the unit's shared folder as CSV and GeoTIFF, with a README giving the CRS, units and provenance of each. If you use them outside this unit, check the licence file — AfSIS and CHIRPS are openly licensed with attribution requirements, and the Kenyan borehole data is redistributed for teaching only.

EGS 2405 Geostatistics · Course reader, edition 1.0 · Jomo Kenyatta University of Agriculture and Technology · September 2026. Corrections and suggestions are welcome and will be credited in the next edition.