astro py:docs

- Index
- Modules

### Navigation

- next »
- « previous |
- Astropy v5.2.1 »
- Astronomical Coordinate Systems (`astropy.coordinates`) »
- Separations, Offsets, Catalog Matching, and Related Functionality

# 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 the `astropy.coordinates.BaseCoordinateFrame.separation()` or `astropy.coordinates.SkyCoord.separation()` methods, 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
    0.36208800460262563
    >>> sep.hour
    1.3830742984029318
    >>> sep.arcminute
    1244.7668685626384
    >>> sep.arcsecond
    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, `astropy.coordinates.BaseCoordinateFrame.separation_3d()` or `astropy.coordinates.SkyCoord.separation_3d()` methods 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:

    >>> matches = catalog[idx]
    >>> (matches.separation_3d(c) == d3d).all()
    True
    >>> 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)
    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(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]