Note
Click here to download the full example code or to run this example in your browser via Binder
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.
SkyCoord As an example of creating a SkyCoord to represent an ICRS (Right ascension [RA], Declination [Dec]) sky position:
It can also handle arrays (many ways to instantiate a SkyCoord)
SkyCoord can also handle 3D positions, just add a distance
So now cartesian coordinates are available
print('r = ({0}, {1}, {2})'.format(c.cartesian.x, c.cartesian.y, c.cartesian.z))
Out:
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")
print(crab)
Out:
<SkyCoord (ICRS): (ra, dec) in deg
(83.63308333, 22.0145)>
let’s consider 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')
print(vla_coord)
Out:
(-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
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_altaz)
print("Crab's Altitude = {0.alt:}".format(crab_altaz))
Out:
<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
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))
Out:
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))
Out:
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))
Out:
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
Out:
<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
Out:
<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
Out:
<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)>
Total running time of the script: ( 0 minutes 2.177 seconds)