使用 `scipy.interpolate.griddata` 进行非常慢的插值

Question

我scipy.interpolate.griddata在尝试将“几乎”有规律的网格数据插入到地图坐标中时遇到了极其缓慢的性能，以便可以绘制地图和数据，matplotlib.pyplot.imshow因为matplotlib.pyplot.pcolormesh它花费的时间太长并且alpha在其他方面表现不佳。

最好展示一个例子（输入文件可以在这里下载）：

import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import griddata

map_extent = (34.4, 36.2, 30.6, 33.4)
# data corners:
lon = np.array([[34.5,        34.83806236],
                [35.74547079, 36.1173923]])
lat = np.array([[30.8,        33.29936152],
                [30.67890411, 33.17826563]])

# load saved files
topo = np.load('topo.npy')
lons = np.load('lons.npy')
lats = np.load('lats.npy')
data = np.load('data.npy')

# get max res of data
dlon = abs(np.array(np.gradient(lons))).max()
dlat = abs(np.array(np.gradient(lats))).max()

# interpolate the data to the extent of the map
loni,lati = np.meshgrid(np.arange(map_extent[0], map_extent[1]+dlon, dlon),
                        np.arange(map_extent[2], map_extent[3]+dlat, dlat))
zi = griddata((lons.flatten(),lats.flatten()),
              data.flatten(), (loni,lati), method='linear')

绘图：

fig, (ax1,ax2) = plt.subplots(1,2)
ax1.axis(map_extent)
ax1.imshow(topo,extent=extent,cmap='Greys')

ax2.axis(map_extent)
ax2.imshow(topo,extent=extent,cmap='Greys')

ax1.imshow(zi, vmax=0.1, extent=extent, alpha=0.5, origin='lower')
ax1.plot(lon[0],lat[0], '--k', lw=3, zorder=10)
ax1.plot(lon[-1],lat[-1], '--k', lw=3, zorder=10)
ax1.plot(lon.T[0],lat.T[0], '--k', lw=3, zorder=10)
ax1.plot(lon.T[-1],lat.T[-1], '--k', lw=3, zorder=10)


ax2.pcolormesh(lons,lats,data, alpha=0.5)
ax2.plot(lon[0],lat[0], '--k', lw=3, zorder=10)
ax2.plot(lon[-1],lat[-1], '--k', lw=3, zorder=10)
ax2.plot(lon.T[0],lat.T[0], '--k', lw=3, zorder=10)
ax2.plot(lon.T[-1],lat.T[-1], '--k', lw=3, zorder=10)

结果：

请注意，这不能通过简单地使用仿射变换旋转数据来完成。

该griddata调用将每次通话超过80秒。我真实的数据和pcolormesh需要更长的时间（超过2分钟！）。我在这里查看了 Jaimi 的回答和 Joe Kington 的回答，但我无法找到让它为我工作的方法。

我所有的数据集都完全相同lons，lats所以基本上我需要将这些映射一次到地图的坐标，并对数据本身应用相同的转换。问题是我该怎么做？

Answer 1

在长时间忍受极其缓慢的性能之后，scipy.interpolate.griddata我决定放弃griddata使用OpenCV进行图像转换。具体来说，透视变换。

所以对于上面的例子，上面问题中的那个，你可以在这里获取输入文件，这是一段需要 1.1 ms 的代码，而不是 692 ms 在上面的例子中需要重新网格部分。

import cv2
new_data = data.T[::-1]

# calculate the pixel coordinates of the
# computational domain corners in the data array
w,e,s,n = map_extent
dx = float(e-w)/new_data.shape[1]
dy = float(n-s)/new_data.shape[0]
x = (lon.ravel()-w)/dx
y = (n-lat.ravel())/dy

computational_domain_corners = np.float32(zip(x,y))

data_array_corners = np.float32([[0,new_data.shape[0]],
                                 [0,0],
                                 [new_data.shape[1],new_data.shape[0]],
                                 [new_data.shape[1],0]])

# Compute the transformation matrix which places
# the corners of the data array at the corners of
# the computational domain in data array pixel coordinates
tranformation_matrix = cv2.getPerspectiveTransform(data_array_corners,
                                                   computational_domain_corners)

# Make the transformation making the final array the same shape
# as the data array, cubic interpolate the data placing NaN's
# outside the new array geometry
mapped_data = cv2.warpPerspective(new_data,tranformation_matrix,
                                  (new_data.shape[1],new_data.shape[0]),
                                  flags=2,
                                  borderMode=0,
                                  borderValue=np.nan)

我看到此解决方案的唯一缺点是数据中的轻微偏移，如附加图像中的非重叠轮廓所示。黑色的重新网格数据轮廓（可能更准确）和“jet”色标中的 warpPerspective 数据轮廓。

目前，我对性能优势的差异感到满意，我希望这个解决方案也能帮助其他人。

有人（不是我......）应该找到一种方法来提高 griddata 的性能 :) 享受吧！