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

In[1]:

# Python standard-library
from urllib.parse import urlencode
from urllib.request import urlretrieve

# Third-party dependencies
from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.table import Table
import numpy as np
from IPython.display import Image


In[2]:

# Set up matplotlib and use a nicer set of plot parameters
import matplotlib.pyplot as plt
%matplotlib inline


## 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[3]:

hcg7_center = SkyCoord(9.81625*u.deg, 0.88806*u.deg, frame='icrs')
hcg7_center


Out[3]:

<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[4]:

SkyCoord('0h39m15.9s', '0d53m17.016s', frame='icrs')


Out[4]:

<SkyCoord (ICRS): (ra, dec) in deg
(9.81625, 0.88806)>


In[5]:

SkyCoord('0:39:15.9 0:53:17.016', unit=(u.hour, u.deg), frame='icrs')


Out[5]:

<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.

If 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 <https://julien.danjou.info/blog/2013/guide-python-static-class-abstract-methods>__.

In[6]:

hcg7_center = SkyCoord.from_name('HCG 7')
hcg7_center


Out[6]:

<SkyCoord (ICRS): (ra, dec) in deg
(9.81625, 0.88806)>


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 <http://docs.astropy.org/en/stable/units/index.html>__ objects (actually, a subclass called Angle <http://docs.astropy.org/en/stable/api/astropy.coordinates.Angle.html>__, which in turn is subclassed by Latitude <http://docs.astropy.org/en/stable/api/astropy.coordinates.Latitude.html>__ and Longitude <http://docs.astropy.org/en/stable/api/astropy.coordinates.Longitude.html>__). 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[7]:

type(hcg7_center.ra), type(hcg7_center.dec)


Out[7]:

(astropy.coordinates.angles.Longitude, astropy.coordinates.angles.Latitude)


In[8]:

hcg7_center.dec


Out[8]:

$0^\circ53{}^\prime17.016{}^{\prime\prime}$

In[9]:

hcg7_center.ra


Out[9]:

$9^\circ48{}^\prime58.5{}^{\prime\prime}$

In[10]:

hcg7_center.ra.hour


Out[10]:

0.6544166666666668


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[11]:

impix = 1024
imsize = 12*u.arcmin
cutoutbaseurl = 'http://skyservice.pha.jhu.edu/DR12/ImgCutout/getjpeg.aspx'
query_string = urlencode(dict(ra=hcg7_center.ra.deg,
dec=hcg7_center.dec.deg,
width=impix, height=impix,
scale=imsize.to(u.arcsec).value/impix))
url = cutoutbaseurl + '?' + query_string

urlretrieve(url, 'HCG7_SDSS_cutout.jpg')


Out[11]:

('HCG7_SDSS_cutout.jpg', <http.client.HTTPMessage at 0x7ff5cf81f048>)


Now lets take a look at the image.

In[12]:

Image('HCG7_SDSS_cutout.jpg')


Out[12]:

Very pretty!

### Exercises¶

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[None]:

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 http://archive.stsci.edu/cgi-bin/dss_search?f=GIF&ra=RA&dec=DEC)

In[None]:

## 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[13]:

from astroquery.sdss import SDSS
spectro=True,
photoobj_fields=['ra','dec','u','g','r','i','z'])


Out[13]:

/home/circleci/project/venv/lib/python3.6/site-packages/astroquery/sdss/__init__.py:29: UserWarning: Experimental: SDSS has not yet been refactored to have its API match the rest of astroquery (but it's nearly there).
warnings.warn("Experimental: SDSS has not yet been refactored to have its API "


astroquery queries gives us back an astropy.table.Table object <http://docs.astropy.org/en/stable/table/index.html>__. 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 Table.read and use the copy of this table included in the tutorial.)

In[14]:

sdss.write('HCG7_SDSS_photo.dat', format='ascii')


Out[14]:

WARNING: AstropyDeprecationWarning: HCG7_SDSS_photo.dat already exists. Automatically overwriting ASCII files is deprecated. Use the argument 'overwrite=True' in the future. [astropy.io.ascii.ui]


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[15]:

sdss = Table.read('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[16]:

twomass = Table.read('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[17]:

sdss # just to see an example of the format


Out[17]:

Table length=679
float64float64float64float64float64float64float64
9.80588680002290.86416190139736914.9643813.2688812.3802911.9426711.59935
9.491213075462920.74076541700878321.943619.6844518.2597817.5334117.07395
9.814441447426250.94258449363529619.5578417.8927817.2430816.9455216.73516
10.0778276825780.92386587860973619.6258319.3954119.2246619.0558918.77663
9.511939969534981.0857385368605722.6077521.8582520.0075519.1001718.70797
9.550194344498781.01915084951718.7004217.8445917.5772717.4841517.41936
10.00592372609530.98576150885775619.9912519.4045719.1993319.1178919.12818
9.837334475500820.84735815689209321.8971921.3275521.0155421.1413120.89685
10.10756829750370.80201251479586519.8583617.9077917.1151316.8146716.64346
9.759581991658110.7240297461704522.4751419.4800317.9986916.9916616.41765
.....................
9.921562046363480.6771381563455725.038522.8634621.6558520.9129620.6877
10.04822960184670.91653085245531220.3106419.3228518.8668818.5563218.40811
10.1358136613291.1797166320310822.0879520.1580818.4718217.8788517.51656
10.03215369554560.94156277443788222.3757120.0096818.7104918.0987617.74478
9.90021905419111.0916596634240618.2996716.9620516.5370616.4193216.36538
9.863046491255770.86467706305473122.7427221.7234621.0894620.8056220.52088
9.987558930446770.98487044260131819.7621718.8578118.4351718.0940218.0476
9.726191945520211.0649273515845221.9062419.3556317.971216.6404615.9429
10.00896441373031.0507495306401120.660218.9604918.1452217.7840617.53901
9.813631859453180.81751275690135123.4714620.9206319.21218.5175118.19212

In[18]:

twomass # just to see an example of the format


Out[18]:

degdegarcsecmagmagmagmagmagmagdegdegarcsecmagmagmagmagmagmagarcsecdeg
str16float64float64float64float64float64int64float64float64int64float64float64int64float64int64float64int64float64float64float64float64float64float64float64str1float64float64
00402069+005250810.0862180.8807989.413.8350.068013.010.086012.5880.08900.8700.823518.6213.6320.08812.7440.10412.3980.1050972.12061191.538952
00395984+01035459.999351.0651412.912.9250.035012.1830.042011.890.06700.8350.74035.912.4690.04811.910.06611.5220.0870916.92763645.951861
00401849+004944810.0770620.829136.014.9180.086014.1130.107013.7140.10300.6-151.09011.3514.6310.12113.9530.16913.5250.1610962.489231102.73149
00395277+00571249.9699070.9534725.314.7020.049014.2480.069013.8990.09500.6-600.44-5010.5914.620.14414.150.29613.730.20601.13644466.93659
00401864+004724510.0777040.7901437.615.5850.134115.0030.18114.0490.14210.5300.463014.4814.9770.13814.8550.30313.6530.1801004.982128110.53147
00393485+00513559.8952190.85988239.311.4150.031310.7550.044310.5140.06830.6-300.7-6092.2911.4150.01810.1550.0549.9760.0850301.813395109.639102
00392964+01034959.8735261.06376910.914.4630.065013.6180.067013.2580.09100.4550.286020.3514.20.08613.3630.09113.1010.1330665.30141518.051526
00403343+004907910.1392930.8188655.015.4840.150------13.970.13701.0901.09010.0515.0350.18314.7250.013.6540.18901189.207905102.088788
00393319+00355059.8883050.59738111.513.1560.033012.5090.043012.0730.05900.6-550.52-4021.6413.0260.0412.2470.04611.9780.06501078.11027166.0785
.................................................................................
00391798+00415889.8249360.6996876.115.6850.168014.890.191014.0030.15501.0901.09011.415.6770.31214.4150.22613.5680.190678.863209177.360117
00384796+00345729.6998580.5825785.114.9250.077014.2240.114013.5360.07901.0901.09010.214.8390.13314.1110.19213.4610.13701176.842625200.856597
00390392+00505799.7663450.8494195.014.8950.07014.2380.087013.8340.1101.0901.09010.0514.7060.10714.0330.13213.750.1870227.201453232.24689
00391339+00515089.8057970.86413552.810.3620.01409.6310.01709.3340.02400.3-150.4-1575.0210.2790.0159.5270.0169.2470.023093.990015203.598476
00391786+00544589.8244180.91274327.911.0820.016010.3840.022010.1470.03200.550.7542.7510.9140.01810.2510.02110.0310.03093.59655518.308033
00385879+00572699.7449710.9574785.015.5350.122014.7960.145014.2780.16501.0901.09010.0515.5350.12214.6230.22714.1470.2690358.163568314.246475
00391879+00533089.8283030.89190915.413.0440.047012.4120.063012.0770.09400.8600.746523.6212.7550.04812.2830.07211.7130.096045.54456272.287562
00391213+01024089.800551.0446915.015.5680.126015.0470.181014.3560.17601.0901.09010.0515.2950.18115.0470.18114.0670.250566.696375354.276982
00383990+01044429.6662681.0789685.315.2550.108014.2320.121013.8730.11301.0901.09010.4415.1510.1813.8120.14913.5520.1550873.946372321.851314
00384916+00502129.7048720.8392445.115.0750.088014.6510.17013.8040.10101.0901.09010.215.0530.15914.6510.1713.6820.1710437.740484246.331036

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[19]:

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[20]:

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[21]:

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


Out[21]:

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[22]:

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


Out[22]:

(<Angle [ 9.75040523,  9.94247716,  9.51932457, 10.10700162,  9.91149052,
10.07348473,  9.63040458,  9.72564952, 10.12905391, 10.13845705,
10.14555822, 10.14727304,  9.98444481,  9.86790027,  9.60942396,
9.81159962,  9.81474753,  9.93849517,  9.84686176,  9.81952013,
10.13719007,  9.77876948,  9.83794283] deg>,
<Angle [0.68425352, 0.75630516, 1.01135303, 0.80980643, 0.79894519,
0.89597708, 0.57416757, 0.84498122, 0.69965054, 1.19019674,
1.0687527 , 0.6926511 , 0.84684462, 1.16005717, 0.79050233,
1.19149118, 0.90539845, 0.88624344, 0.77435883, 0.69998442,
0.89731188, 0.5900144 , 1.10508309] 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[23]:

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[24]:

plt.hist(d2d_sim.arcsec, bins='auto', histtype='step', label='Simulated', linestyle='dashed')
plt.hist(d2d_sdss.arcsec, bins='auto', histtype='step', label='2MASS')
plt.xlabel('separation [arcsec]')
plt.legend(loc=0)
plt.tight_layout()


Out[24]:

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[25]:

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.xlabel('r-K')
plt.ylabel('r')
plt.xlim(2.5, 4)
plt.ylim(18, 12) #mags go backwards!

plt.subplot(1, 2, 2)
plt.scatter(rKcolor, rmag)
plt.xlabel('r-K')
plt.ylabel('g-r')
plt.xlim(2.5, 4)

plt.tight_layout()


Out[25]:

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.

### Exercises¶

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[None]:

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[None]:

## 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[26]:

hcg7_center.galactic


Out[26]:

<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[27]:

from astropy.coordinates import Galactic
hcg7_center.transform_to(Galactic())


Out[27]:

<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[28]:

hcg7_center.galactic.ra  # should fail because galactic coordinates are l/b not RA/Dec


Out[28]:

AttributeErrorTraceback (most recent call last)

<ipython-input-28-d7bc134707f6> in <module>()
----> 1 hcg7_center.galactic.ra  # should fail because galactic coordinates are l/b not RA/Dec

~/project/venv/lib/python3.6/site-packages/astropy/coordinates/sky_coordinate.py in __getattr__(self, attr)
693         # Fail
694         raise AttributeError("'{0}' object has no attribute '{1}'"
--> 695                              .format(self.__class__.__name__, attr))
696
697     def __setattr__(self, attr, val):

AttributeError: 'SkyCoord' object has no attribute 'ra'


In[29]:

hcg7_center.galactic.b


Out[29]:

$-61^\circ49{}^\prime51.581{}^{\prime\prime}$

### 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[30]:

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[31]:

from astropy.coordinates import AltAz

aa = AltAz(location=observing_location, obstime=observing_time)
aa


Out[31]:

<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[32]:

hcg7_center.transform_to(aa)


Out[32]:

<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.19392036, 55.0624736)>


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[33]:

# 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)]')
plt.ylim(0.9,3)
plt.tight_layout()


Out[33]:

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[34]:

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')
plt.tight_layout()


Out[34]:

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

### Exercises¶

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[None]:

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[None]:

## Wrap-up¶

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.