Reprojecting a Helioprojective Map to Helioprojective Radial

In this example we use reproject to transform a helioprojective Map to helioprojective radial.

import matplotlib.pyplot as plt

import astropy.units as u

import sunpy.data.sample
import sunpy.map
from sunpy.coordinates import HelioprojectiveRadial
from sunpy.map.header_helper import make_hpr_header
from sunpy.visualization import show_hpr_impact_angle

Let’s start with an AIA map in the helioprojective Cartesian (Helioprojective) frame.

aia_map = sunpy.map.Map(sunpy.data.sample.AIA_193_IMAGE)

We now plot the map and overlay grid lines for position angle in the helioprojective radial (HelioprojectiveRadial) frame.

fig = plt.figure()
ax = fig.add_subplot(projection=aia_map)
aia_map.plot(axes=ax)

overlay = ax.get_coords_overlay(
    HelioprojectiveRadial(
        obstime=aia_map.reference_date,
        observer=aia_map.observer_coordinate
    )
)

overlay['delta'].add_tickable_gridline('limb', 960*u.arcsec - 90*u.deg)

overlay['psi'].grid(color='red', linestyle='dashed')
overlay['psi'].set_ticks(spacing=30*u.deg)
overlay['psi'].set_ticks_visible(False)
overlay['psi'].set_ticklabel(color='red')
overlay['psi'].set_ticklabel_position(('limb',))
overlay['psi'].set_axislabel('Helioprojective Position Angle', color='red')
overlay['psi'].set_axislabel_position('r')

overlay['delta'].set_ticks_visible(False)
overlay['delta'].set_ticklabel_visible(False)
overlay['delta'].set_axislabel('')

In order to reproject this map to a helioprojective radial frame, we use sunpy.map.make_hpr_header() to create an appropriate FITS WCS header.

hpr_header = make_hpr_header(
    aia_map.observer_coordinate,
    (1000, 360),
    1.8*u.arcsec
)

With the new header, reproject the data into the new coordinate frame. The reproject_to() defaults to using the fast reproject.reproject_interp() algorithm, but a different algorithm can be specified (e.g., reproject.reproject_adaptive()).

outmap = aia_map.reproject_to(hpr_header)

Now we plot the reprojected map. First we add an overlay for the original helioprojective Cartesian frame, but disable its associated ticks and labels.

fig = plt.figure()
ax = fig.add_subplot(projection=outmap)
outmap.plot(axes=ax)

ax.coords[0].set_ticks(spacing=30*u.deg)
ax.set_aspect(0.3)

overlay = ax.get_coords_overlay(aia_map.coordinate_frame)
overlay.grid(color='blue', linestyle='dashed')
for coord in overlay:
    coord.set_ticks(spacing=500*u.arcsec)
    coord.set_ticks_visible(False)
    coord.set_ticklabel_visible(False)
    coord.set_axislabel('')

Due to the FITS WCS machinery, the impact-angle component of helioprojective radial is instead stored as a declination by subtracting 90 degrees, and the vertical axis being declination can be unfamiliar. To change the declination tick labels to impact-angle tick labels, we use show_hpr_impact_angle() to modify the rendering of the tick labels on the appropriate axis (ax.coords[1]).

show_hpr_impact_angle(ax.coords[1])
ax.coords[1].set_format_unit('arcsec')

plt.show()

🏷 Tags: Map, Reproject, AIA, Coordinates