My Kriging Fit Fails or Returns a Flat Surface

Problem statement

Kriging fails in four distinct ways, and only one of them raises an exception:

• a singular matrix, which raises
• a flat surface, which does not
• a fitted range or sill that is nonsense, which does not
• variances that bear no relation to the actual error, which does not

The last one is measurable. On a real dataset, kriging's own variance implied a standard deviation of 67.1 m where the actual RMSE against known truth was 40.9 m — an overstatement of 64%.

Quick answer

Check the variogram before blaming the kriging:

print(f"nugget {nugget:,.0f}  sill {sill:,.0f}  range {rng:,.0f} m")
print(f"data variance {values.var():,.0f}")
print(f"largest lag   {lags.max():,.0f} m")

if nugget / (nugget + sill) > 0.9:
    print("! nearly all nugget — no spatial structure, kriging will be flat")
if rng > lags.max():
    print("! range beyond the data — treat it as a lower bound")
if rng < lags.min():
    print("! range below the shortest lag — every sample is independent")

nugget 0  sill 42,519  range 3,762 m
data variance 38,896
largest lag   6,569 m

A fitted range shorter than the closest sample pair means the model believes no two samples are related, and kriging degenerates to the global mean everywhere.

Duplicates, pure nugget, implausible fits, negatives

1. Singular matrix: duplicate coordinates

Two samples at the same location make the kriging matrix singular, because two rows are identical.

from scipy.spatial import cKDTree
distances, _ = cKDTree(points).query(points, k=2)
duplicates = int((distances[:, 1] < 1e-6).sum())
print(f"duplicate coordinates: {duplicates}")

duplicate coordinates: 0

Aggregate them — mean, median, or the most recent — before kriging. Using a pseudo-inverse hides the problem and produces weights that are arbitrary between the duplicated pair.

Near-duplicates cause the same problem more subtly: a matrix that is not singular but is badly conditioned, giving weights of ±10⁶ that cancel.

2. Flat surface: a pure nugget model

If the fitted nugget is most of the sill, the model says the field has no spatial structure. Kriging then correctly returns the mean everywhere.

That is not a bug, it is a finding: either the variable genuinely has no structure at the sampled scale, or the sampling is too sparse to see it.

The check is the empirical variogram. If it is flat from the shortest lag, no interpolation method will help.

3. Implausible range or sill

nugget      0 m²
sill   42,519 m²    data variance 38,896 m²
range   3,762 m     largest lag 6,569 m

Two checks. A range beyond the largest lag means the optimiser placed the plateau where no data constrains it — treat it as a lower bound.

A sill above the data variance is normal when the range is a large fraction of the study extent, as here: 3.8 km inside a 10 km window. It is not a fitting error and removing a trend does not fix it.

4. Miscalibrated variances

mean kriging variance   4,498 m²  ->  predicted sd 67.1 m
actual RMSE                                        40.9 m

Kriging variance depends only on the sample geometry and the variogram model — never on the measured values. So it inherits every error in the variogram and cannot detect that the model is wrong.

The diagnostic is the variance of the standardised leave-one-out residuals, which should be about 1:

LOO standardised residual variance: 0.288

Well below 1 means the variances are inflated, which is exactly the 64% overstatement above.

5. Negative predictions from a non-negative variable

Kriging weights can be negative — measured, every one of 110,889 grid targets had at least one negative weight. That lets the surface leave the range of the data, and a concentration or rainfall surface can come out negative.

Transform before kriging (log, logit) and back-transform after, or constrain the weights to be non-negative.

A pre-flight check and the residual test

The data used on this page

The surface is real: the Copernicus DEM GLO-30 over a 10 × 10 km window centred on the summit of Snowdon (Yr Wyddfa), warped to a 30 m grid in EPSG:32630 (x 422,880–432,870, y 5,875,440–5,885,430; 333 × 333 cells, 99.8 km²). The 500 samples are cell centres drawn at random with a fixed seed; 12,000 further cells, none of them samples, are held out as test_points/test_values. grid_points/truth are every cell centre and its elevation, and truth_fn reads the DEM bilinearly at any location, which the error tables score against. idw is the function from Spatial interpolation explained, and surface is its power-2, 12-neighbour surface on that 30 m grid. The setup also fits the spherical variogram (nugget, sill — the partial sill — and rng, from 15 lag bins in lags) and defines krige, ordinary kriging on the 24 nearest samples, which returns predictions, kriging variances and, optionally, the weights.

