如何从3D NumPy数组中按2D数组的Z值提取对应点？

你的问题核心是索引数组的形状无法广播——你用的levels是展平后的一维数组（721200个元素），但y和x分别是(601,)和(1200,)的一维数组，这三个数组的形状无法对齐，导致索引失败。

不需要循环，用NumPy的矢量化索引就能高效解决，关键是让三个索引数组的形状完全一致。这里有两种简洁的实现方式：

方法一：使用np.indices生成坐标索引 #

np.indices可以直接生成和你的2D数组形状一致的y、x坐标索引矩阵，这样三个索引（z值、y坐标、x坐标）的形状完全匹配，直接就能完成索引：

import numpy as np

# 假设你的数组已经定义：array2d (601,1200)，oldArray3d (137,601,1200)
y_indices, x_indices = np.indices(array2d.shape)

# 直接索引，得到的newArray2d形状就是(601,1200)
newArray2d = oldArray3d[array2d, y_indices, x_indices]

原理说明 #

np.indices(array2d.shape)会生成两个数组：

• y_indices：形状(601,1200)，每个位置的值是该点的y坐标（0到600）
• x_indices：形状(601,1200)，每个位置的值是该点的x坐标（0到1199）

当你用oldArray3d[array2d, y_indices, x_indices]时，NumPy会对每个位置(i,j)，取oldArray3d[array2d[i,j], i, j]，正好就是你需要的每个(y,x)点对应z层的值。

方法二：使用np.meshgrid生成坐标索引 #

如果习惯用meshgrid，记得加上indexing='ij'参数来匹配数组的(y,x)维度顺序：

y_coords = np.arange(601)
x_coords = np.arange(1200)

# 生成和array2d形状一致的坐标矩阵
y_indices, x_indices = np.meshgrid(y_coords, x_coords, indexing='ij')

# 索引并得到结果
newArray2d = oldArray3d[array2d, y_indices, x_indices]

这两种方法都是完全矢量化的，效率远高于循环，而且代码简洁易读。

内容的提问来源于stack exchange，提问作者Sev