Example script for creating composite RGB images from MGN-filtered AIA observations

mgn-rgb-example.py

"""
Create a composite RGB image using MGN-processed images from SDO/AIA 211, 193, and 171
"""
import astropy.units as u
from astropy.coordinates import SkyCoord
from astropy.visualization import make_lupton_rgb,ImageNormalize
import numpy as np
import matplotlib.pyplot as plt
from sunpy.net import Fido, attrs as a
from sunpy.map import Map
from aiapy.calibrate import register, update_pointing, normalize_exposure
from sunkit_image.enhance import mgn


q = Fido.search(
    a.Instrument.aia,
    a.Time('2012-09-24T13:00:00', '2012-09-24T13:00:15'),
    a.Wavelength(171*u.angstrom,211*u.angstrom),
    a.Physobs.intensity,
)
f = Fido.fetch(q)

m_l1 = Map(sorted(f))
m_l15 = [normalize_exposure(register(update_pointing(m))) for m in m_l1]

blc = SkyCoord(Tx=-350*u.arcsec, Ty=-200*u.arcsec, frame=m_l15[0].coordinate_frame)
trc = SkyCoord(Tx=155*u.arcsec, Ty=250*u.arcsec, frame=m_l15[0].coordinate_frame)
m_l15_cutout = [m.submap(blc, top_right=trc) for m in m_l15]

m_l15_cutout_mgn = [Map(mgn(m.data, k=5), m.meta) for m in m_l15_cutout]

im_mgn_rgb = make_lupton_rgb(
    m_l15_cutout_mgn[2].data,
    m_l15_cutout_mgn[1].data,
    m_l15_cutout_mgn[0].data,
    Q=0,
    stretch=1,
)

fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111, projection=m_l15_cutout_mgn[0].wcs)
im = ax.imshow(im_mgn_rgb)
lon,lat = ax.coords
lon.set_axislabel('Helioprojective Longitude')
lat.set_axislabel('Helioprojective Latitude')
fig.savefig('mgn-rgb-211-193-171.png', dpi=200, bbox_inches='tight')