Coords 2: Transforming between coordinate systems

Authors

Erik Tollerud, Kelle Cruz, Stephen Pardy, Stephanie T. Douglas

Learning Goals

• Create astropy.coordinates.SkyCoord objects

• Transform to different coordinate systems on the sky

• Transform to altitude/azimuth coordinates from a specific observing site

Keywords

coordinates, units, observational astronomy

Summary

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.

Imports

# Third-party dependencies
from astropy import units as u
from astropy.coordinates import SkyCoord
import numpy as np

# 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

Section 0: Quickstart

Note: If you already worked through Coords 1 you can feel free to skip to Section 1.

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:

hcg7_center = SkyCoord(9.81625*u.deg, 0.88806*u.deg, frame='icrs')  # using degrees directly
print(hcg7_center)

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

hcg7_center = SkyCoord('0h39m15.9s', '0d53m17.016s', frame='icrs')  # passing in string format
print(hcg7_center)

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

We can get the right ascension and declination components of the object directly by accessing those attributes.

print(hcg7_center.ra)
print(hcg7_center.dec)

9d48m58.5s
0d53m17.016s

Section 1:

Introducing frame transformations

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:

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.

hcg7_center.galactic

<SkyCoord (Galactic): (l, b) in deg
    (116.47556813, -61.83099472)>

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() <http://docs.astropy.org/en/stable/api/astropy.coordinates.SkyCoord.html#astropy.coordinates.SkyCoord.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.

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.

hcg7_center.transform_to('galactic')

We can transform to many coordinate frames and equinoxes.

These coordinates are available by default:

• ICRS

• FK5

• FK4

• FK4NoETerms

• Galactic

• Galactocentric

• Supergalactic

• AltAz

• GCRS

• CIRS

• ITRS

• HCRS

• PrecessedGeocentric

• GeocentricTrueEcliptic

• BarycentricTrueEcliptic

• HeliocentricTrueEcliptic

• SkyOffsetFrame

• GalacticLSR

• LSR

• BaseEclipticFrame

• BaseRADecFrame

Let’s focus on just a few of these. We can try FK5 coordinates next:

hcg7_center_fk5 = hcg7_center.transform_to('fk5')
print(hcg7_center_fk5)

<SkyCoord (FK5: equinox=J2000.000): (ra, dec) in deg
    (9.81625645, 0.88806155)>

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.

from astropy.coordinates import FK5
hcg7_center_fk5.transform_to(FK5(equinox='J1975'))  # precess to a different equinox

Beware: Changing frames also changes some of the attributes of the object, but usually in a way that makes sense. The following code should fail.

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:

print(hcg7_center.galactic.l, hcg7_center.galactic.b)

116d28m32.0453s -61d49m51.581s

Section 2:

