Separations, Catalog Matching, and Related Functionality¶
astropy.coordinates
contains commonly-used tools for comparing or
matching coordinate objects. Of particular importance are those for
determining separations between coordinates and those for matching a
coordinate (or coordinates) to a catalog. These are mainly implemented
as methods on the coordinate objects.
Separations¶
The on-sky separation is easily computed with the
astropy.coordinates.BaseCoordinateFrame.separation()
or
astropy.coordinates.SkyCoord.separation()
methods,
which computes the great-circle distance (not the small-angle
approximation):
>>> import numpy as np
>>> from astropy import units as u
>>> from astropy.coordinates import SkyCoord
>>> c1 = SkyCoord('5h23m34.5s', '-69d45m22s', frame='icrs')
>>> c2 = SkyCoord('0h52m44.8s', '-72d49m43s', frame='fk5')
>>> sep = c1.separation(c2)
>>> sep
<Angle 20.74611447604398 deg>
The returned object is an Angle
instance, so it
is straightforward to access the angle in any of several equivalent angular
units:
>>> sep.radian
0.36208800460262575
>>> sep.hour
1.3830742984029323
>>> sep.arcminute
1244.7668685626388
>>> sep.arcsecond
74686.01211375833
Also note that the two input coordinates were not in the same frame - one is automatically converted to match the other, ensuring that even though they are in different frames, the separation is determined consistently.
In addition to the on-sky separation described above,
astropy.coordinates.BaseCoordinateFrame.separation_3d()
or
astropy.coordinates.SkyCoord.separation_3d()
methods will
determine the 3D distance between two coordinates that have distance
defined:
>>> from astropy.coordinates import SkyCoord
>>> c1 = SkyCoord('5h23m34.5s', '-69d45m22s', distance=70*u.kpc, frame='icrs')
>>> c2 = SkyCoord('0h52m44.8s', '-72d49m43s', distance=80*u.kpc, frame='icrs')
>>> sep = c1.separation_3d(c2)
>>> sep
<Distance 28.743988157814094 kpc>
Matching Catalogs¶
coordinates
supports leverages the coordinate framework to make it
straightforward to find the closest coordinates in a catalog to a desired set
of other coordinates. For example, assuming ra1
/dec1
and
ra2
/dec2
are numpy arrays loaded from some file:
>>> from astropy.coordinates import SkyCoord
>>> from astropy import units as u
>>> c = SkyCoord(ra=ra1*u.degree, dec=dec1*u.degree)
>>> catalog = SkyCoord(ra=ra2*u.degree, dec=dec2*u.degree)
>>> idx, d2d, d3d = c.match_to_catalog_sky(catalog)
You can also find the nearest 3d matches, different from the on-sky
separation shown above only when the coordinates were initialized with
a distance
:
>>> c = SkyCoord(ra=ra1*u.degree, dec=dec1*u.degree, distance=distance1*u.kpc)
>>> catalog = SkyCoord(ra=ra2*u.degree, dec=dec2*u.degree, distance=distance2*u.kpc)
>>> idx, d2d, d3d = c.match_to_catalog_3d(catalog)
Now idx
are indices into catalog
that are the closest objects to each
of the coordinates in c
, d2d
are the on-sky distances between them, and
d3d
are the 3-dimensional distances. Because coordinate objects support
indexing, idx
enables easy access to the matched set of coordinates in
the catalog:
>>> matches = catalog[idx]
>>> (matches.separation_3d(c) == d3d).all()
True
>>> dra = (matches.ra - c.ra).arcmin
>>> ddec = (matches.dec - c.dec).arcmin
This functionality can also be accessed from the
match_coordinates_sky()
and
match_coordinates_3d()
functions. These
will work on either SkyCoord
objects or the lower-level frame classes:
>>> from astropy.coordinates import match_coordinates_sky
>>> idx, d2d, d3d = match_coordinates_sky(c, catalog)
>>> idx, d2d, d3d = match_coordinates_sky(c.frame, catalog.frame)
Searching Around Coordinates¶
Closely-related functionality can be used to search for all coordinates within
a certain distance (either 3D distance or on-sky) of another set of coordinates.
The search_around_*
methods (and functions) provide this functionality,
with an interface very similar to match_coordinates_*
:
>>> idxc, idxcatalog, d2d, d3d = catalog.search_around_sky(c, 1*u.deg)
>>> np.all(d2d < 1*u.deg)
True
>>> idxc, idxcatalog, d2d, d3d = catalog.search_around_3d(c, 1*u.kpc)
>>> np.all(d3d < 1*u.kpc)
True
The key difference for these methods is that there can be multiple (or no)
matches in catalog
around any locations in c
. Hence, indices into both
c
and catalog
are returned instead of just indices into catalog
.
These can then be indexed back into the two SkyCoord
objects, or, for that
matter, any array with the same order:
>>> np.all(c[idxc].separation(catalog[idxcatalog]) == d2d)
True
>>> np.all(c[idxc].separation_3d(catalog[idxcatalog]) == d3d)
True
>>> print catalog_objectnames[idxcatalog]
['NGC 1234' 'NGC 4567' ...]
Note, though, that this dual-indexing means that search_around_*
does not
work well if one of the coordinates is a scalar, because the returned index
would not make sense for a scalar:
>>> scalarc = SkyCoord(1*u.deg, 2*u.deg)
>>> idxscalarc, idxcatalog, d2d, d3d = catalog.search_around_sky(scalarc, 1*u.deg) # THIS DOESN'T ACTUALLY WORK
>>> scalarc[idxscalarc]
IndexError: 0-d arrays can't be indexed
As a result (and because the search_around_*
algorithm is inefficient in
the scalar case, anyway), the best approach for this scenario is to instead
use the separation*
methods:
>>> d2d = scalarc.separation(catalog)
>>> catalogmsk = d2d < 1*u.deg
>>> d3d = scalarc.separation_3d(catalog)
>>> catalog3dmsk = d3d < 1*u.kpc
The resulting catalogmsk
or catalog3dmsk
variables are boolean arrays
rather than arrays of indices, but in practice they usually can be used in
the same way as idxcatalog
from the above examples. If you definitely do
need indices instead of boolean masks, you can do:
>>> idxcatalog = np.where(catalogmsk)[0]
>>> idxcatalog3d = np.where(catalog3dmsk)[0]