在不规则网格上插值

问题描述 投票：0回答：7

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

有什么灵活的、最好是快速的方法可以做到这一点吗？ 我忍不住认为答案是显而易见的......谢谢！！

python numpy scipy interpolation
7个回答
12
投票

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

9
投票

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()
1
投票

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

你说你的输入数据来自三极数据。如何构建这些数据主要分为三种情况。

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

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

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

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

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

0
投票

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

0
投票

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

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

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

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

0
投票

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

从描述来看：

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

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

0
投票
2019年，

一种球面网格之间的替代双线性插值方法发布了一个很酷的方法，我在这里实现了四点，但使其成为通用插值函数现在已在我的待办事项列表中。 这个问题是关于数学StackExchange。

from matplotlib import pyplot as plt from matplotlib.animation import FuncAnimation import numpy as np P1 = np.array((0, 0, 1), dtype=np.float64) P2 = np.array(0(1, 0, 2), dtype=np.float64) P3 = np.array((1, 1, 2), dtype=np.float64) P4 = np.array((0, 1, 1), dtype=np.float64) points = np.column_stack((P1, P2, P3, P4)) x = np.linspace(points[0].min() - 0.2, points[0].max() + 0.2, num=400) y = np.linspace(points[1].min() - 0.2, points[1].max() + 0.2+0.5, num=400) X, Y = np.meshgrid(x, y) def irreg_grid_interp(P1, P2, P3, P4): x1, y1, f1 = P1 x2, y2, f2 = P2 x3, y3, f3 = P3 x4, y4, f4 = P4 D0 = np.array( [ [1, x1, y1, x1 * y1], [1, x2, y2, x2 * y2], [1, x3, y3, x3 * y3], [1, x4, y4, x4 * y4], ] ) D1 = np.array( [ [f1, x1, y1, x1 * y1], [f2, x2, y2, x2 * y2], [f3, x3, y3, x3 * y3], [f4, x4, y4, x4 * y4], ] ) D2 = [ [1, f1, y1, x1 * y1], [1, f2, y2, x2 * y2], [1, f3, y3, x3 * y3], [1, f4, y4, x4 * y4], ] D3 = np.array( [ [1, x1, f1, x1 * y1], [1, x2, f2, x2 * y2], [1, x3, f3, x3 * y3], [1, x4, f4, x4 * y4], ] ) D4 = np.array([[1, x1, y1, f1], [1, x2, y2, f2], [1, x3, y3, f3], [1, x4, y4, f4]]) a = np.linalg.det(D1) b = np.linalg.det(D2) c = np.linalg.det(D3) d = np.linalg.det(D4) return (a + b * X + c * Y + d * X * Y) / np.linalg.det(D0) F = irreg_grid_interp(P1, P2, P3, P4) extent = x.min(), x.max(), y.min(), y.max() im = plt.imshow(F.reshape(x.size, y.size)[::-1, :], extent=extent) pts = plt.scatter( *points[:2], c=plt.cm.viridis((points[2] - F.min()) / (F.max() - F.min())), ec="k" ) plt.box(False) plt.gca().xaxis.set_visible(False) plt.gca().yaxis.set_visible(False) theta = np.linspace(0, 2*np.pi, num=200) def update(i): t = theta[i] plt.title(rf"$\theta={np.rad2deg(t):.1f}$°") P3[0] = 0.5+np.cos(t)/2 P3[1] = 1+np.sin(t)/2 P4[0] = 0.5-np.cos(t)/2 P4[1] = 1-np.sin(t)/2 points = np.column_stack((P1, P2, P3, P4)) F[:] = irreg_grid_interp(P1, P2, P3, P4) pts.set_offsets(points[:2].T) im.set_data(F[::-1]) anim = FuncAnimation(plt.gcf(), update, frames=theta.size, interval=50) plt.gca().set_aspect('auto') # anim.save("movie.gif", dpi=100) plt.show()
最新问题
© www.soinside.com 2019 - 2024. All rights reserved.