Kriging（高斯过程回归）误差过高且插值可视化异常，空间插值结果不符预期求助

Hey there! Let’s work through the issues you’re having with your Kriging interpolation—no color in your output map and sky-high prediction errors are both fixable with a few key adjustments to your code. Let’s break it down step by step:

1. Fix Your Grid Generation (The Root Cause of Most Issues) #

Looking at how you’re creating your interpolation grid:

xx, yy = np.mgrid[6.67:6.69:6.72j, 50.86:50.88:50.92j]

This has two critical problems:

• The grid range doesn’t cover your data: Your sample data has longitudes around 6.70 and latitudes around 50.91, but your grid only spans 6.67–6.69 (longitude) and 50.86–50.88 (latitude). You’re interpolating in a completely empty area with no data points to support predictions—hence the massive errors and blank-looking map.
• The j syntax is misused: In np.mgrid, start:stop:NUMBERj means "generate NUMBER points between start and stop". 6.72j would create 672 points, but since your start/stop range is tiny, this is overkill and irrelevant.

Fix this by generating a grid that actually covers your data’s geographic range:

# Get the min/max of your data, add a small buffer
lon_min, lon_max = data.longitude.min() - 0.01, data.longitude.max() + 0.01
lat_min, lat_max = data.latitude.min() - 0.01, data.latitude.max() + 0.01

# Generate a 100x100 grid (adjust the number of points as needed)
xx, yy = np.mgrid[lon_min:lon_max:100j, lat_min:lat_max:100j]

2. Align Your Map Axes to Real Coordinates #

Your current matshow call uses pixel indices for axes instead of actual lat/long, which makes your set_xlim/set_ylim ineffective. Add the extent parameter to matshow to map the grid to real geographic coordinates:

# For the interpolation map
art = axes[0].matshow(field.T, origin='lower', cmap='plasma', extent=[lon_min, lon_max, lat_min, lat_max])

# For the error map
art = axes[1].matshow(s2.T, origin='lower', cmap='YlGn_r', extent=[lon_min, lon_max, lat_min, lat_max])

This will also make your data points (+k/+w) align correctly with the interpolation grid.

3. Tune Your Variogram for Better Kriging #

Your Variogram’s maxlag=5 is way too large for your data—your coordinates are in decimal degrees, where a difference of 5 would span hundreds of kilometers, way beyond your local dataset. This leads to a poorly fitted variogram, which ruins Kriging predictions.

Adjust the maxlag to a meaningful value based on your data’s spatial extent:

from scipy.spatial.distance import pdist

# Calculate the maximum distance between any two data points
max_dist = pdist(data[['longitude', 'latitude']]).max()

# Set maxlag to half that distance (a common rule of thumb)
V = Variogram(
    data[['longitude', 'latitude']].values, 
    data.subsidence, 
    normalize=False, 
    maxlag=max_dist/2,  # This is key
    n_lags=15
)

You can also add a radius parameter to your OrdinaryKriging call to limit predictions to points within the variogram’s effective range (0.05, per your output):

ok = OrdinaryKriging(
    V, 
    min_points=5, 
    max_points=15, 
    mode='exact',
    radius=V.describe()['effective_range']  # Only use nearby points
)

4. Clean Up Redundant Code #

You’re reading your CSV twice (data = pd.read_csv("kerpencsv.csv") appears twice)—remove the duplicate line to keep things clean.

Full Corrected Code Snippet #

Here’s how your final code should look with these fixes:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from pprint import pprint
from scipy.spatial.distance import pdist
plt.style.use('ggplot')
from skgstat import Variogram, OrdinaryKriging

# Load data once
data = pd.read_csv("kerpencsv.csv")
print("Loaded %d rows and %d columns" % data.shape)
print(data.head())

# Plot data points
fig, ax = plt.subplots(1, 1, figsize=(9, 9))
art = ax.scatter(data.longitude, data.latitude, s=50, c=data.subsidence, cmap='plasma')
plt.colorbar(art)
plt.show()

# Calculate meaningful maxlag for variogram
max_dist = pdist(data[['longitude', 'latitude']]).max()
V = Variogram(
    data[['longitude', 'latitude']].values, 
    data.subsidence, 
    normalize=False, 
    maxlag=max_dist/2, 
    n_lags=15
)
fig = V.plot()
print('Sample variance: %.2f Variogram sill: %.2f' % (data.subsidence.var(), V.describe()['sill']))
pprint(V.describe())
print(V)

# Generate proper interpolation grid
lon_min, lon_max = data.longitude.min() - 0.01, data.longitude.max() + 0.01
lat_min, lat_max = data.latitude.min() - 0.01, data.latitude.max() + 0.01
xx, yy = np.mgrid[lon_min:lon_max:100j, lat_min:lat_max:100j]

# Kriging with radius constraint
ok = OrdinaryKriging(V, min_points=5, max_points=15, mode='exact', radius=V.describe()['effective_range'])
field = ok.transform(xx.flatten(), yy.flatten()).reshape(xx.shape)
s2 = ok.sigma.reshape(xx.shape)

# Plot results with correct axes
fig, axes = plt.subplots(1, 2, figsize=(20, 10))

# Interpolation map
art = axes[0].matshow(field.T, origin='lower', cmap='plasma', extent=[lon_min, lon_max, lat_min, lat_max])
axes[0].set_title('Kriging Interpolation')
axes[0].plot(data.longitude, data.latitude, '+k')
plt.colorbar(art, ax=axes[0])

# Error map
art = axes[1].matshow(s2.T, origin='lower', cmap='YlGn_r', extent=[lon_min, lon_max, lat_min, lat_max])
axes[1].set_title('Kriging Prediction Error')
plt.colorbar(art, ax=axes[1])
axes[1].plot(data.longitude, data.latitude, '+w')

plt.show()

Why This Works #

• Your grid now covers the area where your data lives, so Kriging has actual points to interpolate between—this will fix the "no color" issue and drastically reduce prediction errors.
• The variogram is now fitted to relevant spatial scales, ensuring the Kriging model uses meaningful spatial correlation.
• The extent parameter makes your map axes match real coordinates, so you can properly visualize how the interpolation aligns with your data points.

内容的提问来源于stack exchange，提问作者simonericmoon