Interpolating a 3d array in Python. How to avoid for loops?

Question

I have an array which I want to interpolate over the 1st axes. At the moment I am doing it like this example:

import numpy as np
from scipy.interpolate import interp1d

array = np.random.randint(0, 9, size=(100, 100, 100))
new_array = np.zeros((1000, 100, 100))
x = np.arange(0, 100, 1)
x_new = np.arange(0, 100, 0.1)

for i in x:
    for j in x:
        f = interp1d(x, array[:, i, j])
        new_array[:, i, j] = f(xnew)

The data I use represents 10 years of 5-day averaged values for each latitude and longitude in a domain. I want to create an array of daily values.

I have also tried using splines. I don't really know how they work but it was not much faster.

Is there a way to do this without using for loops? If for loops must be used, are there other ways to speed this up?

Thank you in advance for any suggestions.

pv. · Accepted Answer · 2011-10-15 20:40:55Z

7

You can specify an axis argument to interp1d:

import numpy as np
from scipy.interpolate import interp1d
array = np.random.randint(0, 9, size=(100, 100, 100))
x = np.linspace(0, 100, 100)
x_new = np.linspace(0, 100, 1000)
new_array = interp1d(x, array, axis=0)(x_new)
new_array.shape # -> (1000, 100, 100)

answered Oct 15, 2011 at 20:40

pv.

35.4k9 gold badges62 silver badges51 bronze badges

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

2 Comments

Joe Kington Over a year ago

Good point! If the OP wants truly 1D interpolation (and not bilinear), then this is the way to go.

nicholaschris Over a year ago

This also works well. Thank you! Interesting to note that (in this case at least) this method results in the mean of the interpolated array being closer to the mean of the original array.

Joe Kington · Accepted Answer · 2011-10-13 15:02:47Z

6

Because you're interpolating regularly-gridded data, have a look at using scipy.ndimage.map_coordinates.

As a quick example:

import numpy as np
import scipy.ndimage as ndimage

interp_factor = 10
nx, ny, nz = 100, 100, 100
array = np.random.randint(0, 9, size=(nx, ny, nz))

# If you're not familiar with mgrid: 
# http://docs.scipy.org/doc/numpy/reference/generated/numpy.mgrid.html
new_indicies = np.mgrid[0:nx:interp_factor*nx*1j, 0:ny, 0:nz]

# order=1 indicates bilinear interpolation. Default is 3 (cubic interpolation)
# We're also indicating the output array's dtype should be the same as the 
# original array's. Otherwise, a new float array would be created.
interp_array = ndimage.map_coordinates(array, new_indicies, 
                                       order=1, output=array.dtype)
interp_array = interp_array.reshape((interp_factor * nx, ny, nz))

answered Oct 13, 2011 at 15:02

Joe Kington

287k73 gold badges621 silver badges474 bronze badges

2 Comments

nicholaschris Over a year ago

Great thanks, it looks like it will work. Will it work with masked arrays?

nicholaschris Over a year ago

EDIT: Great thanks, it looks like it works well. I am using it with a masked array as the array to be interpolated. Will that complicate matters? If I set the output= array.dtype there is a strange result but if I leave the output as default it seems to work fine.

Collectives™ on Stack Overflow

Interpolating a 3d array in Python. How to avoid for loops?

2 Answers 2

2 Comments

2 Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

2 Comments

2 Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related