hpr_tests.py

import astropy.units as u
from astropy.coordinates import SkyCoord, SkyOffsetFrame
import sunpy.coordinates
import sunpy.map
from sunpy.data.sample import AIA_171_IMAGE
import matplotlib.pyplot as plt
import numpy as np
plt.ion()

aia = sunpy.map.Map(AIA_171_IMAGE)

origin = SkyCoord(0, 90, unit=u.deg, frame=aia.coordinate_frame)
sof = SkyOffsetFrame(origin=origin)

intensities = []
coords = []
ring_lats = -90*u.deg + np.linspace(1000, 1200, 5, endpoint=True) * u.arcsec
for lat in ring_lats:
    coord = SkyCoord(np.linspace(0, 360, 360), lat, unit=u.deg, frame=sof)
    intensity, coord = sunpy.map.extract_along_coord(aia, coord)
    intensities.append(intensity)
    coords.append(coord)


prop_cycle = plt.rcParams['axes.prop_cycle']
colors = prop_cycle.by_key()['color']

fig = plt.figure()
ax = fig.add_subplot(1, 2, 1, projection=aia)
im = aia.plot(axes=ax, clip_interval=(1, 99)*u.percent, annotate=False)
tx, ty = ax.coords
tx.set_ticks_visible(False)
ty.set_ticks_visible(False)
tx.set_ticklabel_visible(False)
ty.set_ticklabel_visible(False)

for i, coord in enumerate(coords):
    ax.plot_coord(coord, color=colors[i%len(colors)], alpha=0.8, linestyle="--")
overlay = ax.get_coords_overlay(sof)

lon, lat = overlay
lon_color = colors[5]
lon.grid(color=lon_color)
lon.set_ticks(spacing=15*u.deg)
lon.set_ticks_position('all')
lon.set_ticklabel_position('all')
lon.set_axislabel_position('b')
lon.set_ticklabel(color=lon_color)
lon.set_axislabel("Helioprojective Radial Longitude", color=lon_color)

lat_color = "blue"
# lat.grid(color=lat_color)
lat.set_ticks_position('')
lat.set_ticklabel_position('')
lat.set_axislabel_position('')
lat.set_ticklabel(color=lat_color)
lat.set_axislabel("Helioprojective Radial Latitude", color=lat_color)

ax1 = fig.add_subplot(1, 2, 2, projection="polar")

for i, (intensity, coord) in enumerate(zip(intensities, coords)):
    off_coord = coord.transform_to(sof)
    lat_offset = (90*u.deg + off_coord[0].lat).to_string(unit=u.arcsec, decimal=True, format="latex")
    ax1.plot(off_coord.lon.to(u.rad), intensity, color=colors[i%len(colors)], label=f'{lat_offset}\"', linewidth=0.75)
ax1.set_xlabel("Helioprojective Radial Longitude [deg]")
ax1.set_ylabel(f"Intensity [{aia.unit}]", labelpad=30)
ax1.set_theta_zero_location("N")
ax1.set_theta_direction(-1)
ax1.legend()