1

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).

Community
  • 1
  • 1
Hooked
  • 84,485
  • 43
  • 192
  • 261
  • You want the Nxdim array but tricking numpy into not actually allocating the full array? This is not possible.You would have to use a tool designed for regular grids, but this I guess does not exist for higher dimensions in scipy. – seberg Aug 16 '12 at 15:41
  • @Sebastian, you can simulate larger arrays from smaller ones. For example if `x.shape,y.shape = (n,m)`, you can create a broadcasted array of `f(x,y) = x+y` *not* by taking `[X,Y] = np.meshgrid(x,y); S=X+Y` but rather `S=x+y[:,np.newaxis]`. See the linked question for more details. – Hooked Aug 16 '12 at 15:51
  • yes, but it seems you want to simulate more elements then are actually inside the xi,zi,yi arrays. This is (with stride_tricks) in theory possible. However since the resulting array is supposed to be 2 dimensional, it is not possible to construct it for this case, and even if it is unlikely scipy would not just create a copy anyways lateron. – seberg Aug 16 '12 at 15:55
  • @Sebastian The shape of `xi,yi,zi -> (n,m,l)` dictates that I want an object `points` of shape `(n*m*l,3)`. This is exactly what `meshgrid` does for 2D `xi,yi -> (n,m)` into `(n*m,2)`. At some point scipy will create an internal copy anyways for the interpolation but I'm trying to avoid a needless copy at this step. – Hooked Aug 16 '12 at 16:01
  • @Hooked, that's not what `meshgrid` does. Your title is a little misleading. I've given an answer which does what you want as far as producing the mesh, but I don't think you can keep it this way and turn it into something of shape `(n*m*l,3)`. – user545424 Aug 16 '12 at 22:50

3 Answers3

3
>>> 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)])
user545424
  • 15,713
  • 11
  • 56
  • 70
  • 1
    The only way to avoid any extra _temporary_ arrays (that are not views) for the (m*n*l,3) points array, is with the use of xx.flat (instead of ravel or such) for insertion into the uninitialized points array. – seberg Aug 21 '12 at 20:25
  • When call `ndmesh`; what is the argument goes in ? – ElleryL Dec 26 '18 at 01:17
  • The arrays you want to broadcast. In the example I gave it would be `ndmesh(xi,yi,zi)`. – user545424 Dec 27 '18 at 16:43
1

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)

pv.
  • 33,875
  • 8
  • 55
  • 49
0

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.

seberg
  • 8,785
  • 2
  • 31
  • 30