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')
Plotting:
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)
Result:
Note, this can not be done by simply rotating the data with affine transformations.
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.
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?
See Question&Answers more detail:
os