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 1

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:

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_ to construct an 'open mesh' from the coordinates that you want to interpolate data at.