🔧 Fix Interpolation & surfaces Spatial statistics Interpolation & Surfaces · step 7 of 13 Updated

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.3 m where the actual RMSE against known truth was 40.2 m — an overstatement of 68%.

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")

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.

Four kriging failures: singular matrix from duplicates, flat surface from a pure nugget, an implausible range, and miscalibrated variances.
Only the first raises. The rest produce a surface that looks like output.

Step-by-step solution

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())

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   51,039 m²    data variance 41,950 m²
range   4,537 m     largest lag 4,875 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: 4.5 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,536 m²  ->  predicted sd 67.3 m
actual RMSE                                        40.2 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.323

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

5. Negative predictions from a non-negative variable

Kriging weights can be negative — measured, every one of 116,137 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.

Standardised leave-one-out residuals with a variance of 0.323 rather than 1, indicating inflated kriging variances.
The standardised residual variance tests the variogram, not the interpolation. It is the diagnostic almost nobody runs.

Code examples

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

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
  LOO RMSE 40.29
  mean standardised residual -0.003  (want ~0)
  variance of standardised    0.323  (want ~1)
  ! variances inflated by about 1.8x — 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)

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.

Explanation

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 68% overstatement came largely from a sill 22% above the data variance, which the variance scales with directly.

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 116,137 grid targets had at least one negative weight, and the resulting surface reached 51.6 m where the lowest sample was 61.2 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.

A decision path from a flat kriged surface through the empirical variogram to no structure, a fitting failure, a trend, or a fine variogram.
Read the empirical variogram before the fitted one — the fit can only report what the points show.

Edge cases or notes

  • 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 — 68% too wide in one measurement.
  • Standardised LOO residual variance should be about 1. 0.323 means inflated variances.
  • Kriging weights can be negative, so a non-negative variable can go negative.

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.3 m where the actual RMSE was 40.2 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.323 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.