使用NumPy.where查找数组坐标遇阻及跨数组匹配问题求助

针对你遇到的浮点数精度匹配失败、数组形状不匹配引发警告这两个核心问题，以下是具体解决方案：

一、先处理浮点数精度问题 #

浮点数的存储特性导致直接用==精确匹配极易出错，推荐两种可靠处理方式：

• 四舍五入后匹配：对数组保留指定小数位数后再比较（你已尝试的方案，适合精度明确的场景）
• 容差范围匹配：用np.allclose设置合理容差，判断坐标是否在误差范围内（更稳妥，避免四舍五入的截断误差）

二、解决数组形状不匹配的比较问题 #

由于coordsRAS是(30435615,3)的大数组，points是更小的(M,3)数组，直接==会触发广播错误，以下三种方案按需选择：

方案1：广播+容差比较（适合小规模points） #

通过扩展数组维度实现广播，逐点判断是否匹配：

import numpy as np

# 设置精度容差，根据数据实际情况调整，比如1e-8
tolerance = 1e-8

# 扩展维度实现广播，得到(N,M)的布尔数组，每行表示coordsRAS当前点是否匹配points中的某一个
matches = np.allclose(coordsRAS[:, np.newaxis, :], points[:, :3], atol=tolerance)
# 对每行取逻辑或，得到(N,)的布尔数组，标记coordsRAS中存在于points的点
result = np.any(matches, axis=1)

如果用四舍五入方式，代码类似：

coords_rounded = np.round(coordsRAS, 8)
points_rounded = np.round(points[:, :3], 8)

matches = np.all(coords_rounded[:, np.newaxis, :] == points_rounded, axis=2)
result = np.any(matches, axis=1)

方案2：结构化数组+np.isin（适合中等规模points） #

将多维坐标转为结构化数组，利用np.isin直接匹配：

# 先处理精度，再转换为结构化数组
coords_rounded = np.round(coordsRAS, 8)
points_rounded = np.round(points[:, :3], 8)

coords_struct = coords_rounded.view('f8,f8,f8').reshape(-1)
points_struct = points_rounded.view('f8,f8,f8').reshape(-1)

# 直接判断每个坐标是否存在于points中
result = np.isin(coords_struct, points_struct)

方案3：KDTree高效查找（适合超大规模数组） #

针对3000多万行的coordsRAS，广播比较效率极低，用KDTree空间索引可大幅提速：

from scipy.spatial import KDTree

tolerance = 1e-8
# 构建points的KDTree索引
tree = KDTree(points[:, :3])
# 查询每个coordsRAS点的最近邻，返回距离和对应索引
distances, indices = tree.query(coordsRAS, k=1, distance_upper_bound=tolerance)
# 距离小于等于容差即为匹配点
result = distances <= tolerance

注意事项 #

• 容差数值需根据原始数据的精度调整，比如数据保留8位小数时，用1e-8即可
• 当points规模较大时，KDTree方案的时间复杂度优势会非常明显，避免O(N*M)的低效计算

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