matplotlib.mlab.griddata非常慢，输入有效数据却返回nan数组

Question

我正在尝试将一个不规则网格的数据集（原始卫星数据）与一些经纬度对应到一个规则网格的经纬度，这个规则网格是通过 basemap.makegrid() 生成的。我使用了 matplotlib.mlab.griddata，并且安装了 mpl_toolkits.natgrid。下面是我在 ipython 中用 whos 命令查看到的一些变量和它们的统计信息：

Variable   Type       Data/Info
-------------------------------
datalat    ndarray    666x1081: 719946 elems, type `float32`, 2879784 bytes (2 Mb)
datalon    ndarray    666x1081: 719946 elems, type `float32`, 2879784 bytes (2 Mb)
gridlat    ndarray    1200x1000: 1200000 elems, type `float64`, 9600000 bytes (9 Mb)
gridlon    ndarray    1200x1000: 1200000 elems, type `float64`, 9600000 bytes (9 Mb)
var        ndarray    666x1081: 719946 elems, type `float32`, 2879784 bytes (2 Mb)

In [11]: var.min()
Out[11]: -30.0

In [12]: var.max()
Out[12]: 30.0

In [13]: datalat.min()
Out[13]: 27.339874

In [14]: datalat.max()
Out[14]: 47.05302

In [15]: datalon.min()
Out[15]: -137.55658

In [16]: datalon.max()
Out[16]: -108.41629

In [17]: gridlat.min()
Out[17]: 30.394031556984299

In [18]: gridlat.max()
Out[18]: 44.237140350357713

In [19]: gridlon.min()
Out[19]: -136.17646180595321

In [20]: gridlon.max()
Out[20]: -113.82353819404671

datalat 和 datalon 是原始数据的坐标。

gridlat 和 gridlon 是我想要插值到的坐标。

var 包含了实际的数据。

使用这些变量，当我调用 griddata(datalon, datalat, var, gridlon, gridlat) 时，处理时间最长可以达到20分钟，并且返回的结果是一个包含 nan 的数组。从数据来看，经纬度似乎是正确的，原始坐标与新区域有部分重叠，还有一些数据点位于新区域之外。有没有人有什么建议？这些 nan 值让我觉得我可能做错了什么……

Answer 1

如果你的数据是以网格的形式排列的，比如在点 (datalon[i], datalat[j]) 这个位置的数据就是 data[i,j]，那么你可以使用 scipy.interpolate.RectBivariateSpline 这个工具，而不是用 griddata。不过，有些专门针对地理数据的库可能会提供更多的功能。

Answer 2

很可能，griddata这个东西太复杂了。它是为了处理随机采样的数据而设计的。而你的数据几乎肯定是规则采样的，只是它的网格和你想要的输出网格不一样。

你可以考虑一种更简单的方法，比如使用仿射变换，或者在小块区域上进行一系列的仿射变换，特别是如果地球的地形或曲率会影响你的结果的话。

还有一些现成的解决方案可能会对你有帮助，比如GDAL就是一个很好的例子。

此外，这类问题在地理信息系统（GIS）中也经常被讨论。你可以看看这个链接：

https://gis.stackexchange.com/questions/10430/changing-image-projection-using-python

Answer 3

看起来，mlab.griddata这个功能可能会对你的输出数据施加一些额外的限制，而这些限制可能并不是必要的。虽然输入的位置可以是任何地方，但输出的位置必须是规则的网格。因为你的例子是在经纬度空间中，所以你选择的地图投影可能会导致这个规则被违反（也就是说，在x/y坐标下的规则网格，在经纬度下可能就不是规则网格了）。

你可以尝试使用来自SciPy的interpolate.griddata作为替代方案。不过，你需要把你的位置信息合并成一个单独的数组，因为这个函数的调用方式不同：大概是这样的

import scipy.interpolate
data_locations = np.vstack(datalon.ravel(), datalat.ravel()).T
grid_locations = np.vstack(gridlon.ravel(), gridlat.ravel()).T
grid_data      = scipy.interpolate.griddata(data_locations, val.ravel(),
                                            grid_locations, method='nearest')

用于最近邻插值。这会把位置放进一个有两列的数组，分别对应你的两个维度。你可能还想在你地图投影的变换空间中进行插值。