在Python中将大型不规则网格插值到另一个不规则网格上

Question

我正在尝试用Python将复杂值从一个不规则网格插值到另一个不规则网格。这个网格是二维的，总共有103,113个数据点。我使用的是Python 2.6.6，Scipy 0.7.2，Numpy 1.3.0，Matplotlib 0.99.3。

在Matlab中，使用griddata大约可以在5秒内完成这个操作。

BnGRID2  = griddata(R_GRID1,Z_GRID1,BnGRID1,R_GRID2,Z_GRID2) (MATLAB)

（注意，所有数组的大小都是201 x 513）

但是，如果我尝试使用matplotlib.mlab.griddata，就会出现内存错误，即使我只处理实部也会这样：

mlab.griddata(R_GRID1.flatten(),Z_GRID1.flatten(),num.real(BnGRID1.flatten()),R_GRID2.flatten(),Z_GRID2.flatten())

如果我尝试使用interp2d，就会出现段错误，Python会直接退出：

a = interp.interp2d(R_GRID1,Z_GRID1,num.real(BnGRID1))

我尝试使用KDTree，这似乎可以正常工作，但相比Matlab的几秒钟，这个方法需要几分钟，不过我还没有深入研究这个选项。

我在想有没有人能给我一些建议，如何能像Matlab那样快速完成这个任务？我注意到Scipy的新版本也有griddata，大家知道这个能处理大规模不规则网格吗？

Answer 1

Scipy的griddata功能似乎能够处理这么大数据集而没有问题：

import numpy as np
import scipy.interpolate

# old grid
x, y = np.mgrid[0:1:201j, 0:1:513j]
z = np.sin(x*20) * (1j + np.cos(y*3))**2   # some data

# new grid
x2, y2 = np.mgrid[0.1:0.9:201j, 0.1:0.9:513j]

# interpolate onto the new grid
z2 = scipy.interpolate.griddata((x.ravel(), y.ravel()), z.ravel(), (x2, y2), method='cubic')

在一台老旧的AMD Athlon电脑上，使用griddata大约需要5秒钟。

如果你的数据是按照网格排列的（也就是说，值z[i,j]对应的坐标是(x[i], y[j]）），那么你可以通过使用scipy.interpolate.RectBivariateSpline来提高速度。

z3 = (scipy.interpolate.RectBivariateSpline(x[:,0], y[0,:], z.real)(x2[:,0], y2[0,:])
 + 1j*scipy.interpolate.RectBivariateSpline(x[:,0], y[0,:], z.imag)(x2[:,0], y2[0,:]))

这个方法只需要0.05秒。它快得多，因为即使你的网格间距不规则，只要网格是矩形的，就可以使用更高效的算法。