如何匹配三个不完全相等的数组中的坐标?

Question

我有六个数组，其中包含天文图像中检测到的源的坐标。每个数组在观测值的特定过滤器 (r，i，z) 中具有源的正确赤经或赤纬。仓位不会 100% 相等(浮点数)。 我发现可以找到整数的匹配对，可以检查相等性(如何在两个坐标数组中找到相同的 x，y 坐标？ 但是如何找到浮点数的匹配对呢？我的想法是检查两个位置之间的距离是否小于距离 dist=0.5/3600(所以 0.5 角秒)。然后我认为这两个来源是同一个对象。一旦找到三个过滤器中的匹配对，我想将有关该源(幅度)的一些信息附加到列表中。 这是我尝试过的代码：

# open the relevant catalogs
cat_r = pyfits.open("%s/%s_R_for_calib.cat" %(cluster,cluster))[2].data
cat_i = pyfits.open("%s/%s_I_for_calib.cat" %(cluster,cluster))[2].data
cat_z = pyfits.open("%s/%s_Z_for_calib.cat" %(cluster,cluster))[2].data
mask_r = cat_r['FLAGS'] == 0
mask_i = cat_i['FLAGS'] == 0
mask_z = cat_z['FLAGS'] == 0
# get the information needed for the analysis
cat_r_MAG               = cat_r.field("MAG_AUTO")[mask_r]
cat_i_MAG               = cat_i.field("MAG_AUTO")[mask_i]
cat_z_MAG               = cat_z.field("MAG_AUTO")[mask_z]
cat_r_MAG_ERR           = cat_r.field("MAGERR_AUTO")[mask_r]
cat_i_MAG_ERR           = cat_i.field("MAGERR_AUTO")[mask_i]
cat_z_MAG_ERR           = cat_z.field("MAGERR_AUTO")[mask_z]
cat_r_ra = cat_r.field("ALPHA_J2000")[mask_r]
cat_r_dec = cat_r.field("DELTA_J2000")[mask_r]
cat_i_ra = cat_i.field("ALPHA_J2000")[mask_i]
cat_i_dec = cat_i.field("DELTA_J2000")[mask_i]
cat_z_ra = cat_z.field("ALPHA_J2000")[mask_z]
cat_z_dec = cat_z.field("DELTA_J2000")[mask_z]
dist = 0.5 / 3600
matches = []
for a in range(0,len(cat_r_ra)):
for b in range(0,len(cat_i_ra)):
for c in range(0,len(cat_z_ra)):
if np.sqrt((cat_r_ra[a] - cat_i_ra[b])**2 + (cat_r_dec[a] - cat_i_dec[b])**2) < dist and np.sqrt((cat_r_ra[a] - cat_z_ra[c])**2 + (cat_r_dec[a] - cat_z_dec[c])**2 < dist):
match = [cat_r_ra[a], cat_r_dec[a], cat_r_MAG[a], cat_r_MAG_ERR[a], cat_i_MAG[b], cat_i_MAG_ERR[b], cat_z_MAG[c], cat_z_MAG_ERR[c]] 
matches.append(match)
matches = np.array(matches)

问题是，运行需要很长时间。我认为 for 循环和 if 语句的组合会减慢它的速度。但是我对python编程仍然相当陌生，所以我不知道更有效的方法。python 中是否有一种有效的方法来找到匹配的源，达到一定的精度(在本例中为 0.5/3600)？

感谢您的帮助！

Answer 1

您可以尝试事先对数组进行舍入。使用您链接的线程中建议的方法，这将为您提供：

import numpy as np
A = np.array([[ 2.31,  1, 4.33],
[ 4.46,  3, 10.75]])
B = np.array([[ 2.34,  3,  4.34],
[ 4.47, 11, 10.78]])
roundedA = np.around(A, decimals=1)
roundedB = np.around(B, decimals=1)
print('%sn%sn' % (roundedA, roundedB))
print(set(zip(*roundedA)) & set(zip(*roundedB)))

输出：

[[  2.3   1.    4.3]
[  4.5   3.   10.8]]
[[  2.3   3.    4.3]
[  4.5  11.   10.8]]
{(4.2999999999999998, 10.800000000000001), (2.2999999999999998, 4.5)}

当心浮点错误。