import os

import numpy as np
import rasterio
from rasterio.enums import Resampling
from rasterio.transform import from_origin
from rasterio.warp import reproject
from scipy.interpolate import RegularGridInterpolator

bounds = (422_880.0, 5_875_440.0, 432_870.0, 5_885_430.0)
transform = from_origin(bounds[0], bounds[3], 30, 30)
shape = (333, 333)
if not os.path.exists("snowdon_dem_utm30.tif"):
    COP = ("https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N53_00_{0}_00_DEM/"
           "Copernicus_DSM_COG_10_N53_00_{0}_00_DEM.tif")
    dem = np.full(shape, np.nan, "float32")
    for tile in ("W005", "W004"):
        part = np.full(shape, np.nan, "float32")
        with rasterio.open(COP.format(tile)) as src:
            reproject(rasterio.band(src, 1), part, dst_transform=transform, dst_crs="EPSG:32630",
                      dst_nodata=np.nan, resampling=Resampling.bilinear)
        dem = np.where(np.isnan(dem), part, dem)
    with rasterio.open("snowdon_dem_utm30.tif", "w", driver="GTiff", width=shape[1], height=shape[0], count=1,
                       dtype="float32", crs="EPSG:32630", transform=transform, nodata=np.nan) as dst:
        dst.write(dem, 1)
dem = rasterio.open("snowdon_dem_utm30.tif").read(1)

xs = bounds[0] + (np.arange(shape[1]) + 0.5) * 30
ys = bounds[3] - (np.arange(shape[0]) + 0.5) * 30
gx, gy = np.meshgrid(xs, ys)
grid_points = np.column_stack([gx.ravel(), gy.ravel()])
truth = dem.ravel().astype("float64")

rng = np.random.default_rng(0)
order = rng.permutation(truth.size)
sample_idx, test_idx = order[:500], order[500:12_500]
points, values = grid_points[sample_idx], truth[sample_idx]
test_points, test_values = grid_points[test_idx], truth[test_idx]
_bilinear = RegularGridInterpolator((ys[::-1], xs), dem[::-1], bounds_error=False, fill_value=None)
truth_fn = lambda xy: _bilinear(np.column_stack([xy[:, 1], xy[:, 0]]))
from scipy.spatial import cKDTree


def idw(points, values, targets, power=2, k=12):
    tree = cKDTree(points)
    distance, index = tree.query(targets, k=min(k, len(points)))
    distance, index = distance.reshape(len(targets), -1), index.reshape(len(targets), -1)
    weights = 1.0 / np.maximum(distance, 1e-9) ** power
    predicted = (weights * values[index]).sum(1) / weights.sum(1)
    exact = distance[:, 0] < 1e-9
    predicted[exact] = values[index[exact, 0]]
    return predicted, distance[:, 0]


surface = idw(points, values, grid_points)[0].reshape(shape)

from scipy.optimize import curve_fit
from scipy.spatial.distance import pdist

_d, _g = pdist(points), 0.5 * pdist(values[:, None], "sqeuclidean")
_edges = np.linspace(0, _d.max() / 2, 16)
lags = np.array([_d[(_d >= a) & (_d < b)].mean() for a, b in zip(_edges[:-1], _edges[1:])])
semivariance = np.array([_g[(_d >= a) & (_d < b)].mean() for a, b in zip(_edges[:-1], _edges[1:])])
_sph = lambda h, n, s, r: n + s * np.where(h < r, 1.5 * h / r - 0.5 * (h / r) ** 3, 1.0)
(nugget, sill, rng), _ = curve_fit(_sph, lags, semivariance, p0=[0, semivariance.max(), _d.max() / 4],
                                  bounds=([0, 0, 1], [semivariance.max(), 3 * semivariance.max(), _d.max()]),
                                  maxfev=20000)


