为什么 SciPy 的插值网格数据对于不同大小的样本返回不同的结果？

Question

我有一个网格，其间隔均匀的点为 20x20 厘米，我想对存储在 Pandas Dataframe 中的值（

scipy.interpolate.griddata(method="cubic")

）进行插值，其中包含“x”、“y”和“value”列。我使用了 DF 的一部分进行测试，将我的脚本应用于整个 DF 后，我注意到插值结果有所不同。差别不大，但足以造成我进一步处理后整体结果的差异，这使得用样本切片快速测试参数几乎是不可能的......

样本部分位于整个数据集的点的中间，两次迭代（样本部分和整个数据集）的网格以及使用的参数都是相同的。用于样本部分的点和值与该区域的整个数据集中的点和值相同。我预计样本部分的边缘在插值后会给出不同的结果，但中间的区域不会。这种行为也会出现在不同的插值方法中，例如

linear

或使用

nearest

的

scipy.interpolate.griddata

。是否需要更改或添加一些内容来防止这种行为（例如不同的库或方法）？

这是我的代码：

import numpy as np
import pandas as pd
import geopandas as gpd
from scipy import interpolate

# Import data
df_grid = pd.read_table("whole_dataset.csv", sep=",", header=0)
df = pd.read_table("sample_section.csv", sep=",", header=0) or df = pd.read_table("whole_dataset.csv", sep=",", header=0)

# Define the regular grid for interpolation
x_min = df_grid['x'].min()
x_max = df_grid['x'].max()
y_min = df_grid['y'].min()
y_max = df_grid['y'].max()
spacing = 0.2
x_grid = np.arange(x_min, x_max + spacing, spacing)
y_grid = np.arange(y_min, y_max + spacing, spacing)

# Create a mesh grid from the regular grid
X, Y = np.meshgrid(x_grid, y_grid)

# Interpolate values using 'cubic' method
points = df[['x', 'y']].values
values = df['value'].values
interpolated_values = interpolate.griddata(points, values, (X, Y), method='cubic')

# Convert the interpolated values to a DataFrame with 'x' and 'y' coordinates
data = pd.DataFrame({
    'x': X.ravel(),
    'y': Y.ravel(),
    'value': interpolated_values.ravel()
})

# Drop all Null values
data = data.dropna().reset_index(drop=True)

gdf = gpd.GeoDataFrame(data, geometry=gpd.points_from_xy(data['x'], data['y']))
gdf_combined.to_file("data.geojson", driver='GeoJSON', crs="EPSG:32633")
        

这是我得到的结果的差异：

更新：

这是两个 CSV 文件

whole_dataset.csv

和

sample_section.csv

(Google Drive)。使用

QGIS

提取 sample_section，以便从数据集的空间中心获取一个区域，然后保存为 CSV。坐标的CRS是WGS 84 / UTM zone 33N - EPSG:32633。

GIF 也是来自 QGIS 的屏幕截图，其中插值网格的两层用相同的符号覆盖（QGIS 样式文件 (qml) 也在 Google Drive 中）

Answer 1

我有一个网格，其间隔均匀的点为 20x20 厘米，我想对存储在 Pandas Dataframe 中的值（

scipy.interpolate.griddata(method="cubic")

）进行插值，其中包含“x”、“y”和“value”列。

我认为这里有两个问题：一个是这里使用griddata的理论问题，另一个是数值稳定性问题。

对来自直线网格的输入数据使用 griddata

当 SciPy 插入这些数据时，它首先确定哪些点是相关的。它通过查找每个点之间的三角形，并通过查找三角形上的三个点并在它们之间进行插值来在点之间进行插值来实现此目的。

griddata()

不允许您直接检查此三角剖分，但您可以调用它所包装的类，并直接获取信息。

interpolator_sample = CloughTocher2DInterpolator(points, values, fill_value=np.nan, rescale=rescale)
interpolated_values = interpolator_sample((X, Y))
scipy.spatial.delaunay_plot_2d(interpolator_sample.tri)
plt.xlim([516207.788138, 516232.388138])
plt.ylim([5.533969e+06, 5.533987e+06])

剧情：

您可以将同样的事情应用于完整的数据集：

对此我想指出三点。

1. 它们不是相同的三角形。 （如果您对此不相信，我建议在两个图上使用“在选项卡中打开图像”并在它们之间翻转。）当样本数据集插值器和完整数据集插值器在数据点之间插值时，它们是在不同数据点之间插值。

2. 三角剖分包含许多太小的三角形。 （它们在此图上不可见，但有角度为 0、0、180 的三角形连接了许多这些点。）Delaunay 三角剖分应该避免创建任何角度接近零的三角形，但上面的三角剖分包含许多非常窄的三角形。

3. 输入点位于旋转网格上。这意味着该数据有多个可能的三角剖分，这些三角剖分同样有效。

   这个问题很好地解释了原因：interpolate.griddata 结果不一致

   我相信即使网格像这样旋转/倾斜，这个论点仍然成立。这意味着在添加更多数据时无法使

   griddata()

   具有确定性，并且任何使用它的方法都注定要失败。

数值稳定性

我尝试了一些想法来改进这一点，同时使用下图衡量成功。

data_merged = data_interp_sample.merge(data_interp_full, on=['x', 'y'], suffixes=['_sample', '_full'])
ax = data_merged.plot.scatter(x='value_sample', y='value_full', alpha=0.1)
ax.axline((0, 0), slope=1, color='r')

该图显示了样本插值和完整数据集插值的所有常见值。红线沿着 y = x，它代表理想情况：每个点在两次插值中都完全相同。

剧情：

我尝试了几种方法来改进这一点，但一件事似乎出人意料地有用，那就是减去每个 XY 值的平均值，然后在插值后将其加回来。

df_grid = pd.read_table("whole_dataset.csv", sep=",", header=0)
df = pd.read_table("sample_section.csv", sep=",", header=0)
df_grid_mean_x = df_grid['x'].mean()
df_grid_mean_y = df_grid['y'].mean()
df['x'] = df['x'] - df_grid_mean_x
df_grid['x'] = df_grid['x'] - df_grid_mean_x
df['y'] = df['y'] - df_grid_mean_y
df_grid['y'] = df_grid['y'] - df_grid_mean_y

# ... snip unchanged code ...

# Convert the interpolated values to a DataFrame with 'x' and 'y' coordinates
data = pd.DataFrame({
    'x': X.ravel() + df_grid_mean_x,
    'y': Y.ravel() + df_grid_mean_y,
    'value': interpolated_values.ravel()
})

这应该没有什么区别 - 即使输入坐标被任何常量平移，德劳尼三角剖分也应该完全相同。然而，它似乎产生了更合理的三角测量。我相信这可能是由 Qhull 内部的浮点不精确误差引起的，减去平均值可以使其更精确地计算三角形角度。这可以避免创建小三角形。

因此，这个三角测量在样本和完整之间更加相似，但并不完全一致。