python - interpolate.griddata 结果不一致

Question

我正在将一些代码从 Matlab 转换为 Python，发现我从 scipy.interpolate.griddata 得到的结果与从 Matlab scatteredInterpolant 得到的结果不同。经过大量研究和实验，我发现 scipy.interpolate.griddata 的插值结果似乎取决于所提供数据集的大小。似乎存在导致插值改变的阈值。这是一个错误吗？或者有人可以解释会导致这种情况的算法。这是演示问题的代码。

import numpy as np
from scipy import interpolate

# This code provides a simple example showing that the interpolated value 
# for the same location changes depending on the size of the input data set.

# Results of this example show that the interpolated value changes 
# at repeat 10 and 300.

def compute_missing_value(data):
    """Compute the missing value example function."""

    # Indices for valid x, y, and z data
    # In this example x and y are simply the column and row indices
    valid_rows, valid_cols = np.where(np.isnan(data) == False)
    valid_data = data[np.isnan(data) == False]

    interpolated_value = interpolate.griddata(np.array((valid_rows, 
                   valid_cols)).T, valid_data, (2, 2), method='linear')


    print('Size=', data.shape,'  Value:', interpolated_value)


# Sample data
data = np.array([[0.2154, 0.1456, 0.1058, 0.1918],
                 [-0.0398, 0.2238, -0.0576, 0.3841],
                 [0.2485, 0.2644, 0.2639, 0.1345],
                 [0.2161, 0.1913, 0.2036, 0.1462],
                 [0.0540, 0.3310, 0.3674, 0.2862]])

# Larger data sets are created by tiling the original data.
# The location of the invalid data to be interpolated is maintained at 2,2
repeat_list =[1, 9, 10, 11, 30, 100, 300]
for repeat in repeat_list:
    new_data = np.tile(data, (1, repeat))
    new_data[2,2] = np.nan
    compute_missing_value(new_data)

结果是:

Size= (5, 4)   Value: 0.07300000000000001  
Size= (5, 36)   Value: 0.07300000000000001  
Size= (5, 40)   Value: 0.19945000000000002  
Size= (5, 44)   Value: 0.07300000000000001  
Size= (5, 120)   Value: 0.07300000000000001  
Size= (5, 400)   Value: 0.07300000000000001  
Size= (5, 1200)   Value: 0.19945000000000002

Answer 1

Jaime's answer 描述了 scipy.interpolate.griddata 如何使用 Delaunay 三角剖分对值进行插值:

[When] you make a call to scipy.interpolate.griddata:

1. First, a call to sp.spatial.qhull.Delaunay is made to triangulate the irregular grid coordinates.
2. Then, for each point in the new grid, the triangulation is searched to find in which triangle ... does it lay.
3. The barycentric coordinates of each new grid point with respect to the vertices of the enclosing simplex are computed.
4. An interpolated values is computed for that grid point, using the barycentric coordinates, and the values of the function at the vertices of the enclosing simplex.

pv. explains 德劳内正方形网格生成的三角剖分不是唯一的。由于要点根据三角测量进行线性插值，你可以得到不同的结果取决于生成的特定 Delaunay 三角剖分。

这是您的脚本的修改版本，它绘制了所用的 Delaunay 三角剖分:

import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt
import scipy.spatial as spatial
import matplotlib.collections as mcoll

def compute_missing_value(data):
    """Compute the missing value example function."""

    mask = ~np.isnan(data)
    valid_rows, valid_cols = np.where(mask)
    valid_data = data[mask]
    interpolated_value = interpolate.griddata(
        (valid_cols, valid_rows), valid_data, (2, 2), method='linear')

    print('Size: {:<12s} Value: {:<.4f}'.format(
        str(data.shape), interpolated_value))

    points = np.column_stack((valid_cols, valid_rows))

    tess = spatial.Delaunay(points)
    tri = tess.simplices 
    verts = tess.points[tri]
    lc = mcoll.LineCollection(
        verts, colors='black', linewidth=2, zorder=5)
    fig, ax = plt.subplots(figsize=(6, 6))
    ax.add_collection(lc)
    
    ax.plot(valid_cols, valid_rows, 'ko')
    ax.set(xlim=(0, 3), ylim=(0, 3))
    plt.title('Size: {:<12s} Value: {:<.4f}'.format(
        str(data.shape), interpolated_value))

    for label, x, y in zip(valid_data, valid_cols, valid_rows):
        plt.annotate(
            label,
            xy=(x, y), xycoords='data',
            xytext = (-20, -40), textcoords = 'offset points',
            horizontalalignment = 'center',
            verticalalignment = 'bottom',
            bbox = dict(
                boxstyle='round,pad=0.5', fc='yellow', alpha=0.5),
            arrowprops = dict(arrowstyle='->', connectionstyle='arc3,rad=0'))

    plt.show()


# Sample data
orig_data = np.array([[0.2154, 0.1456, 0.1058, 0.1918],
                 [-0.0398, 0.2238, -0.0576, 0.3841],
                 [0.2485, 0.2644, 0.2639, 0.1345],
                 [0.2161, 0.1913, 0.2036, 0.1462],
                 [0.0540, 0.3310, 0.3674, 0.2862]])

repeat_list =[1, 4]
for repeat in repeat_list:
    print('{}: '.format(repeat), end='')
    new_data = np.tile(orig_data, (1, repeat))
    new_data[2,2] = np.nan
    compute_missing_value(new_data)

如您所见，两个插值 0.1995 和 0.073 是 (A,C) 或 (B) 的平均值,D)(使用 pv.'s notation ):

In [159]: (0.2644+0.1345)/2
Out[159]: 0.19945000000000002

In [160]: (0.2036-0.0576)/2
Out[160]: 0.07300000000000001