def krige(points, values, targets, params, k=24, return_weights=False):
    """Local ordinary kriging with a spherical model; prediction, variance (and weights)."""
    n_, s_, r_ = params
    gamma = lambda h: np.where(h > 0, n_ + s_ * np.where(h < r_, 1.5 * h / r_ - 0.5 * (h / r_) ** 3, 1.0), 0.0)
    k = min(k, len(points))
    dist, idx = cKDTree(points).query(targets, k=k)
    dist, idx = dist.reshape(len(targets), -1), idx.reshape(len(targets), -1)
    P = points[idx]                                                   # (t, k, 2)
    A = np.ones((len(targets), k + 1, k + 1))
    A[:, :k, :k] = gamma(np.linalg.norm(P[:, :, None, :] - P[:, None, :, :], axis=-1))
    A[:, k, k] = 0.0
    b = np.ones((len(targets), k + 1))
    b[:, :k] = gamma(dist)
    solution = np.linalg.solve(A, b[..., None])[..., 0]
    weights, mu = solution[:, :k], solution[:, k]
    predicted = (weights * values[idx]).sum(1)
    variance = (weights * b[:, :k]).sum(1) + mu
    return (predicted, variance, weights) if return_weights else (predicted, variance)
print(f"DEM {np.nanmin(dem):.0f} to {np.nanmax(dem):.0f} m; {len(points)} samples, {len(test_points):,} held out")

DEM 52 to 1072 m; 500 samples, 12,000 held out

Example 1 — a pre-flight check before kriging

import numpy as np
from scipy.spatial import cKDTree


def kriging_preflight(points, values, params, lags, extent):
    """Everything that makes kriging fail, checked before it runs."""
    nugget, sill, rng = params
    problems = []

    distances, _ = cKDTree(points).query(points, k=2)
    duplicates = int((distances[:, 1] < 1e-6).sum())
    close = int((distances[:, 1] < rng / 1000).sum())
    if duplicates:
        problems.append(f"{duplicates} duplicate coordinates — the matrix "
                        "will be singular")
    if close > duplicates:
        problems.append(f"{close - duplicates} pairs closer than range/1000 — "
                        "the matrix will be ill-conditioned")

    nugget_share = nugget / max(nugget + sill, 1e-12)
    print(f"  nugget {nugget:,.0f} ({nugget_share:.1%} of the sill)")
    print(f"  sill   {sill:,.0f}   data variance {values.var():,.0f} "
          f"({sill / max(values.var(), 1e-12):.2f}x)")
    print(f"  range  {rng:,.0f} m   largest lag {lags.max():,.0f} m   "
          f"extent {extent:,.0f} m")

    if nugget_share > 0.9:
        problems.append("nearly all nugget — the surface will be flat")
    if rng > lags.max():
        problems.append("range beyond the largest lag — a lower bound only")
    if rng < distances[:, 1].min():
        problems.append("range below the closest sample pair — every sample "
                        "is independent and kriging returns the mean")
    if rng > 0.25 * extent:
        print("  note: the range is a large fraction of the extent, so the "
              "sill legitimately exceeds the data variance")

    for problem in problems:
        print(f"  ! {problem}")
    return problems


extent = bounds[2] - bounds[0]
_ = kriging_preflight(points, values, (nugget, sill, rng), lags, extent)

  nugget 0 (0.0% of the sill)
  sill   42,519   data variance 38,896 (1.09x)
  range  3,762 m   largest lag 6,569 m   extent 9,990 m
  note: the range is a large fraction of the extent, so the sill legitimately exceeds the data variance

Running this before the kriging turns four silent failures into four messages. Three of them cannot be diagnosed from the output surface at all.

Example 2 — the standardised residual test

import numpy as np


def variogram_calibration(points, values, params, krige, k=24):
    """Leave-one-out standardised residuals test the variogram, not the fit."""
    n = len(points)
    predicted = np.empty(n)
    variance = np.empty(n)

    for i in range(n):
        keep = np.ones(n, bool)
        keep[i] = False
        p, v = krige(points[keep], values[keep], points[i:i + 1], params, k=k)
        predicted[i], variance[i] = p[0], v[0]

    residual = values - predicted
    standardised = residual / np.sqrt(np.maximum(variance, 1e-12))

    rmse = float(np.sqrt(np.mean(residual ** 2)))
    print(f"  LOO RMSE {rmse:.2f}")
    print(f"  mean standardised residual {standardised.mean():+.3f}  (want ~0)")
    print(f"  variance of standardised   {standardised.var():.3f}  (want ~1)")

    if standardised.var() < 0.7:
        factor = 1 / np.sqrt(standardised.var())
        print(f"  ! variances inflated by about {factor:.1f}x — the sill is "
              "probably overestimated")
    elif standardised.var() > 1.4:
        print("  ! variances too small — the model is overconfident")
    return standardised


