开发者

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

开发者 https://www.devze.com 2023-04-12 15:27 出处:网络
I have an array which I want to interpolate over the 1st axes. At the moment I am doing it like this example:

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.interpol开发者_如何学Pythonate 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.


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)


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))
0

精彩评论

暂无评论...
验证码 换一张
取 消

关注公众号