在不规则网格上进行插值

Question

所以，我有三个 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

我是否正确地认为您的数据网格看起来像这样（红色是旧数据，蓝色是新的插值数据）？

替代文本 http:// www.geekops.co.uk/photos/0000-00-02%20%28Forum%20images%29/DataSeparation.png

这可能是一种有点暴力的方法，但是渲染现有数据怎么样？作为位图（opengl 将通过配置正确的选项为您进行简单的颜色插值，并且您可以将数据渲染为三角形，这应该相当快）。然后，您可以在新点的位置对像素进行采样。

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

Answer 2

有一个名为 BIVAR 的 FORTRAN 库，它非常有用适合这个问题。通过一些修改，您可以使用 f2py 使其在 python 中可用。

从描述来看：

BIVAR 是一个 FORTRAN90 库，可对分散的双变量数据进行插值，作者：Hiroshi Akima。

BIVAR 接受一组散布在 2D 中的 (X,Y) 数据点以及相关的 Z 数据值，并且能够构建与给定数据一致的平滑插值函数 Z(X,Y)，并且可以在平面上的其他点进行评估。

Answer 3

尝试结合反距离加权和
scipy.spatial.KDTree
描述于SO
反距离-weighted-idw-interpolation-with-python。
Kd 树
在 2d 3d ... 中工作良好，反距离加权平滑且局部，
并且 k= 最近邻居的数量可以改变以权衡速度/准确性。

Answer 4

Roger Veciana 有一个 很好的反距离示例我 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 5

这里有很多选项，哪一个最好取决于您的数据......
但是我不知道适合您的现成解决方案

您说您的输入数据来自三极数据。对于如何构建这些数据，存在三种主要情况。

1. 从三极空间中的 3d 网格采样，投影回 2d LAT、LON 数据。
2. 从三极空间中的二维网格采样，投影为二维 LAT LON 数据。
3. 三极空间中的非结构化数据投影到 2d LAT LON 数据

其中最简单的是 2。不是在 LAT LON 空间中插值，而是“只是”将您的点转换回源空间并在那里插值。

适用于 1 和 2 的另一个选项是搜索从三极空间映射的单元以覆盖您的样本点。 （您可以使用 BSP 或网格类型结构来加速此搜索）选择一个单元格，然后在其中进行插值。

最后还有一堆非结构化插值选项..但它们往往很慢。
我个人最喜欢的是使用最近 N 个点的线性插值，找到这 N 个点可以再次通过网格或 BSP 来完成。另一个不错的选择是对非结构化点进行 Delauney 三角剖分，并在生成的三角网格上进行插值。

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

Answer 6

我建议您看一下 GRASS（开源 GIS 软件包）插值功能 (http://grass.ibiblio.org/gdp/html_grass62/v.surf.bspline.html）。它不是用 python 编写的，但您可以重新实现它或与 C 代码交互。