numpy 3D meshgrid only as a view

Question

I have a cubic grid defined by the spacing xi,yi,zi:

xi,yi,zi = [linspace(ox,ox+s*d,s) for ox,s,d in zip(origin,size,delta)]

I also have set of scalar values W onto that grid. W.shape() == size. I'd like to use scipy's linear interpolation, which requires as input:

class scipy.interpolate.LinearNDInterpolator(points, values):

Parameters :

points : ndarray of floats, shape (npoints, ndims) Data point coordinates.

values : ndarray of float or complex, shape (npoints, ...) Data values.

How do I create a fake set of points (via magical broadcasting) from xi,yi,zi? Right now I'm creating an intermediate array to feed to the interpolation function - is there a better way?

Related Question: Numpy meshgrid in 3D. The answers in this post actually create the grid - I only want to simulate it as input to another function (pure numpy solution preferred).

Answer 1

>>> xi, yi, zi = [np.arange(3) for i in range(3)]
>>> xx, yy, zz = np.broadcast_arrays(xi,yi[:,np.newaxis],zi[:,np.newaxis,np.newaxis])
>>> xx.shape
(3, 3, 3)
>>> xx.strides
(0, 0, 8)

You can see it didn't create new copies since the strides are 0 in the first two dimensions.

I wrote a n dimensional version of this also:

def ndmesh(*args):
   args = map(np.asarray,args)
   return np.broadcast_arrays(*[x[(slice(None),)+(None,)*i] for i, x in enumerate(args)])

Answer 2

You can construct the necessary points array in a similar way as explained in the other answers:

xx, yy, zz = np.broadcast_arrays(xi[:,None,None], yi[None,:,None], zi[None,None,:])
points = (xx.ravel(), yy.ravel(), zz.ravel())
ip = LinearNDInterpolator(points, data.ravel())

However, if you have a regular grid, then using LinearNDInterpolator is most likely not the best choice, since it is designed for scattered data interpolation. It constructs a Delaunay triangulation of the data points, but in this case the original data has already a very regular structure that would be more efficient to make use of.

Since your grid is rectangular, you can build up the interpolation as a tensor product of three 1-D interpolations. Scipy doesn't have this built-in (so far), but it's fairly easy to do, see this thread: http://mail.scipy.org/pipermail/scipy-user/2012-June/032314.html (use e.g. interp1d instead of pchip to get 1-D interpolation)

Answer 3

I do not believe there is any way you can pass something to LinearNDInterpolator short of a full copy (as there are no functions for regular grids in three dimensions too). So the only place to avoid creating full arrays would be during creation of this points array, I do not know how you do it right now, so maybe it is already efficient in this regard, but I guess its likely not worth the trouble to avoid this.

Other then np.mgrid+reshape maybe something like this might be an option (not to hard to write for n-dimensions too):

# Create broadcastest versions of xi, yi and zi
# np.broadcast_arrays does not allocate the full arrays
xi, yi, zi = np.broadcast_arrays(xi[:,None,None], yi[:,None,None], zi[:,None,None])

# then you could use .flat to fill a point array:
points = np.empty((xi.size, 3), dtype=xi.dtype)
points[:,0] = xi.flat
points[:,1] = yi.flat
points[:,2] = zi.flat

Opposed to the .repeat function, the temporary arrays created here are not larger then the original xi, etc. arrays.