I am using Astropy.coordinates to match two astronomical catalogues with RA,DEC coordinates.
I can find the list of nearest neighbours by following the astropy documentation(link to astropy.coordinates documentation) and doing:
from astropy.coordinates import SkyCoord
from astropy import units as u
cat1 = SkyCoord(ra=ra1*u.degree, dec=dec1*u.degree)
cat2 = SkyCoord(ra=ra2*u.degree, dec=dec2*u.degree)
idx, d2d, d3d = cat1.match_to_catalog_sky(cat2)
Where ra1, ra2, dec1, dec2 are vectors containing the coordinates in the catalogues 1 and 2.
The result idx gives, for each object in the catalogue 1, the id of nearest match in the catalogue 2. d2d gives the 2d separation between the matches, and d3d gives the 3d separation between the matches.
Therefore, to select matches between a desired matching radius, for example, using a 1" radii, I can do:
matched=idx[np.argwhere(d2d<1.*u.arcsec)[0]]
Now, in order to chose what is the appropriate radii for this last step, I would like to examine what is the distance d2d between each source in cat1and their second-nearest-neighbour.
Does anyone knows how can I do this matching process while also recording the second neighbours?
Note that match_to_catalog_sky takes a keyword nthneighbor which defaults to 1 (i.e., nearest) but can be set to other values. Changing this should do it.
It seems like you may also want to pay attention to search_around_sky method using which you can find all matches up to a separation limit.