Using astropy.coordinates to Match Catalogs and Plan Observations

In this tutorial, we will explore how the astropy.coordinates package and related astropy functionality can be used to help in planning observations or other exercises focused on large coordinate catalogs.

Note that this tutorial requires astropy v1.0 or higher.

In [1]:
# this cell imports functions from the python standard library in a way that works with both python 2 and 3
    # Python 3.x
    from urllib.parse import urlencode
    from urllib.request import urlretrieve
except ImportError:
    # Python 2.x
    from urllib import urlencode
    from urllib import urlretrieve
In [2]:
import numpy as np
import IPython.display
In [3]:
%matplotlib inline
from matplotlib import pyplot as plt
In [4]:
from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.table import Table

Describing on-sky locations with coordinates

Let's start by considering a field around the picturesque Hickson Compact Group 7. To do anything with this, we need to get an object that represents the coordinates of the center of this group.

In Astropy, the most common object you'll work with for coordinates is SkyCoord. A SkyCoord can be created most easily directly from angles as shown below. It's also wise to explicitly specify the frame your coordinates are in, although this is not strictly necessary because the default is ICRS.

(If you're not sure what ICRS is, it's basically safe to think of it as an approximation to an equatorial system at the J2000 equinox).

In [5]:
hcg7_center = SkyCoord(9.81625*u.deg, 0.88806*u.deg, frame='icrs')
<SkyCoord (ICRS): (ra, dec) in deg
    ( 9.81625,  0.88806)>

SkyCoord will also accept string-formatted coordinates either as separate strings for ra/dec or a single string. You'll have to give units, though, if they aren't part of the string itself.

In [6]:
SkyCoord('0h39m15.9s', '0d53m17.016s', frame='icrs')
<SkyCoord (ICRS): (ra, dec) in deg
    ( 9.81625,  0.88806)>
In [7]:
SkyCoord('0:39:15.9 0:53:17.016', unit=(u.hour, u.deg), frame='icrs')
<SkyCoord (ICRS): (ra, dec) in deg
    ( 9.81625,  0.88806)>

If the object you're interested in is in SESAME, you can also look it up directly from its name using the SkyCoord.from_name() class method1. Note that this requires an internet connection. It's safe to skip if you don't have one, because we defined it above explicitly.

1If you don't know what a class method is, think of it like an alternative constructor for a `SkyCoord` object -- calling `SkyCoord.from_name()` with a name gives you a new `SkyCoord` object. For more detailed background on what class methods are and when they're useful, see [this page](

In [8]:
hcg7_center = SkyCoord.from_name('HCG 7')
timeout                                   Traceback (most recent call last)
/Users/erik/miniconda3/envs/astropy-tutorials/lib/python3.5/site-packages/astropy/coordinates/ in get_icrs_coordinates(name)
    137             # Retrieve ascii name resolve data from CDS
--> 138             resp = urllib.request.urlopen(url, timeout=data.conf.remote_timeout)
    139             resp_data =

/Users/erik/miniconda3/envs/astropy-tutorials/lib/python3.5/urllib/ in urlopen(url, data, timeout, cafile, capath, cadefault, context)
    162         opener = _opener
--> 163     return, data, timeout)

/Users/erik/miniconda3/envs/astropy-tutorials/lib/python3.5/urllib/ in open(self, fullurl, data, timeout)
--> 466         response = self._open(req, data)

/Users/erik/miniconda3/envs/astropy-tutorials/lib/python3.5/urllib/ in _open(self, req, data)
    483         result = self._call_chain(self.handle_open, protocol, protocol +
--> 484                                   '_open', req)
    485         if result:

/Users/erik/miniconda3/envs/astropy-tutorials/lib/python3.5/urllib/ in _call_chain(self, chain, kind, meth_name, *args)
    443             func = getattr(handler, meth_name)
--> 444             result = func(*args)
    445             if result is not None:

/Users/erik/miniconda3/envs/astropy-tutorials/lib/python3.5/urllib/ in http_open(self, req)
   1281     def http_open(self, req):
-> 1282         return self.do_open(http.client.HTTPConnection, req)

/Users/erik/miniconda3/envs/astropy-tutorials/lib/python3.5/urllib/ in do_open(self, http_class, req, **http_conn_args)
   1256                 raise URLError(err)
-> 1257             r = h.getresponse()
   1258         except:

