Separations, Offsets, 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.
In the examples below, we will assume that the following imports have already been executed:
>>> import astropy.units as u
>>> from astropy.coordinates import SkyCoord
Separations#
The on-sky separation can be computed with
separation()
,
which computes the great-circle distance (not the small-angle approximation):
>>> c1 = SkyCoord('5h23m34.5s', '-69d45m22s', frame='icrs')
>>> c2 = SkyCoord('0h52m44.8s', '-72d49m43s', frame='fk5')
>>> sep = c1.separation(c2)
>>> sep
<Angle 20.74611448 deg>
The returned object is an Angle
instance, so it
is possible to access the angle in any of several equivalent angular
units:
>>> sep.radian
np.float64(0.36208800460262563)
>>> sep.hour
np.float64(1.3830742984029318)
>>> sep.arcminute
np.float64(1244.7668685626384)
>>> sep.arcsecond
np.float64(74686.0121137583)
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,
separation_3d()
will
determine the 3D distance between two coordinates that have distance
defined:
>>> 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.74398816 kpc>
Offsets#
Closely related to angular separations are offsets between coordinates. The key
distinction for offsets is generally the concept of a “from” and “to” coordinate
rather than the single scalar angular offset of a separation.
coordinates
contains conveniences to compute some of the common
offsets encountered in astronomy.
The first piece of such functionality is the
position_angle()
method.
This method computes the position angle between one
SkyCoord
instance and another (passed as the argument) following the
astronomy convention (positive angles East of North):
>>> c1 = SkyCoord(1*u.deg, 1*u.deg, frame='icrs')
>>> c2 = SkyCoord(2*u.deg, 2*u.deg, frame='icrs')
>>> c1.position_angle(c2).to(u.deg)
<Angle 44.97818294 deg>
The combination of separation()
and position_angle()
thus give a
set of directional offsets.
To do the inverse operation — determining the new
“destination” coordinate given a separation and position angle — the
directional_offset_by()
method is provided:
>>> c1 = SkyCoord(1*u.deg, 1*u.deg, frame='icrs')
>>> position_angle = 45 * u.deg
>>> separation = 1.414 * u.deg
>>> c1.directional_offset_by(position_angle, separation)
<SkyCoord (ICRS): (ra, dec) in deg
(2.0004075, 1.99964588)>
This technique is also useful for computing the midpoint (or indeed any point) between two coordinates in a way that accounts for spherical geometry (i.e., instead of averaging the RAs/Decs separately):
>>> coord1 = SkyCoord(0*u.deg, 0*u.deg, frame='icrs')
>>> coord2 = SkyCoord(1*u.deg, 1*u.deg, frame='icrs')
>>> pa = coord1.position_angle(coord2)
>>> sep = coord1.separation(coord2)
>>> coord1.directional_offset_by(pa, sep/2)
<SkyCoord (ICRS): (ra, dec) in deg
(0.49996192, 0.50001904)>
There is also a spherical_offsets_to()
method for computing angular offsets (e.g., small shifts like you might give a
telescope operator to move from a bright star to a fainter target):
>>> bright_star = SkyCoord('8h50m59.75s', '+11d39m22.15s', frame='icrs')
>>> faint_galaxy = SkyCoord('8h50m47.92s', '+11d39m32.74s', frame='icrs')
>>> dra, ddec = bright_star.spherical_offsets_to(faint_galaxy)
>>> dra.to(u.arcsec)
<Angle -173.78873354 arcsec>
>>> ddec.to(u.arcsec)
<Angle 10.60510342 arcsec>
The conceptual inverse of
spherical_offsets_to()
is also available as
a method on any SkyCoord
object:
spherical_offsets_by()
, which accepts two
angular offsets (in longitude and latitude) and returns the coordinates at the
offset location:
>>> target_star = SkyCoord(86.75309*u.deg, -31.5633*u.deg, frame='icrs')
>>> target_star.spherical_offsets_by(1.3*u.arcmin, -0.7*u.arcmin)
<SkyCoord (ICRS): (ra, dec) in deg
(86.77852168, -31.57496415)>
“Sky Offset” Frames#
To extend the concept of spherical offsets, coordinates
has
a frame class SkyOffsetFrame
which creates distinct frames that are centered on a specific point.
These are known as “sky offset frames,” as they are a convenient way to create
a frame centered on an arbitrary position on the sky suitable for computing
positional offsets (e.g., for astrometry):
>>> from astropy.coordinates import SkyOffsetFrame, ICRS
>>> center = ICRS(10*u.deg, 45*u.deg)
>>> center.transform_to(SkyOffsetFrame(origin=center))
<SkyOffsetICRS Coordinate (rotation=0.0 deg, origin=<ICRS Coordinate: (ra, dec) in deg
(10., 45.)>): (lon, lat) in deg
(0., 0.)>
>>> target = ICRS(11*u.deg, 46*u.deg)
>>> target.transform_to(SkyOffsetFrame(origin=center))
<SkyOffsetICRS Coordinate (rotation=0.0 deg, origin=<ICRS Coordinate: (ra, dec) in deg
(10., 45.)>): (lon, lat) in deg
(0.69474685, 1.00428706)>
Alternatively, the convenience method
skyoffset_frame()
lets you create a sky
offset frame from an existing SkyCoord
:
>>> center = SkyCoord(10*u.deg, 45*u.deg)
>>> aframe = center.skyoffset_frame()
>>> target.transform_to(aframe)
<SkyOffsetICRS Coordinate (rotation=0.0 deg, origin=<ICRS Coordinate: (ra, dec) in deg
(10., 45.)>): (lon, lat) in deg
(0.69474685, 1.00428706)>
>>> other = SkyCoord(9*u.deg, 44*u.deg, frame='fk5')
>>> other.transform_to(aframe)
<SkyCoord (SkyOffsetICRS: rotation=0.0 deg, origin=<ICRS Coordinate: (ra, dec) in deg
(10., 45.)>): (lon, lat) in deg
(-0.71943945, -0.99556216)>
Note
While sky offset frames appear to be all the same class, this not the
case: the sky offset frame for each different type of frame for origin
is
actually a distinct class. E.g., SkyOffsetFrame(origin=ICRS(...))
yields an object of class SkyOffsetICRS
, not SkyOffsetFrame
.
While this is not important for most uses of this class, it is important for
things like type-checking, because something like
SkyOffsetFrame(origin=ICRS(...)).__class__ is SkyOffsetFrame
will
not be True
, as it would be for most classes.
This same frame is also useful as a tool for defining frames that are relative to a specific, known object useful for hierarchical physical systems like galaxy groups. For example, objects around M31 are sometimes shown in a coordinate frame aligned with standard ICRA RA/Dec, but on M31:
>>> m31 = SkyCoord(10.6847083*u.deg, 41.26875*u.deg, frame='icrs')
>>> ngc147 = SkyCoord(8.3005*u.deg, 48.5087389*u.deg, frame='icrs')
>>> ngc147_inm31 = ngc147.transform_to(m31.skyoffset_frame())
>>> xi, eta = ngc147_inm31.lon, ngc147_inm31.lat
>>> xi
<Longitude -1.59206948 deg>
>>> eta
<Latitude 7.26183757 deg>
Note
Currently, distance information in the origin
of a
SkyOffsetFrame
is not
used to compute any part of the transform. The origin
is only used for
on-sky rotation. This may change in the future, however.
Matching Catalogs#
coordinates
leverages the coordinate framework to make it
possible 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:
>>> 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)
The distances returned d3d
are 3-dimensional distances.
Unless both source (c
) and catalog (catalog
) coordinates have
associated distances, this quantity assumes that all sources are at a distance
of 1 (dimensionless).
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:
>>> d3d
<Quantity [1335.55538257] kpc>
>>> matches = catalog[idx]
>>> matches.separation_3d(c)
<Distance [1335.55538257] kpc>
>>> dra, ddec = c.spherical_offsets_to(matches)
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)
It is possible to impose a separation constraint (e.g., the maximum separation to be
considered a match) by creating a boolean mask with d2d
or d3d
. For example:
>>> max_sep = 1.0 * u.arcsec
>>> idx, d2d, d3d = c.match_to_catalog_3d(catalog)
>>> sep_constraint = d2d < max_sep
>>> c_matches = c[sep_constraint]
>>> catalog_matches = catalog[idx[sep_constraint]]
Now, c_matches
and catalog_matches
are the matched sources in c
and catalog
, respectively, which are separated by less than 1 arcsecond.
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_*
:
>>> import numpy as np
>>> idxc, idxcatalog, d2d, d3d = catalog.search_around_sky(c, 1*u.deg)
>>> np.all(d2d < 1*u.deg)
np.True_
>>> idxc, idxcatalog, d2d, d3d = catalog.search_around_3d(c, 1*u.kpc)
>>> np.all(d3d < 1*u.kpc)
np.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)
np.True_
>>> np.all(c[idxc].separation_3d(catalog[idxcatalog]) == d3d)
np.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(ra=1*u.deg, dec=2*u.deg, distance=distance1*u.kpc)
>>> idxscalarc, idxcatalog, d2d, d3d = catalog.search_around_sky(scalarc, 1*u.deg)
ValueError: One of the inputs to search_around_sky is a scalar.
As a result (and because the search_around_*
algorithm is inefficient in
the scalar case), 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]