Here we go!
Hello folks. This notebook was prepared with the aim of enlightening you all with what's been done so far in my GSoC project: Re-implementation of sunpy.wcs as sunpy.coordinates. A lot has changed in the past few months until Astropy hit a somewhat stable period. We will be going through how to use the coordinates framework in its present form. Quick Recap: Astropy released their APE5 model for desigining coordinate frames, based on which the SunPy coordinates framework was envisaged. We have four frame classes so far -:
HelioGraphicStonyhurst
HelioGraphicCarrington
HelioCentric
HelioProjective
Each one of them subclasses ~astropy.coordinates.BaseCoordinateFrame
, except the second, which subclasses from the first (to re-use some details in the constructor of the first). APE5 also defines how to implement a transformation framework for your frame classes, and we have that too. We shall be using the ~astropy.coordinates.SkyCoord
class as the high-level interface to our frames.
SkyCoord
provides an easy-to-use interface to pass in data to our frames. Note that once you have defined a SkyCoord
or frame object, they are immutable. Frames can have two kinds of attributes mainly - attributes which are relevant to the representation of the frame in coordinate space, and FrameAttribute instances.from sunpy.coordinates import frames
from astropy.coordinates import * # Sorry for * import, heh
from astropy import units as u
sc = SkyCoord(1*u.deg, 1*u.deg, 3*u.km, frame='heliographicstonyhurst', dateobs='2011/01/01T00:00:45')
print(sc)
print(sc.frame)
print(sc.hlon) # An attribute of the representation of frame.
print(sc.dateobs) # A FrameAttribute derived attribute of the frame.
<SkyCoord (HelioGraphicStonyhurst): B0=0.0 deg, dateobs=2011-01-01 00:00:45, L0=0.0 deg, hlon=1.0 deg, hlat=1.0 deg, rad=3.0 km> <HelioGraphicStonyhurst Coordinate: B0=0.0 deg, dateobs=2011-01-01 00:00:45, L0=0.0 deg, hlon=1.0 deg, hlat=1.0 deg, rad=3.0 km> 1d00m00s 2011-01-01 00:00:45
~astropy.coordinates.BaseRepresentation
object must be fed to it. For example, HelioGraphicStonyhurst
expects parameters to be in ~sunpy.coordinates.SphericalWrap180Representation
, which is an offshoot of SphericalRepresentation
that allows the longitude parameter to have a range of -180 to 180 degrees. Frame instances can, of course, be represented in different forms.sc2 = SkyCoord(CartesianRepresentation(1*u.km, 1*u.km, 1*u.km), frame='heliographicstonyhurst', dateobs='2011/01/01T00:00:45')
print(sc2) # note how it prints in default representation.
print(sc2.represent_as(CartesianRepresentation)) # close to our original parameters.
print(sc2.represent_as(CylindricalRepresentation)) # why? because I can.
sc_setrepr = SkyCoord(1*u.deg, 1*u.deg, frame='heliographicstonyhurst', dateobs='2011/01/01T00:00:45')
sc_setrepr.representation = 'cartesian' # this changes the representation in which the SkyCoord object is printed, but does *not* change the underlying default rep.
# Doing sc.x/sc.y/sc.z would yield an error.
print(sc_setrepr)
<SkyCoord (HelioGraphicStonyhurst): B0=0.0 deg, dateobs=2011-01-01 00:00:45, L0=0.0 deg, hlon=45.0 deg, hlat=35.2643896828 deg, rad=1.73205080757 km> (1.0000000000000002, 1.0, 1.0) km (1.41421356237 km, 0.785398163397 rad, 1.0 km) <SkyCoord (HelioGraphicStonyhurst): B0=0.0 deg, dateobs=2011-01-01 00:00:45, L0=0.0 deg, x=695296.157659 km, y=12136.4395763 km, z=12138.2882964 km>
DynamicMatrixTransform
, StaticMatrixTransform
or FunctionTransform
. One may go through the Astropy Docs for more information.print(sc.transform_to('heliographiccarrington'))
print(sc.transform_to('heliocentric'))
print(sc.transform_to('helioprojective'))
print(sc.transform_to('heliographicstonyhurst'))
import numpy.testing as npt
sc3 = sc.transform_to('heliocentric').transform_to('helioprojective').transform_to('heliographicstonyhurst')
print('\n')
print((sc.hlon, sc3.hlon))
print((sc.hlat, sc3.hlat))
print((sc.rad, sc3.rad))
# Values are close.
<SkyCoord (HelioGraphicCarrington): B0=0.0 deg, dateobs=2011-01-01 00:00:45, L0=0.0 deg, hlon=-104.698696917 deg, hlat=1.0 deg, rad=3.0 km> <SkyCoord (HelioCentric): D0=149597870.7 km, B0=0.0 deg, dateobs=2011-01-01 00:00:45, L0=0.0 deg, x=0.0523492450538 km, y=0.0523572193119 km, z=2.99908624053 km> <SkyCoord (HelioProjective): L0=0.0 deg, B0=0.0 deg, dateobs=2011-01-01 00:00:45, D0=149597870.7 km, Tx=7.21788823206e-05 arcsec, Ty=7.21898771885e-05 arcsec, distance=149597867.701 km> <SkyCoord (HelioGraphicStonyhurst): B0=0.0 deg, dateobs=2011-01-01 00:00:45, L0=0.0 deg, hlon=1.0 deg, hlat=1.0 deg, rad=3.0 km> (<Longitude180 1.0 deg>, <Longitude180 0.9999999992842148 deg>) (<Latitude 1.0 deg>, <Latitude 0.9999999992844331 deg>) (<Distance 3.0 km>, <Distance 3.000000002146484 km>)
HelioGraphicStonyhurst
defaults its rad
representation attribute to the solar radius if no value is specified. Similarly, its L0
and B0
FrameAttributes are defaulted to zero degrees. HelioProjective
frames also make use of this defaulting mechanism for their distance
attributes. Specifically, HelioProjective
frames can also carry a zeta
property by which distance
is calculated. Only one of them is specified at a time.sc = SkyCoord(1*u.deg, 1*u.deg, frame='heliographicstonyhurst', dateobs='2011/01/01T00:00:45')
print(sc) # note how rad and L0,B0 magically appear.
from sunpy.coordinates.frames import HelioProjective
frame = SkyCoord(HelioProjective(1*u.deg, 1*u.deg, zeta=1*u.km, dateobs='2012/02/02T00:00:50'))
print(frame) # unfortunately, as zeta is a frame property object, SkyCoord cannot recognize it firsthand.
<SkyCoord (HelioGraphicStonyhurst): B0=0.0 deg, dateobs=2011-01-01 00:00:45, L0=0.0 deg, hlon=1.0 deg, hlat=1.0 deg, rad=695508.0 km> <SkyCoord (HelioProjective): L0=0.0 deg, B0=0.0 deg, dateobs=2012-02-02 00:00:50, D0=149597870.7 km, Tx=3600.0 arcsec, Ty=3600.0 arcsec, distance=149597869.7 km>
SkyCoord
has some other interesting methods. They are on display here. Also, see the documentation for more information.sc = SkyCoord([1,2,3]*u.deg, [3,2,1]*u.deg, [1,1,1]*u.km, frame='heliographicstonyhurst', dateobs='2011/01/01T00:00:45') # cool array type inputs
print(sc)
print(sc2.position_angle(sc3)) # calculates the position angle between this SkyCoord and another.
print(sc.separation(sc3)) # calculates on-sky separation. Extra keyword args such as 'dateobs' must be the same in both SkyCoords.
print(sc.separation_3d(sc3)) # calculates the 3D separation. Can not be used on 2D frames. Keywordargs condition applies.
<SkyCoord (HelioGraphicStonyhurst): B0=0.0 deg, dateobs=2011-01-01 00:00:45, L0=0.0 deg, (hlon, hlat, rad) in (deg, deg, km) [(1.0, 3.0, 1.0), (2.0, 2.0, 1.0), (3.0, 1.0, 1.0)]> 4.18879rad [u'2d00m00s' u'1d24m50.2642s' u'1d59m58.9033s'] [ 2.00091355 2.00045669 2.00091327] km