/Users/erik/miniconda3/envs/astropy-tutorials/lib/python3.5/http/ in getresponse(self)
   1196             try:
-> 1197                 response.begin()
   1198             except ConnectionError:

/Users/erik/miniconda3/envs/astropy-tutorials/lib/python3.5/http/ in begin(self)
    296         while True:
--> 297             version, status, reason = self._read_status()
    298             if status != CONTINUE:

/Users/erik/miniconda3/envs/astropy-tutorials/lib/python3.5/http/ in _read_status(self)
    257     def _read_status(self):
--> 258         line = str(self.fp.readline(_MAXLINE + 1), "iso-8859-1")
    259         if len(line) > _MAXLINE:

/Users/erik/miniconda3/envs/astropy-tutorials/lib/python3.5/ in readinto(self, b)
    574             try:
--> 575                 return self._sock.recv_into(b)
    576             except timeout:

timeout: timed out

During handling of the above exception, another exception occurred:

NameResolveError                          Traceback (most recent call last)
<ipython-input-8-a0c6eba392ac> in <module>()
----> 1 hcg7_center = SkyCoord.from_name('HCG 7')
      2 hcg7_center

/Users/erik/miniconda3/envs/astropy-tutorials/lib/python3.5/site-packages/astropy/coordinates/ in from_name(cls, name, frame)
   1284         from .name_resolve import get_icrs_coordinates
-> 1286         icrs_coord = get_icrs_coordinates(name)
   1287         icrs_sky_coord = cls(icrs_coord)
   1288         if frame in ('icrs', icrs_coord.__class__):

/Users/erik/miniconda3/envs/astropy-tutorials/lib/python3.5/site-packages/astropy/coordinates/ in get_icrs_coordinates(name)
    155             raise NameResolveError(
    156                 "Unable to retrieve coordinates for name '{0}'; connection "
--> 157                 "timed out".format(name))
    159     # All Sesame URL's timed out...

NameResolveError: Unable to retrieve coordinates for name 'HCG 7'; connection timed out

This object we just created has various useful ways of accessing the information contained within it. In particular, the ra and dec attributes are specialized Quantity objects (actually, a subclass called Angle, which in turn is subclassed by Latitude and Longitude). These objects store angles and provide pretty representations of those angles, as well as some useful attributes to quickly convert to common angle units:

In [9]:
type(hcg7_center.ra), type(hcg7_center.dec)
(astropy.coordinates.angles.Longitude, astropy.coordinates.angles.Latitude)
In [10]:
In [11]:
In [12]:

Now that we have a SkyCoord object, we can try to use it to access data from the Sloan Digitial Sky Survey (SDSS). Let's start by trying to get a picture using the SDSS image cutout service to make sure HCG7 is in the SDSS footprint and has good image quality.

This requires an internet connection, but if it fails, don't worry: the file is included in the repository so you can just let it use the local file'HCG7_SDSS_cutout.jpg', defined at the top of the cell.

