我在尝试将“几乎”规则网格数据插值到地图坐标上,以便使用
matplotlib.pyplot.imshow
同时绘制地图和数据,但是scipy.interpolate.griddata
的性能非常慢。由于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
,因此基本上我需要将它们映射到地图坐标并将相同的变换应用于数据本身。问题是如何做到这一点?