不规则网格上的插值

Posted 2023-02-16

【中文标题】不规则网格上的插值

【英文标题】：Interpolation over an irregular grid

【发布时间】：2011-03-15 14:49:16

【问题描述】：

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

有没有一种巧妙的，最好是快速的方法来做到这一点？我不禁认为答案很明显......谢谢！

【问题讨论】：

标题中的“不规则网格”让我有点失望。您有一个恰好分布在空间中的点样本，但您没有matplotlib.org/examples/pylab_examples/tripcolor_demo.html 中的网格结构您的数据是分散在一个字段中的点，您可以假设它有些平滑。可以使用 matplotlib.tri matplotlib.org/api/tri_api.html 对不规则或非结构化网格或网格进行插值，以考虑场中的不连续性。

【参考方案1】：

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

【讨论】：

你，我的朋友，是个天才。 KDTree 课程很棒！正是我需要的......

我在使用 vanilla 逆权重时遇到了一些麻烦。发现当样本点在点簇之外时，它有一些严重的伪影。我通过将线性近似（而不是常数近似）拟合到 N 最近邻的加权数据来解决这个问题。这在相同的搜索量下产生了相当不错的结果，只是求解 NxN 线性系统的开销。

@Michael，你的数据是 2d 的，有多分散，Nnear 是什么？你能举一个行为不端的距离和价值观的例子吗？例如距离 1 1 1 1 1 10，值 1 1 1 1 1 10 => 插值 (6 / 5.1) = 1.18。另外，NxN ?在 2d 中，将平面 ax + by + c 拟合到 N 个点（权重为 1/dist）是 numpy.linalg .lstsq Nx3 或 .solve 3x3 。

我的数据是 3D，但即使是 1D 也会出现问题。使用线性数据 (1,1) (2,2),(3,3) 取 N=3，在 2.5 处采样，您会得到大约 2.3 的结果（低估 10%）。如果我们估计为 3.5，情况会更糟，产生接近 2.5 而不是“真实” 3.5 的值。有人会说我们现在是在做外推而不是插值，但如果我们的数据点在 1,2,3,10 那么 1,2,3 仍然是最接近 3.5 的三个点......我们最终会得到同样的结果。这就是我所说的集群之外的值。拟合一条线给出了“正确”的结果——至少对于我的数据而言

【参考方案2】：

如果您对此感兴趣，有一个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() 

【参考方案3】：

这里有很多选项，哪个最好取决于您的数据... 但是我不知道为您提供开箱即用的解决方案

您说您的输入数据来自三极数据。该数据的结构主要有以下三种情况。

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

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

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

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

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

【讨论】：

+1: ..对于 BSP 树的提及，通常把我得到的东西比我管理的更重要:-) 你可以通过将每个 BSP 节点集中在其中一个上来形成 BSP新数据点，然后简单地向下钻取以找到所有相邻点。

不错！共识似乎是我将不得不在这方面工作一点，但这没关系。我很喜欢你对 BSP 技术的建议......非常感谢！

案例 3 的一部分可能是您在非结构化网格上定义了一个数据，其中生成的 Delauney 凸包可能不合适。例如。 matplotlib.org/examples/pylab_examples/tripcolor_demo.html 然后在给定的三角形网格上插值可能会很好：matplotlib.org/api/tri_api.html

【参考方案4】：

我建议您看看 GRASS（一个开源 GIS 包）插值功能 (http://grass.ibiblio.org/gdp/html_grass62/v.surf.bspline.html)。它不在 python 中，但您可以重新实现它或与 C 代码接口。

【讨论】：

嗯，这当然看起来不错，虽然需要重新实现一些工作！我会调查的。谢谢！

无需重新实现，调用即可。使用 SEXTANTE 工具箱查看 QGIS。

【参考方案5】：

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

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

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

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

【讨论】：

网格的想法是正确的，虽然我实际上是在跟踪虚拟粒子穿过网格时的属性，所以蓝点看起来更像是面包屑的痕迹：!mesh希望那张照片有效。图像渲染的想法很有趣——我确实有可用的 PIL，所以我可以试一试。谢谢！

【参考方案6】：

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

来自描述：

BIVAR 是一个 FORTRAN90 库，它对分散的二元数据进行插值，作者是 Hiroshi Akima。

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