使用通用函数的 Python numpy 网格转换

Posted 2023-02-23

【中文标题】使用通用函数的 Python numpy 网格转换

【英文标题】：Python numpy grid transformation using universal functions

【发布时间】：2015-09-08 13:15:25

【问题描述】：

这是我的问题：我使用 numpy 和 python 操作代表 time*(space) 的 432*46*136*136 网格包含在 numpy 数组中。我有一个数组alt，它包含网格点的高度，另一个数组temp 存储网格点的温度。

比较有问题：如果T1和T2是两个结果，T1[t0,z0,x0,y0]和T2[t0,z0,x0,y0]分别代表H1[t0,z0,x0,y0]和H2[t0,z0,x0,y0]米处的温度。但我想比较的是同一高度点的温度，而不是同一网格点。

因此我想修改矩阵的 z 轴来表示高度而不是网格点。我创建了一个函数conv(alt[t,z,x,y])，它为每个高度赋予了一个介于 -20 和 200 之间的数字。这是我的代码：

def interpolation_extended(self,temp,alt):
    [t,z,x,y]=temp.shape
    new=np.zeros([t,220,x,y])
    for l in range(0,t):
       for j in range(0,z):
          for lat in range(0,x):
             for lon in range(0,y):
                new[l,conv(alt[l,j,lat,lon]),lat,lon]=temp[l,j,lat,lon]
    return new

但这肯定需要太多时间，我不能这样做。我尝试使用带有 numpy 的通用函数来编写它：

def interpolation_extended(self,temp,alt):
    [t,z,x,y]=temp.shape
    new=np.zeros([t,220,x,y])
    for j in range(0,z):
       new[:,conv(alt[:,j,:,:]),:,:]=temp[:,j,:,:]
    return new

但这不起作用。您是否知道在 python/numpy 中不使用 4 个嵌套循环来执行此操作？

谢谢

【问题讨论】：

这可能是因为 Python 中的函数调用成本很高。预先计算 conv 坐标，例如fromfunction numpy 方法可能会有所帮助。

z 是否大于 220？如果是这种情况，我可以看到您的代码为什么会崩溃。

conv 是alt 的任何线性组合/函数吗？你能举个例子或发布你的代码conv吗？如果conv可以向量化，可以写成1行。

感谢您的回答。不，z 介于 0 和 46 之间。我的海拔在 -10.000m 和 +100.000m 之间，所以 conv 实际上只是：round(alt[t,z,x,y]/500.)（我弄错了，conv 取值在-20 和 200，并且不在 0 和 220 之间）

【参考方案1】：

我无法真正尝试代码，因为我没有你的矩阵，但这样的事情应该可以完成。

首先，不要将conv 声明为函数，而是获取所有数据的整个高度投影：

conv = np.round(alt / 500.).astype(int)

使用np.round，圆的numpys版本，它通过向量化C中的操作对矩阵的所有元素进行四舍五入，因此，您可以非常快速地获得一个新数组（以C速度）。通过将所有数组移动其最小值（在您的情况下为-20），以下行将高度对齐以从 0 开始：

conv -= conv.min()

上面的行会将您的高度矩阵从 [-20, 200] 转换为 [0, 220]（更适合索引）。

这样，可以通过获取多维索引轻松完成插值：

t, z, y, x = np.indices(temp.shape)

上面的 vectors 包含索引原始矩阵所需的所有索引。然后，您可以通过以下方式创建新矩阵：

new_matrix[t, conv[t, z, y, x], y, x] = temp[t, z, y, x]

完全没有任何循环。

让我知道它是否有效。它可能会给你一些错误，因为我很难在没有数据的情况下对其进行测试，但它应该可以完成这项工作。

---

以下玩具示例运行良好：

A = np.random.randn(3,4,5) # Random 3x4x5 matrix -- your temp matrix
B = np.random.randint(0, 10, 3*4*5).reshape(3,4,5) # your conv matrix with altitudes from 0 to 9
C = np.zeros((3,10,5)) # your new matrix

z, y, x = np.indices(A.shape)
C[z, B[z, y, x], x] = A[z, y, x]

C 包含您的海拔结果。

【讨论】：

效果很好，非常感谢。我什至不知道将索引用作向量的可能性。我不得不改变：conv = np.round(alt / 500.) 到 conv = np.round(alt / 500.).astype(int) 使用 conv 作为索引。

@Arnaud PROST 谢谢！已经编辑好了。矢量索引非常有用。例如：A[[3, 7, 8], :]（对于A 2D 数组）将返回矩阵的第 4 行、第 8 行和第 9 行。使用np.indices 的答案是对多维数组中所有轴的所有元素的索引的扩展。