Transform frames to get to altitude-azimuth (“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 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:

from astropy.coordinates import EarthLocation
from astropy.time import Time

Let’s first see the sky position at Kitt Peak National Observatory in Arizona.

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

kitt_peak = EarthLocation.of_site('Kitt Peak')

We can see the list of observing sites:

EarthLocation.get_site_names()

['',
 '',
 '',
 'ALMA',
 'ATST',
 'Anglo-Australian Observatory',
 'Apache Point',
 'Apache Point Observatory',
 'Atacama Large Millimeter Array',
 'BAO',
 'BBSO',
 'Beijing XingLong Observatory',
 'Black Moshannon Observatory',
 'CHARA',
 'Canada-France-Hawaii Telescope',
 'Catalina Observatory',
 'Cerro Pachon',
 'Cerro Paranal',
 'Cerro Tololo',
 'Cerro Tololo Interamerican Observatory',
 'DCT',
 'DKIST',
 'Discovery Channel Telescope',
 'Dominion Astrophysical Observatory',
 'GBT',
 'Gemini South',
 'Green Bank Telescope',
 'Hale Telescope',
 'Haleakala Observatories',
 'Happy Jack',
 'IAO',
 'JCMT',
 'James Clerk Maxwell Telescope',
 'Jansky Very Large Array',
 'Keck Observatory',
 'Kitt Peak',
 'Kitt Peak National Observatory',
 'La Silla Observatory',
 'Large Binocular Telescope',
 'Las Campanas Observatory',
 'Lick Observatory',
 'Lowell Observatory',
 'MWA',
 'Manastash Ridge Observatory',
 'McDonald Observatory',
 'Medicina',
 'Medicina Dish',
 'Michigan-Dartmouth-MIT Observatory',
 'Mount Graham International Observatory',
 'Mt Graham',
 'Mt. Ekar 182 cm. Telescope',
 'Mt. Stromlo Observatory',
 'Multiple Mirror Telescope',
 'Murchison Widefield Array',
 'NOV',
 'NST',
 'National Observatory of Venezuela',
 'Noto',
 'Observatorio Astronomico Nacional, San Pedro Martir',
 'Observatorio Astronomico Nacional, Tonantzintla',
 'Palomar',
 'Paranal Observatory',
 'Roque de los Muchachos',
 'SAAO',
 'SALT',
 'SPO',
 'SRT',
 'Sac Peak',
 'Sacramento Peak',
 'Siding Spring Observatory',
 'Southern African Large Telescope',
 'Subaru',
 'Subaru Telescope',
 'Sunspot',
 'Sutherland',
 'TUG',
 'UKIRT',
 'United Kingdom Infrared Telescope',
 'Vainu Bappu Observatory',
 'Very Large Array',
 'W. M. Keck Observatory',
 'Whipple',
 'Whipple Observatory',
 'aao',
 'alma',
 'apo',
 'bbso',
 'bmo',
 'cfht',
 'ctio',
 'dao',
 'dct',
 'dkist',
 'ekar',
 'example_site',
 'flwo',
 'gbt',
 'gemini_north',
 'gemini_south',
 'gemn',
 'gems',
 'greenwich',
 'haleakala',
 'iao',
 'irtf',
 'jcmt',
 'keck',
 'kpno',
 'lapalma',
 'lasilla',
 'lbt',
 'lco',
 'lick',
 'lowell',
 'mcdonald',
 'mdm',
 'medicina',
 'mmt',
 'mro',
 'mso',
 'mtbigelow',
 'mwa',
 'mwo',
 'noto',
 'ohp',
 'paranal',
 'salt',
 'sirene',
 'spm',
 'spo',
 'srt',
 'sso',
 'tona',
 'tug',
 'ukirt',
 'vbo',
 'vla']

Let’s check the altitude at 1 AM UTC, which is 6 PM AZ mountain time:

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

from astropy.coordinates import AltAz

aa = AltAz(location=kitt_peak, obstime=observing_time)
print(aa)

<AltAz Frame (obstime=2010-12-21 01:00:00.000, location=(-1994310.09211632, -5037908.60633759, 3357621.75212217) m, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron)>

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

hcg7_center.transform_to(aa)

<SkyCoord (AltAz: obstime=2010-12-21 01:00:00.000, location=(-1994310.09211632, -5037908.60633759, 3357621.75212217) m, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron): (az, alt) in deg
    (149.19392036, 55.0624736)>

To look at just the altitude we can alt attribute:

hcg7_center.transform_to(aa).alt

$55^\circ03{}^\prime44.905{}^{\prime\prime}$

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.

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

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:

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

<SkyCoord (AltAz: obstime=2019-11-25 19:19:56.799292, location=(-1994502.60430614, -5037538.54232911, 3358104.99690298) m, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron): (az, alt) in deg
    (67.17087052, -30.274766)>

Exercises

Exercise 1

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

observing_time = Time('2000-12-21 1:00')

from astropy.coordinates import AltAz

aa = AltAz(location=kitt_peak, obstime=observing_time)
print(aa)

<AltAz Frame (obstime=2000-12-21 01:00:00.000, location=(-1994310.09211632, -5037908.60633759, 3357621.75212217) m, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron)>

hcg7_center.transform_to(aa)

<SkyCoord (AltAz: obstime=2000-12-21 01:00:00.000, location=(-1994310.09211632, -5037908.60633759, 3357621.75212217) m, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron): (az, alt) in deg
    (150.09970846, 55.24312468)>

hcg7_center.transform_to(aa).alt

$55^\circ14{}^\prime35.2488{}^{\prime\prime}$

type(hcg7_center.transform_to(aa).alt)

astropy.coordinates.angles.Latitude

Precession data:

alt 21-12-2010: \(55°03′24.2307′′\) vs \(55°14′14.5696′′\)

from astropy.coordinates import Angle
precess10yr= hcg7_center.transform_to(aa).alt- Angle("55:03:24.2307", unit="degree")

precess10yr

$0^\circ11{}^\prime11.0181{}^{\prime\prime}$

rondgangsectie=precess10yr/Angle("1:00:00", unit="degree")

rondgangsectie

$0.18639393 \; \mathrm{}$

rondgang=180*10*rondgangsectie

rondgang

$335.50907 \; \mathrm{}$

Exercise 2

degree1=Angle("1:00:00", unit="degree")

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.

Exercise 3

Try looking at the altitude of HCG 7 at another observatory.

Wrap-up

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.