title: Very slow interpolation using `scipy.interpolate.griddata`
body: <p>I am experiencing excruciatingly slow performance of <code>scipy.interpolate.griddata</code> when trying to interpolate "almost" regularly gridded data to map coordinates so that both map and data can be plotted with <code>matplotlib.pyplot.imshow</code> because <code>matplotlib.pyplot.pcolormesh</code> is taking too long and does not behave well with <code>alpha</code> among other things.</p>
  
  <p>Best show an example (input files can be downloaded <a href="https://www.dropbox.com/s/m215o0ko304173d/input.zip?dl=0" rel="noreferrer">here</a>):</p>
  
  <pre><code>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')
  </code></pre>
  
  <p>Plotting:</p>
  
  <pre><code>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)
  </code></pre>
  
  <p>Result:</p>
  
  <p><img src="https://i.stack.imgur.com/lh4FX.png" alt="enter image description here"></p>
  
  <p><strong>Note</strong>, this can not be done by simply rotating the data with affine transformations.</p>
  
  <p>The <code>griddata</code> call takes over 80 seconds per call with my real data and <code>pcolormesh</code> takes even longer (over 2 minutes!). I have looked at both Jaimi's answer <a href="https://stackoverflow.com/questions/20915502/speedup-scipy-griddata-for-multiple-interpolations-between-two-irregular-grids">here</a> and Joe Kington's answer <a href="https://stackoverflow.com/questions/16983843/fast-interpolation-of-grid-data">here</a> but I cant figure out a way to get it to work for me.</p>
  
  <p>All my datasets have exactly the same <code>lons</code>, <code>lats</code> 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? </p>