在 Python 中插值 3d 数组 展开

Question

我的问题扩展了此处看到的代码响应：Interpolating a 3d array in Python. How to avoid for loops?。相关原解决方案代码如下：

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)

当 x_new 是一个常量一维数组时，上述方法效果很好，但如果我的 x_new 不是一个常量一维数组，而是取决于纬度/经度维度的索引怎么办在另一个 3-d 数组中。我的 x_new 大小为 355x195x192（时间 x 纬度 x 长），现在我正在通过纬度和经度维度进行双重循环。由于每个纬度/经度对的 x_new 都不同，我怎样才能避免如下所示的循环？不幸的是，我的循环过程需要几个小时......

x_new=(np.argsort(np.argsort(modell, 0), 0).astype(float) + 1) / np.size(modell, 0)
## x_new is shape 355x195x192
## pobs is shape 355x1
## prism_aligned_tmax_sorted  is shape 355x195x192
interp_func = interpolate.interp1d(pobs, prism_aligned_tmax_sorted,axis=0)
tmaxmod = np.empty((355, 195, 192,))
tmaxmod[:] = np.NAN                                    
for latt in range(0, 195):
    for lonn in range(0, 192):
        temp = interp_func(x_new[:,latt,lonn])
        tmaxmod[:,latt,lonn] = temp[:,latt,lonn]

感谢您的所有帮助！

Answer 1

我知道如何摆脱这些循环，但你不会喜欢它。

问题在于interp1d 的这种使用本质上为您提供了一个在一维域上插值的矩阵值函数，即一个F(x) 函数，其中对于每个标量x，您都有一个二维数组形状的@987654324 @。您尝试做的插值是这样的：为您的每对 (lat,lon) 创建一个单独的插值器。这更像是F_(lat,lon)(x)。

这是一个问题的原因是，对于您的用例，您正在为每个查询点计算矩阵值F(x)，然后继续丢弃除单个元素之外的所有矩阵元素（元素[lat,lon] 对应于该对的查询点）。所以你正在做一堆不必要的计算来计算所有那些不相关的函数值。问题是我不确定是否有更有效的方法。

你的用例可以在你背后有适当的记忆来修复。您的循环运行数小时的事实表明这对于您的用例来说实际上是不可能的，但无论如何我会展示这个解决方案。这个想法是将您的 3d 数组转换为 2d 数组，使用此形状进行插值，然后沿插值结果的有效 2d 子空间获取对角线元素。您仍然会为每个查询点计算每个不相关的矩阵元素，但是您无需循环遍历数组，而是能够通过单个索引步骤提取相关的矩阵元素。这样做的代价是创建一个更大的辅助数组，它很可能不适合您的可用 RAM。

无论如何，这里的诀窍是，将您当前的方法与现有方法进行比较：

import numpy as np
from scipy.interpolate import interp1d
arr = np.random.randint(0, 9, size=(3, 4, 5))
x = np.linspace(0, 10, 3)
x_new = np.random.rand(6,4,5)*10

## x is shape 3 
## arr is shape  3x4x5
## x_new is shape  6x4x5

# original, loopy approach
interp_func = interp1d(x, arr, axis=0)
res = np.empty((6, 4, 5))
for lat in range(res.shape[1]):
    for lon in range(res.shape[2]):
        temp = interp_func(x_new[:,lat,lon]) # shape (6,4,5) each iteration
        res[:,lat,lon] = temp[:,lat,lon]

# new, vectorized approach
arr2 = arr.reshape(arr.shape[0],-1) # shape (3,20)
interp_func2 = interp1d(x,arr2,axis=0)
x_new2 = x_new.reshape(x_new.shape[0],-1) # shape (6,20)
temp = interp_func2(x_new2) # shape (6,20,20): 20 larger than original!
s = x_new2.shape[1] # 20, used for fancy indexing ranges
res2 = temp[:,range(s),range(s)].reshape(res.shape) # shape (6,20) -> (6,4,5)

生成的res 和res2 数组完全相同，因此该方法可能有效。但正如我所说，对于形状为 (nx,nlat,nlon) 的查询数组，我们需要一个形状为 (nx,nlat*nlon,nlat*nlon) 的辅助数组，它通常需要大量内存。

---

我能想到的唯一严格的替代方法就是一个一个地执行一维插值：在双循环中定义 nlat*nlon 插值器。这将产生更大的创建插值器的开销，但另一方面，您不会做一堆不必要的工作来计算然后丢弃的插值数组值。

最后，根据您的使用情况，您应该考虑使用多元插值（我在想interpolate.interpnd 或interpolate.griddata）。假设您的函数作为纬度和经度的函数也是平滑的，那么在更高维度上插入完整的数据集可能是有意义的。这样，您需要创建一次插值器，并准确查询您需要的点，而不会妨碍您的方式。

---

如果您最终坚持使用当前的实现，则可以通过将插值轴移动到最后一个位置来大大提高性能。这样，每个向量化操作都作用于连续的内存块（假设默认的 C 内存顺序），这非常符合“一维数组的集合”的理念。所以你应该按照以下方式做一些事情

arr = arr.transpose(1,2,0) # shape (4,5,3)
interp_func = interp1d(x, arr, axis=-1)
...
for lat ...:
    for lon ...:
        res[lat,lon,:] = temp[lat,lon,:] # shape (4,5,6)

如果需要恢复原来的顺序，最后可以用res.transpose(2,0,1)将顺序转回。