CODEWITHSUNDEEP
pythonnumpyscipyinterpolation
Fred
Fred

Reputation: 95

Interpolation of 3D data in Python

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

Upvotes: 4

Views: 7167

Answers (2)

ali_m
ali_m

Reputation: 74154

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.

Upvotes: 5

xnx
xnx

Reputation: 25478

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)

Upvotes: 2

Related Questions