从numpy中另一个多维数组中的二维索引数组创建多维数组

Question

我遇到了一个我正在努力解决的问题。

我有一个形状为 (3,4,5) 的数组，其中包含 3 个时间步长的值，用于 4 个指标和 5 个区域。为了这个例子的目的，假设这个数组是：

values = np.arange(60).reshape(3,4,5)

我还有一张“地图”，它在 5x5 网格上显示了 5 个区域的位置。例如

map = np.array([[1,1,1,2,2],
               [1,1,2,2,2],
               [3,3,3,2,2],
               [3,3,4,4,4],
               [3,3,5,5,4]])

我想要的是最终得到一个数组 (result)，它在 5x5 网格（类似于 map）上显示每个时间步长 (3) 和指标 (4) 上的值。所以

>>>result.shape 
(3,4,5,5)

我希望这是有道理的。

Answer 1

如果我理解正确，values.shape 中的 5 与map.shape 中的(5, 5) 无关，而仅与map 中的值范围有关。首先要做的是使map 中的值遵循传统的python 索引约定：从零开始计数，而不是从一开始。

以下是如何使用形状与值不同的地图：

import numpy as np

m, n, p = 2, 3, 4
values = np.arange(m*n*p).reshape(m, n, p)

map = np.array([
     [0, 0, 0, 3, 3, 3],
     [0, 0, 1, 3, 3, 2],
     [1, 1, 1, 3, 2, 2],
     [2, 1, 2, 2, 2, 1],
     [2, 2, 2, 2, 1, 1]
     ])

q, r = map.shape

map_r = map.reshape(q*r)

# result_r: shape (q*r, m, n)
# If you want to keep the 1-based numbering, index as [:, :, i-1].
result_r = np.array([values[:, :, i] for i in map_r])
result = result_r.transpose(1, 2, 0).reshape(m, n, q, r)

# test
jq, jr = 4, 5
mapval = map[jq, jr]
print(f'map[{jq}, {jr}] = {mapval}')
print(f'result[:, :, {jq}, {jr}] =\n{result[:, :, jq, jr]}')
print(f'values[:, :, {mapval}] =\n{values[:, :, mapval]}')

输出：

map[4, 5] = 1
result[:, :, 4, 5] =
[[ 1  5  9]
 [13 17 21]]
values[:, :, 1] =
[[ 1  5  9]
 [13 17 21]]

顺便说一句：如果您定义最左边的索引来表示更高级别的概念，多维数组通常会更清晰地工作。如果您可以将其表述为“M 个小部件的 N 个列表的列表”，那么自然数组形状将是 (N, M)，而不是 (M,N)。因此，您将拥有values.shape = (p, m, n) 和result.shape (q,r,m,n)（也可能交换m 和n，这取决于您稍后将如何处理数据）。你会避免很多[:, :] 索引和transpose 操作。

Answer 2

这是一个使用 np.pad 的实现：

values = np.arange(60).reshape(3,4,5)

# I renamed map to mapping, because map is a basic function in python

mapping = np.array([[1,1,1,2,2],
                  [1,1,2,2,2],
                  [3,3,3,2,2],
                  [3,3,4,4,4],
                  [3,3,5,5,4]])

v0,v1,v2 = values.shape
m0,m1 = mapping.shape

new = np.zeros((v0,v1,m0,m1))

for a in range(v2):
    q = values[:,:,a].reshape(v0,v1,1)
    r,c = np.where((mapping-1)==a)
    npad = ((0,0), (0,0), (0,len(c)-1))
    new[:,:,r,c] = np.pad(q, pad_width=npad, mode='symmetric')