使用高级索引：

result = reference[map[:,:,0], map[:,:,1]]

(BTW避免使用关键字map作为变量名...）
时间：

import timeit
%timeit test = np.array([[reference[i[0],i[1],:] for i in j] for j in mp])
# 4.53 µs ± 72.5 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
%timeit test = reference[mp[:,:,0], mp[:,:,1]]
# 1.64 µs ± 4.13 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)

在更大的阵列上，增益可能会更大...例如，对于map = np.random.randint(0,3,size=(100,100,2))，我们得到

4.11 ms ± 15.6 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

V.S.

98.5 µs ± 200 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

类似于@Julien的答案，但使用元组索引。你只需要把map_的尺寸按正确的顺序排列：

reference[tuple(np.moveaxis(map_, 2, 0))]

输出量：

array([[[11, 12, 13],
        [27, 28, 29],
        [37, 38, 39]],
       [[34, 35, 36],
        [47, 48, 49],
        [41, 42, 43]]])

也就是说，如果map_是：

array([[[0, 1, 2],
        [2, 3, 3]],
       [[0, 2, 2],
        [1, 2, 0]]])

编辑：朱利安的回答更简洁，应该被接受这个问题。我花了一些时间写这个答案，所以我把它作为1留在这里。我花了很多时间来写，2。其他人可能会在这里面找到一些东西（谁知道呢）
简短的形式在这里。找到下面的长形式：

reference = array([[ [11,12,13], [14,15,16], [17,18,19] ],
                   [ [21,22,23], [24,25,26], [27,28,29] ],
                   [ [31,32,33], [34,35,36], [37,38,39] ],
                   [ [41,42,43], [44,45,46], [47,48,49] ]])
idx = array([
    [[0, 0], [1, 2], [2, 2]],
    [[2, 1], [3, 2], [3, 0]]
])
xd, yd = idx.shape[0], idx.shape[1]
zd = reference.shape[2]
idx = idx.reshape((xd * yd, 2))
indexed_reference = reference[idx[:,0], idx[:, 1]].reshape((xd, yd, zd))
print(indexed_reference)

首先，将你的Map重塑为一个2d数组，以便于索引。注意，我将使用idx而不是map，因为map是python中预定义的函数：

idx = array([
    [[0, 0], [1, 2], [2, 2]],
    [[2, 1], [3, 2], [3, 0]]
])
# idx.shape -> (2, 3, 2), make this (6,2). not reshaped dimensions' product must be constant. 2x3x3 == 6x2
# you can probably also (multiply the first two, last) if you want to make it generic
reshaped_idx = idx.reshape((6,2))
# reshaped_idx = array([[0, 0], [1, 2], [2, 2], [2, 1], [3, 2], [3, 0]])
# now just use the indexes to map the original array
reference = array([[ [11,12,13], [14,15,16], [17,18,19] ],
                   [ [21,22,23], [24,25,26], [27,28,29] ],
                   [ [31,32,33], [34,35,36], [37,38,39] ],
                   [ [41,42,43], [44,45,46], [47,48,49] ]])
indexed_reference = reference[reshaped_idx[:,0], reshaped_idx[:,1]]
# array([
#      [11, 12, 13], [27, 28, 29], [37, 38, 39], 
#      [34, 35, 36], [47, 48, 49], [41, 42, 43]
# ])
# indexed_reference.shape -> (6,3) -> 6x3 = 18
# reshape this to the original idx
indexed_reference_reshaped = indexed_reference.reshape((2,3,3))
# answer: array([
#    [[11, 12, 13],[27, 28, 29], [37, 38, 39]],
#    [[34, 35, 36],[47, 48, 49], [41, 42, 43]]
# ])