Interpolation of 3D data in Python
Question: I have data in 3D. I would like to interpolate this data layer by layer (in
the plane X, Y) because calculating each layer takes a lot of time.

I try to use the interp2D function and loop through the layers but f seems to
apply only to the last value of i0.

    
    
    X = np.linspace(10., 15., 21)
    Y = np.linspace(0.05, 1.85, 10)
    for i0 in np.arange(N):
        f = interpolate.interp2d(X, Y, data[:,:,i0], kind='linear')
    
    Xnew = np.linspace(10., 15., 41)
    Ynew = np.linspace(0.05, 1.85, 20)
    datanew = f(Xnew,Ynew)
    

How can I interpolate each layer of my data?

Thanks


Answer: > I try to use the interp2D function and loop through the layers but f seems
> to apply only to the last value of i0.

You are overwriting the value of your interpolant, `f`, on each iteration of
your `for` loop, so by the time you have finished looping over `i0` values `f`
will correspond only to the last Z-plane of `data`. Using your current
approach, you would need to call `f` inside the `for` loop, e.g.:

    
    
    # some example data for test purposes
    N = 64
    data = np.random.randn(10, 21, N)
    
    X = np.linspace(10., 15., 21)
    Y = np.linspace(0.05, 1.85, 10)
    Xnew = np.linspace(10., 15., 41)
    Ynew = np.linspace(0.05, 1.85, 20)
    
    # initialize output array
    datanew = np.empty((Ynew.shape[0], Xnew.shape[0], N), data.dtype)
    
    for i0 in np.arange(N):
        f = interpolate.interp2d(X, Y, data[:,:,i0], kind='linear')
    
        # fill in this Z-slice
        datanew[:,:,i0] = f(Xnew,Ynew)
    

* * *

You could eliminate the `for` loop by interpolating over all of your Z-planes
simultaneously. One way to do this would be using
[`scipy.interpolate.RegularGridInterpolator`](http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.RegularGridInterpolator.html#scipy.interpolate.RegularGridInterpolator):

    
    
    from scipy.interpolate import RegularGridInterpolator
    
    Z = np.arange(N)
    
    itp = RegularGridInterpolator((Y, X, Z), data, method='linear')
    
    grid = np.ix_(Ynew, Xnew, Z)
    datanew2 = itp(grid)
    

Here I also use
[`np.ix_`](http://docs.scipy.org/doc/numpy-1.10.0/reference/generated/numpy.ix_.html)
to construct an 'open mesh' from the coordinates that you want to interpolate
`data` at.

