Erik Tollerud, Kelle Cruz, Stephen Pardy, Stephanie T. Douglas
astropy.coordinates.SkyCoord objectscoordinates, units, observational astronomy
In this tutorial we demonstrate how to define astronomical coordinates using the astropy.coordinates "frame" classes. We then show how to transform between the different built-in coordinate frames, such as from ICRS (RA, Dec) to Galactic (l, b). Finally, we show how to compute altitude and azimuth from a specific observing site.
In [ ]:
# Third-party dependencies
from astropy import units as u
from astropy.coordinates import SkyCoord
import numpy as np
In [ ]:
# Set up matplotlib and use a nicer set of plot parameters
from astropy.visualization import astropy_mpl_style
import matplotlib.pyplot as plt
plt.style.use(astropy_mpl_style)
%matplotlib inline
In Astropy, the most common object you'll work with for coordinates is SkyCoord. A SkyCoord can most easily be created directly from angles as shown below.
In this tutorial we'll be converting between frames. Let's start in the ICRS frame (which happens to be the default.)
For much of this tutorial we'll work with the Hickson Compact Group 7. We can create an object either by passing the degrees explicitly (using the astropy units library) or by passing in strings. The two coordinates below are equivalent:
In [ ]:
hcg7_center = SkyCoord(9.81625*u.deg, 0.88806*u.deg, frame='icrs') # using degrees directly
print(hcg7_center)
In [ ]:
hcg7_center = SkyCoord('0h39m15.9s', '0d53m17.016s', frame='icrs') # passing in string format
print(hcg7_center)
We can get the right ascension and declination components of the object directly by accessing those attributes.
In [ ]:
print(hcg7_center.ra)
print(hcg7_center.dec)
astropy.coordinates provides many tools to transform between different coordinate systems. For instance, we can use it to transform from ICRS coordinates (in RA and Dec) to Galactic coordinates.
To understand the code in this section, it may help to read over the overview of the astropy coordinates scheme. The key piece to understand is that all coordinates in Astropy are in particular "frames" and we can transform between a specific SkyCoord object in one frame to another. For example, we can transform our previously-defined center of HCG 7 from ICRS to Galactic coordinates:
In [ ]:
hcg7_center = SkyCoord(9.81625*u.deg, 0.88806*u.deg, frame='icrs')
There are three different ways of transforming coordinates. Each has its pros and cons, but all should give you the same result. The first way to transform to other built-in frames is by specifying those attributes. For instance, let's see the location of HCG 7 in Galactic coordinates.
In [ ]:
hcg7_center.galactic
The above is actually a special "quick-access" form that internally does the same as what's in the cell below: it uses the transform_to() method to convert from one frame to another. We can pass in an empty coordinate class to specify what coordinate system to transform into.
In [ ]:
from astropy.coordinates import Galactic # new coordinate baseclass
hcg7_center.transform_to(Galactic())
Finally, we can transform using the transform_to() method and a string with the name of a built-in coordinate system.
In [ ]:
hcg7_center.transform_to('galactic')
We can transform to many coordinate frames and equinoxes.
These coordinates are available by default:
Let's focus on just a few of these. We can try FK5 coordinates next:
In [ ]:
hcg7_center_fk5 = hcg7_center.transform_to('fk5')
print(hcg7_center_fk5)
And, as with the Galactic coordinates, we can acheive the same result by importing the FK5 class from the astropy.coordinates package. This also allows us to change the equinox.
In [ ]:
from astropy.coordinates import FK5
hcg7_center_fk5.transform_to(FK5(equinox='J1975')) # precess to a different equinox
In [ ]:
hcg7_center.galactic.ra # should fail because Galactic coordinates are l/b not RA/Dec
Instead, we now have access to the l and b attributes:
In [ ]:
print(hcg7_center.galactic.l, hcg7_center.galactic.b)
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 altitude-azimuth coordinates, or "AltAz". We first need to specify both where and when we want to try to observe.
We'll need to import a few more specific modules:
In [ ]:
from astropy.coordinates import EarthLocation
from astropy.time import Time
Let's first see the sky position at Kitt Peak National Observatory in Arizona.
In [ ]:
# Kitt Peak, Arizona
kitt_peak = EarthLocation(lat='31d57.5m', lon='-111d35.8m', height=2096*u.m)
For known observing sites we can enter the name directly.
In [ ]:
kitt_peak = EarthLocation.of_site('Kitt Peak')
We can see the list of observing sites:
In [ ]:
EarthLocation.get_site_names()
Let's check the altitude at 1 AM UTC, which is 6 PM AZ mountain time:
In [ ]:
observing_time = Time('2010-12-21 1:00')
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 [ ]:
from astropy.coordinates import AltAz
aa = AltAz(location=kitt_peak, obstime=observing_time)
print(aa)
Now we can transform our ICRS SkyCoord to AltAz to get the location in the sky over Kitt Peak at the requested time.
In [ ]:
hcg7_center.transform_to(aa)
To look at just the altitude we can alt attribute:
In [ ]:
hcg7_center.transform_to(aa).alt
Alright, it's at 55 degrees at 6 PM, but that's pretty early to be observing. We could try various times one at a time to see if the airmass is at a darker time, but we can do better: let's try to create an airmass plot.
In [ ]:
# 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=kitt_peak, 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()
Great! Looks like the lowest airmass is in another hour or so (7 PM). But 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. Let's use it to check if we're in 18-degree twilight or not.
In [ ]:
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()
Looks like it's just below 18 degrees at 7 PM, so you should be good to go!
We can also look at the object altitude at the present time and date:
In [ ]:
now = Time.now()
hcg7_center = SkyCoord(9.81625*u.deg, 0.88806*u.deg, frame='icrs')
kitt_peak_aa = AltAz(location=kitt_peak, obstime=now)
print(hcg7_center.transform_to(kitt_peak_aa))
Try to compute to some arbitrary precision (rather than eyeballing on a plot) when 18 degree twilight or sunrise/sunset hits on that night.
In [ ]:
Try converting the HCG 7 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 or on the list above. One of those will do what you need if you give it the right frame attributes.
In [ ]:
Try looking at the altitude of HCG 7 at another observatory.
In [ ]:
For 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.