在不规则网格上插值

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,但是我得到了一个MemoryError(我的网格大约是400x400).

有没有任何光滑的,最好快速的方式这样做？我不禁想到答案显而易见......谢谢!!

Answer 1

尝试使用 SO 反距离加权-idw -interpolation-with-python中描述 的反距离加权和scipy.spatial.KDTree的组合 . Kd树 在2d 3d中很好地工作...,反距离加权是平滑和局部的,并且k =最近邻居的数量可以变化以折衷速度/准确度.

Answer 2

Roger Veciana i Rovira提供了一个很好的反距离示例，以及一些使用GDAL写入geotiff的代码（如果您对此很感兴趣）。

对于常规网格来说，这是粗略的，但是假设您首先使用pyproj或类似对象将数据投影到像素网格，则始终要注意对数据使用哪种投影。

他的带有test的算法的副本：

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()