避免 for 循环，在 3D numpy 数组的维度上进行插值

Question

我有一个大的 3D numpy 维度数组 (l, n, m)，其元素对应于 x、y 和 z 的 1D 数组，维度分别为 l、n 和 m。我想通过在 x 和 y 的每个组合的 z 值之间进行插值来找到给定 a 值（长度为 b）的元素。这将给出一个尺寸为 (l, n, b) 的输出 3D 数组。我想完全使用 numpy 数组而不是诉诸 for 循环。

例如，如果我的 3D 数组有维度 (2,3,4)：

x = 1 | z = 1 | 2 | 3 | 4 
- - - - - - - - - - - - - - 
y = 1 |[[[ 0,  1,  2,  3],       
y = 2 |  [ 4,  5,  6,  7],
y = 3 |  [ 8,  9, 10, 11]],

x = 2 | z = 1 | 2 | 3 | 4 
- - - - - - - - - - - - -
y = 1 | [[ 12, 13, 14, 15],        
y = 2 |  [ 16, 17, 18, 19],
y = 3 |  [ 20, 21, 22, 23]]]

我想对每一行进行插值 {(x=1, y=1), (x=1, y=2), (x=1, y=3), (x=2, y=1 ), (x=2, y=2), (x=2, y=3)} 对于 a = [1.3, 1.8, 2.34, 2.9, 3.45] 的值给出一个维度为 (2,3, 5):

[[[  0.3,  0.8,  1.34,  1.9,  2.45],
  [  4.3,  4.8,  5.34,  5.9,  6.45],
  [  8.3,  8.8,  9.34,  9.9, 10.45]],

 [[ 12.3, 12.8, 13.34, 13.9, 14.45],
  [ 16.3, 16.8, 17.34, 17.9, 18.45],
  [ 20.3, 20.8, 21.34, 21.9, 22.45]]]

目前我使用 for 循环遍历 x 和 y 的每个组合，并将我的 3D 数组的行输入到 numpy.iterpolate 函数中，并将输出保存到另一个数组中；然而，这对于大型数组来说非常慢。

# array is the 3D array with dimensions (l, n, m)
# x, y and z have length l, n and m respectively
# a is the values at which I wish to interpolate at with length b
# new_array is set up with dimensions (l, n, b) 

new_array = N.zeros(len(x)*len(y)*len(a)).reshape(len(x), len(y), len(a))
for i in range(len(x)):
      for j in range(len(y)):
               new_array[i,j,:] = numpy.interpolate(a, z, array[i,j,:])

任何帮助将不胜感激。

Answer 1

您不需要 for 循环即可通过 scipy.interpolate.griddata 运行数据：

>>> from itertools import product
>>>from scipy.interpolate import griddata

>>> data = np.arange(24).reshape(2, 3, 4)

>>> x = np.arange(1, 3)
>>> y = np.arange(1, 4)
>>> z = np.arange(1, 5)
>>> points = np.array(list(product(x, y, z)))

# This is needed if your x, y and z are not consecutive ints
>>> _, x_idx = np.unique(x, return_inverse=True)
>>> _, y_idx = np.unique(y, return_inverse=True)
>>> _, z_idx = np.unique(z, return_inverse=True)
>>> point_idx = np.array(list(product(x_idx, y_idx, z_idx)))
>>> values = data[point_idx[:, 0], point_idx[:, 1], point_idx[:, 2]]

>>> new_z = np.array( [1.3, 1.8, 2.34, 2.9, 3.45])
>>> new_points = np.array(list(product(x, y, new_z)))
>>> new_values = griddata(points, values, new_points)
>>> new_values.reshape(2, 3, -1)
array([[[  0.3 ,   0.8 ,   1.34,   1.9 ,   2.45],
        [  4.3 ,   4.8 ,   5.34,   5.9 ,   6.45],
        [  8.3 ,   8.8 ,   9.34,   9.9 ,  10.45]],

       [[ 12.3 ,  12.8 ,  13.34,  13.9 ,  14.45],
        [ 16.3 ,  16.8 ,  17.34,  17.9 ,  18.45],
        [ 20.3 ,  20.8 ,  21.34,  21.9 ,  22.45]]])