使用`scipy.interpolate.griddata`的插值非常慢 [英] Very slow interpolation using `scipy.interpolate.griddata`

问题描述

当尝试将几乎"规则网格数据插值到地图坐标时，我遇到了scipy.interpolate.griddata的极慢性能，因为matplotlib.pyplot.pcolormesh花费的时间太长并且没有用matplotlib.pyplot.imshow来绘制地图和数据与alpha一起表现良好.

I am experiencing excruciatingly slow performance of scipy.interpolate.griddata when trying to interpolate "almost" regularly gridded data to map coordinates so that both map and data can be plotted with matplotlib.pyplot.imshow because matplotlib.pyplot.pcolormesh is taking too long and does not behave well with alpha among other things.

最好显示一个示例(可以在):

Best show an example (input files can be downloaded here):

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)

结果:

注意，这不能仅通过仿射变换旋转数据来完成.

Note, this can not be done by simply rotating the data with affine transformations.

使用我的真实数据，每次调用griddata花费的时间超过80秒，而pcolormesh花费的时间甚至更长(超过2分钟！).我已经在这里中查看了这两个答案和Joe Kington的回答此处，但我想不出一种方法来使其工作我.

The griddata call takes over 80 seconds per call with my real data and pcolormesh takes even longer (over 2 minutes!). I have looked at both Jaimi's answer here and Joe Kington's answer here but I cant figure out a way to get it to work for me.

我所有的数据集都具有完全相同的lons，lats，因此基本上我需要将它们一次映射到地图的坐标，并将相同的变换应用于数据本身.问题是我该怎么办?

All my datasets have exactly the same lons, lats so basically I need to map those once to the map's coordinates and apply the same transformation to the data itself. Question is how do I do that?

推荐答案

在长时间忍受了scipy.interpolate.griddata极慢的性能后，我决定放弃griddata来支持使用 OpenCV .具体来说，透视转换.

After a long time of putting up with excruciatingly slow performance of scipy.interpolate.griddata I've decided to ditch griddata in favor of image transformation with OpenCV. Specifically, Perspective Transformation.

对于上面的示例，是上面问题中的示例，您可以为其获取输入文件

So for the above example, the one in the question above, for which you can get the input files here, this is a piece of code which takes 1.1 ms as opposed to the 692 ms which takes the regridding part in the example above.

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)

我看到的这种解决方案的唯一缺点是数据中的轻微偏移，如所附图像中不重叠的轮廓所示.以黑色和经纱重新排列的数据轮廓(可能更准确)以喷射"色标的透视数据轮廓.

The only drawback I see to this solution is a slight offset in the data as illustrated by the non-overlapping contours in the attached image. regridded data contours (probably more accurate) in black and warpPerspective data contours in 'jet' colorscale.

此刻，我在性能优势上与差异相处得很好，我希望这种解决方案也能对其他人有所帮助.

At the moment, I live just fine with the discrepancy at the advantage of performance and I hope this solution helps others as-well.

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

Someone (not me...) should find a way to improve the performance of griddata :) Enjoy!

这篇关于使用`scipy.interpolate.griddata`的插值非常慢的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！