%matplotlib inline

Astropy Coordinates and SunPy Demo#

Written by Steven Christe and presented at the Heliopython meeting on November 13-15, 2018. The purpose of this demo is to show off the AstroPy coordinates framework as well as show how SunPy extends it to add solar coordinate systems.

The astropy coordinates package provides classes for representing a variety of celestial/spatial coordinates and their velocity components, as well as tools for converting between common coordinate systems in a uniform way.

from astropy import units as u
from astropy.coordinates import SkyCoord, AltAz
from astropy.time import Time

SkyCoord As an example of creating a SkyCoord to represent an ICRS (Right ascension [RA], Declination [Dec]) sky position:

c = SkyCoord(ra=10.625*u.degree, dec=41.2*u.degree, frame='icrs')

It can also handle arrays (many ways to instantiate a SkyCoord)

c = SkyCoord(ra=[10, 11, 12, 13]*u.degree, dec=[41, -5, 42, 0]*u.degree)

SkyCoord can also handle 3D positions, just add a distance

c = SkyCoord(ra=10.68458*u.degree, dec=41.26917*u.degree, distance=770*u.kpc)

So now cartesian coordinates are available

print('r = ({0}, {1}, {2})'.format(c.cartesian.x, c.cartesian.y, c.cartesian.z))
r = (568.7128654235231 kpc, 107.3008974042025 kpc, 507.88994291875713 kpc)

Positions of objects Can also register positions of objects or do object lookups

crab = SkyCoord.from_name("Crab")
<SkyCoord (ICRS): (ra, dec) in deg
    (83.63308333, 22.0145)>

Let’s now consider a position in the sky from a specific location on the Earth.

from astropy.coordinates import EarthLocation

Many positions are already availabe such as that of the VLA.

vla_coord = EarthLocation.of_site('vla')
(-1601184.40191992, -5041989.95569235, 3554875.07685646) m

Using a position on the Earth we can calculate Alt/Az, since dkist is missing from the library, let’s add it as a position

dkist = EarthLocation(lat=20.70818*u.deg, lon=-156.2569*u.deg, height=3084*u.m)
utcoffset = -10 * u.hour
midnight = Time('2018-11-14 00:00:00') - utcoffset

We can now get the position of the Crab in the sky as observed from DKIST

crab_altaz = crab.transform_to(AltAz(obstime=midnight,location=dkist))
print("Crab's Altitude = {0.alt:}".format(crab_altaz))
<SkyCoord (AltAz: obstime=2018-11-14 10:00:00.000, location=(-5466027.73432422, -2404324.10015092, 2242293.2433644) m, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron): (az, alt) in deg
    (80.97673524, 55.86134309)>
Crab's Altitude = 55.861343085319625 deg

Let’s now move on to showing how SunPy extends AstroPy coordinates by adding solar coordinate systems.

from sunpy.coordinates import frames
from sunpy.coordinates.sun import earth_distance

SunPy defines HeliographicStonyhurst, HeliographicCarrington, Heliocentric, and Helioprojective. Let’s define the center of the Sun

sun = SkyCoord(0 * u.arcsec, 0 * u.arcsec, obstime=midnight, frame=frames.Helioprojective, observer='Earth')

The position in the sky from the DKIST site is

sun_altaz = sun.transform_to(AltAz(obstime=midnight,location=dkist))
print('Altitude is {0} and Azimuth is {1}'.format(sun_altaz.T.alt, sun_altaz.T.az))
Altitude is -86.69112108195809 deg and Azimuth is 317.3715108577912 deg

As expected the Sun is below the horizon! Let’s consider noon now.

noon = Time('2018-11-14 12:00:00') - utcoffset
sun_altaz = sun.transform_to(AltAz(obstime=noon,location=dkist))
print('Altitude is {0} and Azimuth is {1}'.format(sun_altaz.T.alt, sun_altaz.T.az))
Altitude is 50.83292620533457 deg and Azimuth is 176.42009890469524 deg

Next let’s check this calculation by converting it back to helioprojective. We should get our original input which was the center of the Sun. To go from Altitude/Azimuth to Helioprojective, you will need the distance to the Sun. solar distance. Define distance with SunPy’s almanac.

distance = earth_distance(noon)
b = SkyCoord(az=sun_altaz.T.az, alt=sun_altaz.T.alt, distance=distance, frame=AltAz(obstime=noon,location=dkist), observer='Earth')
sun_helioproj = b.transform_to(frames.Helioprojective)
print('The helioprojective point is {0}, {1}'.format(sun_helioproj.T.Tx, sun_helioproj.T.Ty))
The helioprojective point is -9.287273101974279 arcsec, 1.0377010641434536 arcsec

Let’s now show off how we can convert between Solar coordinates Coordinates. Transform to HeliographicStonyhurst

<SkyCoord (HeliographicStonyhurst: obstime=2018-11-14 10:00:00.000, rsun=695700.0 km): (lon, lat, radius) in (deg, deg, AU)
    (0., 2.94605581, 0.00465047)>

Transform to Heliocentric

<SkyCoord (Heliocentric: obstime=2018-11-14 10:00:00.000, observer=<HeliographicStonyhurst Coordinate for 'Earth'>): (x, y, z) in AU
    (0., 0., 0.00465047)>

Transform to HeliographicCarrington

<SkyCoord (HeliographicCarrington: obstime=2018-11-14 10:00:00.000, rsun=695700.0 km, observer=<HeliographicStonyhurst Coordinate for 'Earth'>): (lon, lat, radius) in (deg, deg, AU)
    (115.41667055, 2.94605581, 0.00465047)>