python - 在 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 不是一个常量一维数组，而是取决于纬度/经度维度的索引怎么办在另一个三维数组中。我的 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 你有一个二维数组形状的 F。您尝试进行的插值是这样的:为您的每个 (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) 的辅助数组，这通常需要大量内存。

---

我能想到的唯一严格的替代方法是逐一执行 1d 插值:在双循环中定义 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) 调换顺序。