重复Scipy的griddata [英] Repeating Scipy's griddata

问题描述

当数据集很多时，使用Scipy的griddata将数据(d)划分为不规则网格(x和y)是费时的.但是，经度和纬度(x和y)始终相同，仅数据(d)在变化.在这种情况下，一旦使用了giddata，如何以不同的重复次数重复该过程以获得更快的结果?

The griding the data (d) in irregular grid (x and y) using Scipy's griddata is timecomsuing when the datasets are many. But, the longitudes and latitudes (x and y) are always same, only the data (d) are changing. In this case, once using the giddata, how to repeat the procedure with different d arrys to achieve faster result?

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的一个数组. 但是我还有其他数百个数组.

The above code works for one array of d. But I have hundreds of other arrays.

推荐答案

griddata与nearest一起使用NearestNDInterpolator结束.这是一个创建迭代器的类，该迭代器使用xi:

griddata with nearest ends up using NearestNDInterpolator. That's a class that creates an iterator, which is called with the xi:

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

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

So you could create your own NearestNDInterpolator and call it with multiple times with different xi.

但是我认为您想更改values.查看我看到的该类的代码

But I think in your case you want to change the values. Looking at the code for that class I see

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

__call__可以:

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

我不知道创建tree与query的相对成本.

I don't know the relative cost of creating the tree versus query.

因此，在使用__call__时应轻松更改values.看起来values可能有多个列，因为它只是在第一维上建立索引.

So it should be easy to change values between uses of __call__. And it looks like values could have multiple columns, since it's just indexing on the 1st dimension.

此插值器非常简单，您可以使用相同的tree想法编写自己的插值器.

This interpolator is simple enough that you could write your own using the same tree idea.

这是一个最近的插值器，可让您对相同的点但不同的z值重复插值.我还没有计时，所以可以节省多少时间

Here's a Nearest Interpolator that lets you repeat the interpolation for the same points, but different z values. I haven't done timings yet to see how much time it saves

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返回最接近的值.

So as long as z has the same shape as before (and as xreg), z[I] will return the nearest value for each xi.

它还可以插补2d数据(例如(225，n)形)

And it can interpolated 2d data as well (e.g. (225,n) shaped)

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

这篇关于重复Scipy的griddata的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！