_ = variogram_calibration(points, values, (nugget, sill, rng), krige)

  LOO RMSE 39.43
  mean standardised residual +0.001  (want ~0)
  variance of standardised   0.288  (want ~1)
  ! variances inflated by about 1.9x — the sill is probably overestimated

This is the test that catches a wrong variogram. The predictions can be excellent while the variances are useless, and no other diagnostic separates the two.

Example 3 — conditioning fixes

import numpy as np


def condition_samples(points, values, tolerance=1e-6, aggregate="median"):
    """Collapse duplicate and near-duplicate locations."""
    from scipy.spatial import cKDTree

    tree = cKDTree(points)
    pairs = tree.query_pairs(r=tolerance, output_type="ndarray")
    if len(pairs) == 0:
        print("  no duplicate locations")
        return points, values

    parent = np.arange(len(points))

    def find(i):
        while parent[i] != i:
            parent[i] = parent[parent[i]]
            i = parent[i]
        return i

    for a, b in pairs:
        ra, rb = find(a), find(b)
        if ra != rb:
            parent[rb] = ra

    groups = {}
    for i in range(len(points)):
        groups.setdefault(find(i), []).append(i)

    reducer = {"median": np.median, "mean": np.mean}[aggregate]
    new_points, new_values = [], []
    for members in groups.values():
        new_points.append(points[members].mean(axis=0))
        new_values.append(reducer(values[members]))

    print(f"  {len(points):,} samples -> {len(new_points):,} after "
          f"collapsing duplicates ({aggregate})")
    return np.array(new_points), np.array(new_values)


_ = condition_samples(points, values)
_ = condition_samples(np.vstack([points, points[:3]]), np.concatenate([values, values[:3] + 1]))

  no duplicate locations
  503 samples -> 500 after collapsing duplicates (median)

Example 4 — the calibration and weight numbers in the steps

pred, var = krige(points, values, test_points, (nugget, sill, rng))
print(f"mean kriging variance {var.mean():,.0f} m²  ->  predicted sd {np.sqrt(var.mean()):.1f} m")
print(f"actual RMSE on the held-out cells          {np.sqrt(np.mean((pred - test_values) ** 2)):.1f} m")
_, _, w = krige(points, values, grid_points, (nugget, sill, rng), return_weights=True)
print(f"grid targets with at least one negative weight: {int((w < 0).any(1).sum()):,} of {len(w):,}")
grid_pred, _ = krige(points, values, grid_points, (nugget, sill, rng))
print(f"kriged surface minimum {grid_pred.min():.1f} m; lowest sample {values.min():.1f} m; "
      f"sill / data variance {sill / values.var():.2f}")
for lag, observed in zip(lags[:3], semivariance[:3]):
    model = nugget + sill * (1.5 * lag / rng - 0.5 * (lag / rng) ** 3)
    print(f"lag bin {lag:5.0f} m: samples {observed:7,.0f} m², model {model:7,.0f} m² ({model / observed:.1f}x)")

mean kriging variance 4,498 m²  ->  predicted sd 67.1 m
actual RMSE on the held-out cells          40.9 m
grid targets with at least one negative weight: 110,889 of 110,889
kriged surface minimum 51.3 m; lowest sample 53.4 m; sill / data variance 1.09
lag bin   300 m: samples   2,766 m², model   5,073 m² (1.8x)
lag bin   703 m: samples  11,744 m², model  11,779 m² (1.0x)
lag bin  1150 m: samples  20,115 m², model  18,884 m² (0.9x)

Collapsing with a union-find handles chains of near-duplicates correctly — three samples each within tolerance of the next collapse to one, which pairwise deduplication would get wrong.

Singular matrices and flat surfaces explained

Why duplicates make the matrix singular

The kriging system's left-hand matrix holds the variogram between every pair of samples. Two samples at the same location have identical rows and identical columns, so the matrix has linearly dependent rows and no unique inverse.

