python - 重复 Scipy 的 griddata

Question

当数据集很多时，使用 Scipy 的 griddata 将数据 (d) 网格化为不规则网格(x 和 y)非常耗时。但是，经度和纬度(x 和 y)始终相同，只有数据 (d) 发生变化。在这种情况下，一旦使用了 giddata，如何使用不同的数组重复该过程以获得更快的结果？

import numpy as np, matplotlib.pyplot as plt
from scipy.interpolate import griddata

x = np.array([110, 112, 114, 115, 119, 120, 122, 124]).astype(float)

y = np.array([60, 61, 63, 67, 68, 70, 75, 81]).astype(float)

d = np.array([4, 6, 5, 3, 2, 1, 7, 9]).astype(float)

ulx, lrx = np.min(x), np.max(x)
uly, lry = np.max(y), np.min(y)
xi = np.linspace(ulx, lrx, 15)
yi = np.linspace(uly, lry, 15)
grided_data = griddata((x, y), d, (xi.reshape(1,-1), yi.reshape(-1,1)), method='nearest',fill_value=0)
plt.imshow(grided_data)
plt.show()

以上代码适用于 d 的一个数组。但我还有数百个其他数组。

Answer 1

griddata 与 nearest 最终使用 NearestNDInterpolator。这是一个创建迭代器的类，使用 xi 调用该迭代器:

elif method == 'nearest':
    ip = NearestNDInterpolator(points, values, rescale=rescale)
    return ip(xi)

因此，您可以创建自己的 NearestNDInterpolator 并使用不同的 xi 多次调用它。

但我认为在您的情况下您想要更改值。查看该类的代码，我发现

    self.tree = cKDTree(self.points)
    self.values = y

__call__ 的作用是:

    dist, i = self.tree.query(xi)
    return self.values[i]

我不知道创建树与查询的相对成本。

因此，在使用 __call__ 之间更改值应该很容易。看起来 values 可以有多个列，因为它只是在第一维上建立索引。

这个插值器非常简单，您可以使用相同的树思想编写自己的插值器。

<小时/>

这是一个最近的插值器，可让您对相同的点但不同的 z 值重复插值。我还没有进行计时，看看它节省了多少时间

class MyNearest(interpolate.NearestNDInterpolator):
    # normal interpolation, but returns the near neighbor indices as well
    def __call__(self, *args):
        xi = interpolate.interpnd._ndim_coords_from_arrays(args, ndim=self.points.shape[1])
        xi = self._check_call_shape(xi)
        xi = self._scale_x(xi)
        dist, i = self.tree.query(xi)
        return i, self.values[i]

def my_griddata(points, values, method='linear', fill_value=np.nan,
             rescale=False):
    points = interpolate.interpnd._ndim_coords_from_arrays(points)

    if points.ndim < 2:
        ndim = points.ndim
    else:
        ndim = points.shape[-1]
    assert(ndim==2)
    # simplified call for 2d 'nearest'
    ip = MyNearest(points, values, rescale=rescale)
    return ip # ip(xi)  # return iterator, not values

ip = my_griddata((xreg, yreg), z,  method='nearest',fill_value=0)
print(ip)
xi = (xi.reshape(1,-1), yi.reshape(-1,1))
I, data = ip(xi)
print(data.shape)
print(I.shape)
print(np.allclose(z[I],data))
z1 = xreg+yreg  # new z data
data = z1[I]  # should show diagonal color bars

因此，只要 z 具有与之前相同的形状(并且与 xreg 相同)，z[I] 将返回最接近的值对于每个 xi。

它也可以插值二维数据(例如(225，n)形状)

z1 = np.array([xreg+yreg, xreg-yreg]).T
print(z1.shape)   # (225,2)
data = z1[I]
print(data.shape)  # (20,20,2)