2011-10-13 69 views
2

我有一個數組,我想插入第一個軸。目前,我正在做這樣的例子:在Python中插值3d數組。如何避免循環?

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) 

我使用的數據代表一個域中的每個緯度和經度10年5天的平均值的。我想創建一個日常值的數組。

我也試過使用樣條線。我真的不知道他們是如何工作的,但速度並不快。

有沒有辦法做到這一點,而不使用for循環? 如果必須使用for循環,還有其他方法可以加速嗎?

非常感謝您的任何建議。

回答

6

您可以指定一個軸參數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) 

+0

好點!如果OP想要真正的一維插值(而不是雙線性),那麼這就是要走的路。 –

+0

這也很好。謝謝!有趣的是(至少在這種情況下)這種方法導致內插數組的平均值更接近原始數組的平均值。 – nicholaschris

5

由於您正在插入有規律的網格數據,請查看使用scipy.ndimage.map_coordinates

作爲一個簡單的例子:

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

十分感謝,它看起來像它會工作。它會與掩碼數組一起工作嗎? – nicholaschris

+0

編輯:非常感謝,它看起來像它的效果很好。我正在使用它作爲要插入的數組的掩碼數組。這會使事情複雜化嗎?如果我設置output = array.dtype有一個奇怪的結果,但如果我離開輸出爲默認它似乎工作正常。 – nicholaschris