In [13]:
impix = 1024
imsize = 12*u.arcmin
cutoutbaseurl = ''
query_string = urlencode(dict(ra=hcg7_center.ra.deg, 
                                     width=impix, height=impix, 
url = cutoutbaseurl + '?' + query_string

# this downloads the image to your disk
urlretrieve(url, 'HCG7_SDSS_cutout.jpg')
('HCG7_SDSS_cutout.jpg', <http.client.HTTPMessage at 0x114788358>)

Now lets take a look at the image.

In [14]:

Very pretty!


Create a SkyCoord of some other astronomical object you find interesting. Using only a single method/function call, get a string with the RA/Dec in the form 'HH:MM:SS.S DD:MM:SS.S'. Check your answer against an academic paper or some web site like SIMBAD that will show you sexigesimal coordinates for the object.

(Hint: SkyCoord.to_string() might be worth reading up on)

In [15]:

Now get an image of that object from the Digitized Sky Survey and download it and/or show it in the notebook. Bonus points if you figure out the (one-line) trick to get it to display in the notebook without ever downloading the file yourself.

(Hint: STScI has an easy-to-access copy of the DSS. The pattern to follow for the web URL is

In [15]:

Using coordinates and table to match and compare catalogs

At the end of the last section, we determined that HCG7 is in the SDSS imaging survey, so that means we can use the cells below to download catalogs of objects directly from the SDSS. Later on, we will match this catalog to another catalog covering the same field, allowing us to make plots using the combination of the two catalogs.

We will access the SDSS SQL database using the astroquery affiliated package. This will require an internet connection and a working install of astroquery. If you don't have these you can just skip down two cells, because the data files are provided with the repository. Depending on your version of astroquery it might also issue a warning, which you should be able to safely ignore.

In [15]:
from astroquery.sdss import SDSS
sdss = SDSS.query_region(coordinates=hcg7_center, radius=20*u.arcmin, 
ImportError                               Traceback (most recent call last)
<ipython-input-15-54321b6f2b92> in <module>()
----> 1 from astroquery.sdss import SDSS
      2 sdss = SDSS.query_region(coordinates=hcg7_center, radius=20*u.arcmin, 
      3                          spectro=True,
      4                          photoobj_fields=['ra','dec','u','g','r','i','z'])

ImportError: No module named 'astroquery'

astroquery queries gives us back an astropy.table.Table object. We could just work with this directly without saving anything to disk if we wanted to. But here we will use the capability to write to disk. That way, if you quit the session and come back later, you don't have to run the query a second time.

(Note that this won't work fail if you skipped the last step. Don't worry, you can just skip to the next cell with and use the copy of this table included in the tutorial.)

In [16]:
sdss.write('HCG7_SDSS_photo.dat', format='ascii')
NameError                                 Traceback (most recent call last)
<ipython-input-16-1d170dff9193> in <module>()
----> 1 sdss.write('HCG7_SDSS_photo.dat', format='ascii')

NameError: name 'sdss' is not defined

If you don't have internet, you can read the table into python by running the cell below. But if you did the astroquery step above, you could skip this, as the table is already in memory as the sdss variable.

In [17]:
sdss ='HCG7_SDSS_photo.dat', format='ascii')

Ok, so we have a catalog of objects we got from the SDSS. Now lets say you have your own catalog of objects in the same field that you want to match to this SDSS catalog. In this case, we will use a catalog extracted from the 2MASS. We first load up this catalog into python.

In [18]:
twomass ='HCG7_2MASS.tbl', format='ascii')

Now to do matching we need SkyCoord objects. We'll have to build these from the tables we loaded, but it turns out that's pretty straightforward: we grab the RA and dec columns from the table and provide them to the SkyCoord constructor. Lets first have a look at the tables to see just what everything is that's in them.

In [19]:
sdss # just to see an example of the format
<Table length=679>
In [20]:
twomass # just to see an example of the format
<Table masked=True length=23>

OK, looks like they both have ra and dec columns, so we should be able to use that to make SkyCoords.

You might first think you need to create a separate SkyCoord for every row in the table, given that up until now all SkyCoords we made were for just a single point. You could do this, but it will make your code much slower. Instead, SkyCoord supports arrays of coordinate values - you just pass in array-like inputs (array Quantitys, lists of strings, Table columns, etc.), and SkyCoord will happily do all of its operations element-wise.

In [21]:
coo_sdss = SkyCoord(sdss['ra']*u.deg, sdss['dec']*u.deg)
coo_twomass = SkyCoord(twomass['ra'], twomass['dec'])

Note a subtle difference here: you had to give units for SDSS but not for 2MASS. This is because the 2MASS table has units associated with the columns, while the SDSS table does not (so you have to put them in manually).

Now we simply use the SkyCoord.match_to_catalog_sky method to match the two catalogs. Note that order matters: we're matching 2MASS to SDSS because there are many more entires in the SDSS, so it seems likely that most 2MASS objects are in SDSS (but not vice versa).

In [22]:
idx_sdss, d2d_sdss, d3d_sdss = coo_twomass.match_to_catalog_sky(coo_sdss)

idx are the indecies into coo_sdss that get the closest matches, while d2d and d3d are the on-sky and real-space distances between the matches. In our case d3d can be ignored because we didn't give a line-of-sight distance, so its value is not particularly useful. But d2d provides a good diagnosis of whether we actually have real matches:

In [23]:
plt.hist(d2d_sdss.arcsec, histtype='step', range=(0,2))
plt.xlabel('separation [arcsec]')

Ok, they're all within an arcsecond that's promising. But are we sure it's not just that anything has matches within an arcescond? Lets check by comparing to a set of random points.

We first create a set of uniformly random points (with size matching coo_twomass) that cover the same range of RA/Decs that are in coo_sdss.

In [24]:
ras_sim = np.random.rand(len(coo_twomass))*coo_sdss.ra.ptp() + coo_sdss.ra.min()
decs_sim = np.random.rand(len(coo_twomass))*coo_sdss.dec.ptp() + coo_sdss.dec.min()
ras_sim, decs_sim
(<Angle [  9.81682259,  9.87369311, 10.09725714, 10.0240289 , 10.04789017,
           9.95885943, 10.07893784,  9.55378732,  9.60853007, 10.12051673,
          10.08499863, 10.04309114,  9.64250554,  9.88556643, 10.13256088,
          10.04541979,  9.78763987, 10.12931866,  9.78014949,  9.548045  ,
           9.71071039,  9.64622839,  9.61040161] deg>,
 <Angle [ 0.82888365, 1.00598845, 1.04969638, 0.68041027, 1.13098551,
          0.83062912, 0.8946654 , 0.80145967, 0.56940723, 0.9385579 ,
          0.60147316, 0.83023458, 0.77181668, 0.85354077, 0.56384891,
          0.58545737, 0.8223226 , 1.07372833, 0.75862946, 0.74663958,
          0.96571573, 1.19241496, 0.61242448] deg>)

Now we create a SkyCoord from these points and match it to coo_sdss just like we did above for 2MASS.

Note that we do not need to explicitly specify units for ras_sim and decs_sim, because they already are unitful Angle objects because they were created from coo_sdss.ra/coo_sdss.dec.

In [25]:
coo_simulated = SkyCoord(ras_sim, decs_sim)  
idx_sim, d2d_sim, d3d_sim = coo_simulated.match_to_catalog_sky(coo_sdss)

Now lets plot up the histogram of separations from our simulated catalog so we can compare to the above results from the real catalog.

In [26]:
plt.hist(d2d_sim.arcsec, histtype='step', color='red', label='Simulated', linestyle='dashed')
plt.hist(d2d_sdss.arcsec, histtype='step', color='blue', label='2MASS')
plt.xlabel('separation [arcsec]')

Alright, great - looks like randomly placed sources should be more like an arcminute away, so we can probably trust that our earlier matches which were within an arcsecond are valid. So with that in mind, we can start computing things like colors that combine the SDSS and 2MASS photometry.

In [27]:
rmag = sdss['r'][idx_sdss]
grcolor = sdss['g'][idx_sdss] - rmag
rKcolor = rmag - twomass['k_m_ext']

plt.subplot(1, 2, 1)
plt.scatter(rKcolor, rmag)
plt.xlim(2.5, 4)
plt.ylim(18, 12) #mags go backwards!

plt.subplot(1, 2, 2)
plt.scatter(rKcolor, rmag)
plt.xlim(2.5, 4)


For more on what matching options are available, check out the separation and matching section of the astropy documentation. Or for more on what you can do with SkyCoord, see its API documentation.


Check that the d2d_sdss variable matches the on-sky separations you get from comaparing the matched coo_sdss entries to coo_twomass.

Hint: You'll likely find the SkyCoord.separation() method useful here.

In [28]:

Compute the physical separation between two (or more) objects in the catalogs. You'll need line-of-sight distances, so a reasonable guess might be the distance to HCG 7, which is about 55 Mpc.

Hint: you'll want to create new SkyCoord objects, but with distance attributes. There's also a SkyCoord method that should do the rest of the work, but you'll have to poke around to figure out what it is.

In [28]:

Transforming between coordinate systems and planning observations

Now lets say something excites you about one of the objects in this catalog, and you want to know if and when you might go about observing it. astropy.coordinates provides tools to enable this, as well.

Introducting frame transformations

To understand the code in this section, it may help to read over the overview of the astropy coordinates scheme. The key bit to understand is that all coordinates in astropy are in particular "frames", and we can transform between a specific SkyCoord object from one frame to another. For example, we can transform our previously-defined center of HCG7 from ICRS to Galactic coordinates:

In [28]:
<SkyCoord (Galactic): (l, b) in deg
    ( 116.47556813, -61.83099472)>

The above is actually a special "quick-access" form which internally does the same as what's in the cell below: uses the transform_to() method to convert from one frame to another.

In [29]:
from astropy.coordinates import Galactic
<SkyCoord (Galactic): (l, b) in deg
    ( 116.47556813, -61.83099472)>

Note that changing frames also changes some of the attributes of the object, but usually in a way that makes sense:

In [30]:
hcg7_center.galactic.ra  # should fail because galactic coordinates are l/b not RA/Dec
AttributeError                            Traceback (most recent call last)
<ipython-input-30-d7bc134707f6> in <module>()
----> 1 hcg7_center.galactic.ra  # should fail because galactic coordinates are l/b not RA/Dec

/Users/erik/miniconda3/envs/astropy-tutorials/lib/python3.5/site-packages/astropy/coordinates/ in __getattr__(self, attr)
    485         # Fail
    486         raise AttributeError("'{0}' object has no attribute '{1}'"
--> 487                              .format(self.__class__.__name__, attr))
    489     def __setattr__(self, attr, val):

AttributeError: 'SkyCoord' object has no attribute 'ra'
In [31]:

Using frame transformations to get to AltAz

To actually do anything with observability we need to convert to a frame local to an on-earth observer. By far the most common choice is horizontal coordinates, or "AltAz" coordinates. We first need to specify both where and when we want to try to observe.

In [32]:
from astropy.coordinates import EarthLocation
from astropy.time import Time

observing_location = EarthLocation(lat='31d57.5m', lon='-111d35.8m', height=2096*u.m)  # Kitt Peak, Arizona
# If you're using astropy v1.1 or later, you can replace the above with this:
#observing_location = EarthLocation.of_site('Kitt Peak')

observing_time = Time('2010-12-21 1:00')  # 1am UTC=6pm AZ mountain time

Now we use these to create an AltAz frame object. Note that this frame has some other information about the atmosphere, which can be used to correct for atmospheric refraction. Here we leave that alone, because the default is to ignore this effect (by setting the pressure to 0).

In [33]:
from astropy.coordinates import AltAz

aa = AltAz(location=observing_location, obstime=observing_time)
<AltAz Frame (obstime=2010-12-21 01:00:00.000, location=(-1994310.09211632, -5037908.606337594, 3357621.752122168) m, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0, obswl=1.0 micron)>

Now we can just transform our ICRS SkyCoord to AltAz to get the location in the sky over Kitt Peak at the requested time.

In [34]:
<SkyCoord (AltAz: obstime=2010-12-21 01:00:00.000, location=(-1994310.09211632, -5037908.606337594, 3357621.752122168) m, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0, obswl=1.0 micron): (az, alt) in deg
    ( 149.19392032,  55.06247359)>

Alright, it's up at 6pm, but that's pretty early to be observing. We could just try various times one at a time to see if the airmass is at a darker time, but we can do better: lets try to create an airmass plot.

In [35]:
# this gives a Time object with an *array* of times
delta_hours = np.linspace(0, 6, 100)*u.hour
full_night_times = observing_time + delta_hours
full_night_aa_frames = AltAz(location=observing_location, obstime=full_night_times)
full_night_aa_coos = hcg7_center.transform_to(full_night_aa_frames)

plt.plot(delta_hours, full_night_aa_coos.secz)
plt.xlabel('Hours from 6pm AZ time')
plt.ylabel('Airmass [Sec(z)]')

Great! Looks like it's at the lowest airmass in another hour or so (7pm). But might that might still be twilight... When should we start observing for proper dark skies? Fortunately, astropy provides a get_sun function that can be used to check this. Lets use it to check if we're in 18-degree twilight or not.

In [36]:
from astropy.coordinates import get_sun

full_night_sun_coos = get_sun(full_night_times).transform_to(full_night_aa_frames)
plt.plot(delta_hours, full_night_sun_coos.alt.deg)
plt.axhline(-18, color='k')
plt.xlabel('Hours from 6pm AZ time')
plt.ylabel('Sun altitude')

Looks like it's just below 18 degrees at 7, so you should be good to go!


Try to actually compute to some arbitrary precision (rather than eye-balling on a plot) when 18 degree twilight or sunrise/sunset hits on that night.

In [37]:

Try converting the HCG7 coordinates to an equatorial frame at some other equinox a while in the past (like J2000). Do you see the precession of the equinoxes?

Hint: To see a diagram of the supported frames look here. One of those will do what you need if you give it the right frame attributes.

In [37]:


For lots more documentation on the many other features of astropy.coordinates, check out its section of the documentation.

You might also be interested in the astroplan affiliated package, which uses the astropy.coordinates to do more advanced versions of the tasks in the last section of this tutorial.