插值过不规则栅格

Question

所以，我有存储纬度，经度，和在网格上的一些属性值3个numpy的阵列 - 也就是说，我有LAT（Y，X），LON（Y，X），和，说温度T（y，x），对于x和y的一些限制。网格不一定是规则的 - 事实上，它是三极的。

然后，我想要将这些属性（温度）值插值到一堆不同的纬度/经度点（存储为lat1（t），lon1（t），约10,000 t ...）实际的网格点。我试过matplotlib.mlab.griddata，但这需要太长时间（毕竟，它并不是真正为我所做的而设计的）。我也尝试scipy.interpolate.interp2d，但我得到一个MemoryError（我的网格大约400x400）。

是否有任何形式的光滑，最好是快速的方式呢？我不禁想到答案是显而易见的......谢谢！

Answer 1

尝试在SO inverse-distance-weighted-idw-interpolation-with-python描述逆距离加权的组合和 scipy.spatial.KDTree 。 Kd-trees 在二维3d ...中很好地工作，反距离加权是平滑和局部的，并且k =最近邻居的数量可以变化以权衡速度/准确度。

Answer 2

我建议你考虑看看草（一种开源GIS包）插补功能（http://grass.ibiblio.org/gdp/html_grass62/v.surf.bspline.html）。它不是在Python中，但你可以重新实现它或与C代码接口。

Answer 3

我是正确的思维你的数据网格是这个样子（红色是旧数据，蓝色是新的内插数据）？

alt text http://www.geekops.co.uk/photos/0000-00-02%20%28Forum%20images%29/DataSeparation.png

这可能是一个稍微强力十岁上下的方法，但有关呈现您的现有数据作为位图（OpenGL的会做的颜色简单的插值为您配置和正确的选项，你可以渲染什么数据应该是相当快的三角形）。然后，您可以在新点的位置采样像素。

或者，你可以在空间上的排序第一组点，然后找到最接近的旧点周围的新点，并根据距离对这些点进行插值。

Answer 4

有一堆的选择这里，哪一个是最好的取决于你的数据... 但是我不知道出的现成的解决方案，你

你说你的输入数据来自三极数据。关于如何构建这些数据，主要有三种情况。

1. 从三维空间中的三维网格采样，投影回第2d个LAT，LON数据。
2. 从三维空间中的2d网格采样，投影到2d LAT LON数据。投影到二维LAT LON数据

最简单的，这些在三极空间

非结构化数据是2而不是插在LAT LON空间“只是”改变你点回源空间和内插在那里。

，对于1和2的工作原理是搜索，从三极空间映射到覆盖你的采样点的细胞的另一种选择。 （您可以使用BSP或网格类型结构来加速此搜索）选取其中一个单元格，然后在其中进行插值。

最后还有非结构化插选项的堆..但他们往往是缓慢的。 我个人最喜欢的是使用最近N个点的线性插值，发现这些N个点可以再次用网格或BSP来完成。另一个不错的选择是Delauney对非结构化点进行三角测量，并在生成的三角形网格上进行插值。

就个人而言，如果我的网格是案例1，我会使用非结构化策略，因为我担心必须处理通过具有重叠投影的单元格进行搜索。选择“正确”的单元格会很困难。

Answer 5

有一个nice inverse distance example by Roger Veciana i Rovira以及一些使用GDAL编写geotiff的代码，如果你进入。

这对于规则网格来说很粗糙，但是假设您将数据首先投影到具有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 6

有一个FORTRAN库调用BIVAR，这是非常适合这个问题。通过一些修改，你可以使用f2py在python中使用它。

从描述：

BIVAR是FORTRAN90库，内插散射二元数据，由Hiroshi阿克玛。 BIVAR接受一组散布在2D中的（X，Y）数据点以及相关的Z数据值，并且能够构造与给定数据一致的平滑插值函数Z（X，Y），并且可以在飞机上的其他点进行评估。