Converting helioprojective to equatorial
Table of Contents
Introduction
Solar physics and nighttime astronomy often use different but related coordinate systems. In solar physics, the helioprojective coordinate system is common. Depending on the subfield of astronomy one might use a variety of coordinate systems, but one of the most common is equatorial. These are both spherical geometries.
However, it turns out to be a bit complicated to convert directly between the two. Astropy and SunPy offer utilities to do the conversion, but I was running into some challenges. I asked Sam Van Kooten for help, and we came up with this solution. It only applies when the observer is the center of the Earth though.
Converting a single coordinate
First, we need to know how to convert a single coordinate.
Astropy has a system of transormations that form a computational graph. SunPy augments this with the common solar coordinate systems. If we try to convert helioprojective to GCRS, we'll be doing several steps of conversions: Helioprojective to Heliocentric to HeliographicStonyhurst to HCRS to ICRS and finally to GCRS. This shouldn't really be required. So, we'll register a new transformation that goes from helioprojective straight to GCRS.
# Ripped from sunpy, modified
# Checks whether times are equal
# We explicitly perform the check in TAI to avoid possible numerical precision differences
# between a time in UTC and the same time after a UTC->TAI->UTC conversion
return
# We also deem the times equal if one is None
return is None or is None
"""
Get the P angle.
Return the position (P) angle for the Sun at a specified time, which is the angle between
geocentric north and solar north as seen from Earth, measured eastward from geocentric north.
The range of P is +/-26.3 degrees.
Parameters
----------
time : {parse_time_types}
Time to use in a parse_time-compatible format
Returns
-------
out : `~astropy.coordinates.Angle`
The position angle
"""
=
# Define the frame where its Z axis is aligned with geocentric north
=
return # noqa: SLF001
= or
# Compute the three angles we need
=
, = ,
= .
# Prepare rotation matrices for each
=
=
=
# Compose the matrices
= @ @
# Extract the input coordinates and negate the HP latitude,
# since it increases in the opposite direction from RA.
=
=
# Apply the tranformation
=
=
# Match the input representation. (If the input was UnitSpherical, meaning there's no
# distance coordinate, this drops the distance coordinate.)
=
# Put the computed coordinates into the output frame
return
The meat of the transformation happens when setting up the matrices in the last function.
They're rotation matrices to go from one to the other using the Sun's right ascension and declination coupled with the solar P angle.
We carefully use the P angle in GCRS instead of SunPy's built-in sun.P, which I believe uses ITRS.
We can then use Astropy's infrastructure for conversions.
=
=
It's quite nifty.
Transforming a WCS from helio to celestial
So let's say you have a world coordinate system in helioprojective. How can you use this to convert to GCRS?
First, you transform the CRVAL using the coordinate transform we just defined. Then, you transform the PC matrix by rotating by the P angle.
"""Extract CROTA from a WCS."""
return *
"""
Calculate a PC matrix given CROTA and CDELT.
Parameters
----------
crota : float
rotation angle from the WCS
cdelt : float
pixel size from the WCS
Returns
-------
np.ndarray
PC matrix
"""
return
"""Calculate the celestial WCS from a helio WCS."""
= == 3
=
=
= . - .
=
= *
= *
=
=
=
return
Since this code was developed for PUNCH, we have a polarization axis that makes this WCS 3D in some cases. So, we handle that here too.
Transforming the opposite direction
You can reverse this too!
# noqa: ANN201, N803, ANN001
"""Convert GCRS to HPC."""
# noqa: TRY003, EM101
= or
# Compute the three angles we need
=
, = ,
= .
# Prepare rotation matrices for each
=
=
=
# Compose the matrices
= @ @
=
# Extract the input coordinates and negate the HP latitude,
# since it increases in the opposite direction from RA.
= # , earth_distance(obstime))
=
# Apply the transformation
=
=
# Match the input representation. (If the input was UnitSpherical, meaning there's no
# distance coordinate, this drops the distance coordinate.)
=
= # , earth_distance(obstime))
=
# Put the computed coordinates into the output frame
return
"""Calculate the helio WCS from a celestial WCS."""
= == 3
# we're at the center of the Earth
=
=
# follow the SunPy tutorial from here
# https://docs.sunpy.org/en/stable/generated/gallery/units_and_coordinates/radec_to_hpc_map.html#sphx-glr-generated-gallery-units-and-coordinates-radec-to-hpc-map-py
=
=
=
=
=
=
= . + .
=
=
=
return ,
Verifying
We have many tests to verify this works properly. We use the Sun's position to go back and forth with known coordinates. We also test that given an arbitrary WCS you can convert to the other system and then convert back and get the original input. All our tests pass.
Conclusions
This code was developed for PUNCH to go back and forth between helioprojective and GCRS. It was used to create some synthetic images that are like what PUNCH will see. Those files include both a helioprojective and equatorial WCS so you can locate things in context of the Sun and the stars.
Here's a sneak peak of what these images look like.
I'll write more about the synthetic images in the future.