title: How to accelerate scipy.map_coordinates for multiple interpolations?
body: <p>I have several values f, g, h that are defined on the same regular grid  (x, y, z) that I want to interpolate onto a new grid (x1, y1, z1). i.e., I have f(x, y, z) , g(x, y, z) , h(x, y, z) and I want to calculate f(x1, y1, z1), g(x1, y1, z1), h(x1, y1, z1).</p>
  <p>I am using <code>scipy.map_coordinates</code> at the moment. However, each interpolation is done separately and the number of points are around 4,000,000, so it is quite slow</p>
  <pre><code>from scipy.ndimage import map_coordinates
  import numpy as np
  
  # examples of f, g, h
  f=np.random.randn(100,50,50)
  g=np.random.randn(100,50,50)
  h=np.random.randn(100,50,50)
  
  # examples of x1, y1, z1
  x1=np.random.rand(4000000)*100
  y1=np.random.rand(4000000)*50
  z1=np.random.rand(4000000)*50
  
  # my solution at the moment
  coords=np.array([x1,y1,z1])
  
  out = np.zeros((3, coords.shape[1]))
  out[0]= map_coordinates(f, coords, order=1)
  out[1]= map_coordinates(g, coords, order=1)
  out[2]= map_coordinates(h, coords, order=1)
  </code></pre>
  <p>Is there a way to accelerate the calculation?</p>