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

Question

这是我的问题:我操纵 432*46*136*136 网格表示 time*(space) 包含在 numpy 和 python 的 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

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

谢谢

Answer 1

我不能真正尝试代码，因为我没有你的矩阵，但像这样的东西应该可以完成这项工作。

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

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

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

conv -= conv.min()

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

有了它，可以通过获取多维索引轻松完成插值:

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

上面的向量包含索引原始矩阵所需的所有索引。然后，您可以通过执行以下操作创建新矩阵:

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 包含您按高度划分的结果。