astro py:docs

- Index
- Modules

### Navigation

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

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