Transforming Coordinates Between SpacePy, Astropy, and SunPy

Written by Matt Wentzel-Long. The purpose of this example is to demonstrate how to pass coordinates between SpacePy, Astropy, and SunPy, and how to compute some simple transformations in each package.

from spacepy.coordinates import Coords
import spacepy.time as spt

import astropy.units as u
from astropy.coordinates import SkyCoord

from sunpy.coordinates import frames

First, create a SpacePy coordinate in the (Cartesian) Geographic Coordinate System (GEO) and attach an observation time to the coordinate. Units are in Earth radii (Re).

coord = Coords([[1,2,4],[1,2,2]], 'GEO', 'car')
coord.ticks = spt.Ticktock(['2002-02-02T12:00:00', '2002-02-02T12:00:00'], 'ISO')

print(coord)

Coords( [[1, 2, 4], [1, 2, 2]] , 'GEO', 'car')

In SpacePy, the convert method can be used to easily convert coordinates into one of the 10 coordinate systems supported. For example, convert the coordinates to the (Cartesian) Solar Magnetic system.

sm = coord.convert('SM','car')
print(sm)

Coords( [[0.814408577493688, 2.648338696371909, 3.6500740839336188], [0.9232104562461743, 2.3057302593531515, 1.6826438793104677]] , 'SM', 'car')

Send the coordinates to an Astropy SkyCoord instance using the SpacePy to_skycoord function. Units are converted to meters.

Note: this must be in the GEO system.

skycoord = coord.to_skycoord()
print(skycoord)

<SkyCoord (ITRS: obstime=[6.96686413e+08 6.96686413e+08]): (x, y, z) in m
    [(6371200., 12742400., 25484800.), (6371200., 12742400., 12742400.)]>

See the Astropy documentation for transforming coordinates. Here is a simple example that transforms the skycoord into the FK5 system.

sky_fk5 = skycoord.transform_to('fk5')
print(sky_fk5)

<SkyCoord (FK5: equinox=J2000.000): (ra, dec, distance) in (deg, deg, m)
    [(136.18233833, 16.72158771, 1.46960941e+11),
     (136.18233946, 16.71682971, 1.46957277e+11)]>

Use the SunPy frames function to transform this coordinate into a Heliogaphic Carrington coordinate.

Note: helioprojective frames require that an observer be defined.

sky_helio = skycoord.transform_to(frames.HeliographicCarrington(observer="earth"))
print(sky_helio)

<SkyCoord (HeliographicCarrington: obstime=[6.96686413e+08 6.96686413e+08], rsun=695700.0 km, observer=<HeliographicStonyhurst Coordinate (obstime=[6.96686413e+08 6.96686413e+08], rsun=695700.0 km): (lon, lat, radius) in (deg, deg, AU)
    [(0., -6.12373415, 0.98554015), (0., -6.12373415, 0.98554015)]>): (lon, lat, radius) in (deg, deg, m)
    [(11.75877738, -6.11476108, 1.47435249e+11),
     (11.75982715, -6.11938557, 1.47431569e+11)]>

See the Astropy Coordinates and SunPy Demo for coordinate transformations in SunPy (TODO: update this link once Jupyter Book gallery replaces the old gallery).

Now, convert the coordinate back into its original form to demonstrate transformations in the other direction, and the loss of precision. First, convert this back to GEO coordinates (ITRS in Astropy).

sun_geo = sky_helio.transform_to('itrs')
print(sun_geo)

<SkyCoord (ITRS: obstime=[6.96686413e+08 6.96686413e+08]): (x, y, z) in m
    [(6371199.99999981, 12742399.99995084, 25484799.99999662),
     (6371199.99999466, 12742400.00000276, 12742400.00001287)]>

Lastly, use the SpacePy from_skycoord function to transform this back into a SpacePy coordinate.

coord = Coords.from_skycoord(sun_geo)
print(coord)

Coords( [[0.9999999999999699, 1.9999999999922848, 3.999999999999469], [0.9999999999991617, 2.000000000000434, 2.0000000000020193]] , 'GEO', 'car')

The observation time is now in Astropy time (APT) (see here).

print(coord.ticks)

Ticktock( [6.96686413e+08 6.96686413e+08], dtype=APT)

You can verify that this is the original observation time by converting it to ISO.

print(coord.ticks.getISO())

['2002-02-02T12:00:00' '2002-02-02T12:00:00']