np.linalg.pinv returns a pseudo-inverse, which produces a solution — one that splits the weight arbitrarily between the duplicated pair. The prediction is usually reasonable and the weights are meaningless.

Near-duplicates are worse in practice: the matrix is invertible but badly conditioned, and the solved weights can be enormous with opposite signs, cancelling to a plausible answer with no numerical stability.

Why a pure nugget gives a flat surface

The nugget is variance present at zero separation. If the fitted model is nearly all nugget, it says two samples one metre apart are as different as two samples ten kilometres apart.

Under that model the optimal prediction anywhere is the global mean, and kriging correctly returns it. The flat surface is the right answer to the model it was given.

Whether the model is right is a separate question. A large nugget can mean measurement error, variation at scales below the sampling, or a variogram fitted to too few pairs at short lags.

Why the variance is not an error estimate

Look at what enters the kriging system. The left-hand matrix holds variogram values between samples; the right-hand vector holds variogram values from samples to the target. Neither contains a measured value.

So the variance is a function of geometry and the variogram model alone. Two datasets with identical sample layouts and completely different values produce identical variance maps.

That has one useful consequence — the variance map can be computed before collecting data, to design a survey — and one awkward one: it cannot detect that the variogram is wrong, and it inherits every error in it.

The measured 64% overstatement is not a sill problem here — the sill is only 9% above the data variance. It is the shape near the origin: in the first lag bin, 300 m, the fitted spherical model gives 5,073 m² where the samples show 2,766 m², 1.8 times as much, and it is the short lags that set each prediction's variance (Example 4).

Why negative weights are normal and consequential

The kriging system minimises expected squared error, with no constraint that weights be positive. Where a sample sits directly behind a nearer one, the system gives the far one a small negative weight — the screen effect.

Measured, all 110,889 grid targets had at least one negative weight, and the resulting surface reached 51.3 m where the lowest sample was 53.4 m.

For elevation that is harmless. For a concentration, a rainfall total or a count, a negative prediction is impossible, and the fix is a transform before kriging or a non-negativity constraint on the weights.

pinv, ranges beyond the lag and LOO variance

• Duplicate coordinates make the matrix singular. Aggregate them.
• pinv hides the problem rather than solving it.
• A pure nugget gives a flat surface, correctly.
• A range beyond the largest lag is a lower bound, not a measurement.
• A sill above the data variance is normal when the range approaches the domain size.
• Kriging variance is geometry, not accuracy — 64% too wide in one measurement.
• Standardised LOO residual variance should be about 1. 0.288 means inflated variances.
• Kriging weights can be negative, so a non-negative variable can go negative.

Internal links

• How to fit a variogram and krige a surface in Python — the working method
• The variogram explained — nugget, sill and range
• IDW, kriging, splines or TIN? Choosing an interpolator — when to use something simpler
• Spatial interpolation explained — the error-versus-distance picture
• How to cross-validate an interpolated surface — validating the predictions
• My interpolated surface extends far beyond the data — masking the output
• My IDW surface is full of bullseyes — the nugget as principled smoothing
• Sample design explained: where to measure — using the variance map to plan

FAQ

Why is my kriging matrix singular?

Almost always duplicate sample coordinates, which produce identical rows. Aggregate them rather than reaching for a pseudo-inverse.

Why is my kriged surface completely flat?

The fitted variogram is nearly all nugget, so the model says there is no spatial structure and the optimal prediction everywhere is the mean.

Why is my fitted sill larger than the data variance?

Usually because the range is a large fraction of the study extent. It is expected, and removing a trend does not fix it.

Why does the kriging variance not match my actual error?

Because it depends only on sample geometry and the variogram, never on the values. In one measurement it implied 67.1 m where the actual RMSE was 40.9 m.

How do I check the variogram is right?

Compute leave-one-out standardised residuals — residual divided by kriging standard deviation. Their variance should be about 1; 0.288 means the variances are inflated.

Why did kriging produce a negative value?

Kriging weights can be negative, so the surface is not bounded by the data. Transform a non-negative variable before kriging, or constrain the weights.

Should I use a pseudo-inverse?

Only as a safety net. If it is doing real work, the neighbourhood contains duplicates and the weights it returns are arbitrary.