0

I have six arrays with coordinates of detected sources in an astronomical image. Each array has the right ascension or declination of the sources in a specific filter of the observation (r,i,z). The positions will not be 100% equal (float numbers). I found out that one can find matching pairs for integer numbers where it is possible to check for equality (How can I find identical x,y coordinates in two arrays of corrdinates?). But how do I find the matching pairs for float numbers? My idea was to check if the distance between two positions is smaller than a distance dist=0.5 / 3600 (so 0.5 arcsec). Then I consider the two sources to be the same object. Once the matching pairs in three filters are found I want to append some information on that source (the magnitudes) to a list. Here is the code that I tried:

# 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)

The problem is, that it takes really long to run. I think the combination of for loops and if statements slows it down. But I am still fairly new to programming in python so I do not know a more efficient way. Is there an efficient way in python to find the matching sources up to some accuracy (in this case 0.5 / 3600)?

Thank you for your help!

Hannah Z.
  • 1
  • You can start by converting both arrays to numpy and then subtracting them; now you just have to find the small entries of the difference array. – perigon Jul 05 '17 at 08:06

1 Answers1

0

You could try rounding your arrays beforehand. Using the method suggested in the thread you linked, this would give you:

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('%s\n%s\n' % (roundedA, roundedB))
print(set(zip(*roundedA)) & set(zip(*roundedB)))

Output:

[[  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)}

Beware of the floating point errors.

Flomp
  • 524
  • 4
  • 17