4

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

Improve this question
2
  • 1
    what is you value for N? Commented Jan 4, 2016 at 12:12
  • The value for N is 64 Commented Jan 4, 2016 at 12:48

2 Answers 2

Reset to default
5

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.

Improve this answer
Sign up to request clarification or add additional context in comments.

Comments

2

At the moment you are creating the interpolation function for each layer and then throwing it away before creating it for the next one. That's why only the last one is still definined when you interpolate onto your new grid.

If i0 is indexing each "layer" of your data and you're interpolating within each layer only, surely you need to create your interpolation function, f, for each layer and update accordingly each time round the loop. Something like the following (note you may need to rollaxis on datanew to get them back in the order you want):

import numpy as np
from scipy import interpolate

N = 64
data = np.random.random((21, 10, 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)
datanew = np.empty((N, 41, 20))
for i0 in np.arange(N):
    f = interpolate.interp2d(X, Y, data[:,:,i0].T, kind='linear')
    datanew[i0] = f(Ynew,Xnew)
datanew = np.rollaxis(datanew, 0, 3)
Improve this answer

2 Comments

Thanks. The only problem with this is that I need to do the interpolation everytime I want to change Xnew and Ynew. For my current case this is not a big deal because my data is relatively small.
You could store a list of interpolation objects and index them f[i0](Ynew,Xnew).

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.