不规则网格上的插值方法

Question

我有三个numpy数组，分别存储了经纬度和一些属性值，这些值是在一个网格上，比如说我有LAT(y,x)、LON(y,x)和温度T(y,x)，这些都是在某个x和y的范围内。这个网格不一定是规则的，实际上，它是三极的。

接下来，我想把这些属性（温度）值插值到一些不同的经纬度点上，这些点存储为lat1(t)、lon1(t)，大约有10,000个这样的点，这些点并不在实际的网格点上。我试过使用matplotlib.mlab.griddata，但这个方法太慢了（毕竟它并不是为我现在的需求设计的）。我也试过scipy.interpolate.interp2d，但遇到了内存错误（我的网格大约是400x400）。

有没有什么简单又快的方法可以做到这一点？我总觉得答案应该很明显……谢谢！！

Answer 1

这里有很多选择，哪种最好取决于你的数据...不过我不知道有没有现成的解决方案适合你。

你提到你的输入数据来自三极数据。这个数据可能有三种主要的结构方式。

1. 从三维网格中采样，然后投影回二维的纬度（LAT）和经度（LON）数据。
2. 从二维网格中采样，然后直接得到二维的纬度和经度数据。
3. 在三极空间中没有结构的数据，投影成二维的纬度和经度数据。

其中最简单的是第二种。与其在纬度和经度的空间中插值，不如“直接”把你的点转换回源空间，然后在那儿进行插值。

对于第一种和第二种情况，另一个选择是寻找从三极空间映射到你的采样点的单元格。（你可以使用BSP或网格结构来加快这个搜索过程）选择其中一个单元格，然后在里面进行插值。

最后，还有很多无结构插值的选项……不过这些通常比较慢。我个人比较喜欢用最近的N个点进行线性插值，找到这N个点的方法也可以用网格或BSP来实现。另一个不错的选择是对无结构的点进行德劳内三角剖分，然后在得到的三角网格上进行插值。

如果我的网格是第一种情况，我会选择无结构的策略，因为我担心需要处理重叠投影的单元格搜索。选择“正确”的单元格会很困难。

Answer 2

这里有一个很不错的例子，来自Roger Veciana i Rovira，讲的是如何利用反距离加权的方法来处理数据。如果你对这方面感兴趣，他还提供了一些使用GDAL库将数据写入地理信息TIFF格式的代码。

当然，这个例子是针对一个常规的网格来说的，但如果你先用pyproj或者其他工具把数据投影到像素网格上，同时注意选择合适的投影方式，那么就可以使用这个方法了。

以下是他的算法和示例脚本的副本：

from math import pow  
from math import sqrt  
import numpy as np  
import matplotlib.pyplot as plt  
  
def pointValue(x,y,power,smoothing,xv,yv,values):  
    nominator=0  
    denominator=0  
    for i in range(0,len(values)):  
        dist = sqrt((x-xv[i])*(x-xv[i])+(y-yv[i])*(y-yv[i])+smoothing*smoothing);  
        #If the point is really close to one of the data points, return the data point value to avoid singularities  
        if(dist<0.0000000001):  
            return values[i]  
        nominator=nominator+(values[i]/pow(dist,power))  
        denominator=denominator+(1/pow(dist,power))  
    #Return NODATA if the denominator is zero  
    if denominator > 0:  
        value = nominator/denominator  
    else:  
        value = -9999  
    return value  
  
def invDist(xv,yv,values,xsize=100,ysize=100,power=2,smoothing=0):  
    valuesGrid = np.zeros((ysize,xsize))  
    for x in range(0,xsize):  
        for y in range(0,ysize):  
            valuesGrid[y][x] = pointValue(x,y,power,smoothing,xv,yv,values)  
    return valuesGrid  
      
  
if __name__ == "__main__":  
    power=1  
    smoothing=20  
  
    #Creating some data, with each coodinate and the values stored in separated lists  
    xv = [10,60,40,70,10,50,20,70,30,60]  
    yv = [10,20,30,30,40,50,60,70,80,90]  
    values = [1,2,2,3,4,6,7,7,8,10]  
      
    #Creating the output grid (100x100, in the example)  
    ti = np.linspace(0, 100, 100)  
    XI, YI = np.meshgrid(ti, ti)  
  
    #Creating the interpolation function and populating the output matrix value  
    ZI = invDist(xv,yv,values,100,100,power,smoothing)  
  
  
    # Plotting the result  
    n = plt.normalize(0.0, 100.0)  
    plt.subplot(1, 1, 1)  
    plt.pcolor(XI, YI, ZI)  
    plt.scatter(xv, yv, 100, values)  
    plt.title('Inv dist interpolation - power: ' + str(power) + ' smoothing: ' + str(smoothing))  
    plt.xlim(0, 100)  
    plt.ylim(0, 100)  
    plt.colorbar()  
  
    plt.show() 

Answer 3

可以试试把反距离加权和scipy.spatial.KDTree结合起来，具体的内容可以参考SO上的反距离加权插值的Python实现。Kd树在二维、三维等场景下都表现得很好，反距离加权方法则是平滑且局部的，而k=最近邻的数量可以调整，以在速度和准确性之间找到一个平衡。