在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是一个常数1-d阵列，但如果我的x_new 也不是什么恒定的1-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给你基本上为x你有一个2D阵列状F每个标量内插在1D域的矩阵值函数，即一个F(x)函数，其中。你要做的插值是这样的：为每一对(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内存顺序），并且这符合“1d数组的集合”哲学。所以，你应该做的沿

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)调换顺序回来。