Made with WCSAxes

sun_zoom_plot.py

import numpy as np
from astropy.wcs import WCS
from astropy.io import fits
from astropy import units as u
import matplotlib.pyplot as plt
from wcsaxes import WCSAxes
from sunpy.cm import get_cmap

## Definition of solar transforms (needed for now) ##

from sunpy.wcs import (convert_hcc_hpc, convert_hg_hcc,
                       convert_hcc_hg, convert_hpc_hcc)

from wcsaxes.transforms import CurvedTransform

class SolarHGtoHPC(CurvedTransform):

    def transform(self, input_coords):
        ilon, ilat = input_coords[:,0], input_coords[:,1]
        olon, olat = convert_hcc_hpc(*convert_hg_hcc(ilon, ilat))
        olon /= 3600.
        olat /= 3600
        return np.vstack([olon, olat]).transpose()

    transform_non_affine = transform

    def inverted(self):
        return SolarHPCtoHG()

class SolarHPCtoHG(CurvedTransform):

    def transform(self, input_coords):
        ilon, ilat = input_coords[:,0], input_coords[:,1]
        olon, olat = convert_hcc_hg(*convert_hpc_hcc(ilon * 3600., ilat * 3600., z=True))
        return np.vstack([olon, olat]).transpose()

    transform_non_affine = transform

    def inverted(self):
        return SolarHGtoHPC()

# Why not?
plt.rcParams['font.family'] = 'Arial'

# Read in solar data
hdu = fits.open('aia.lev1.171A_2012-01-16T17_30_00.34Z.image_lev1.fits')[0]

# Make colorscale figure
fig = plt.figure(figsize=(6,6))

# Create axes
ax = WCSAxes(fig, [0.1, 0.1, 0.8, 0.8], wcs=WCS(hdu.header))
fig.add_axes(ax)

# Show colormap
ax.imshow(np.arcsinh(hdu.data / 50.),
          origin='lower', cmap=get_cmap('sdoaia171'))

# Esthetics

x = ax.coords['hpln']
y = ax.coords['hplt']

x.set_ticks(color='white')
y.set_ticks(color='white')

x.set_ticklabel(size='x-small')
y.set_ticklabel(size='x-small')

x.set_axislabel('Solar-x', weight='bold')
y.set_axislabel('Solar-y', weight='bold')

x.set_ticks_position('bl')
y.set_ticks_position('bl')

ax.coords.grid(color='white', alpha=0.1)

# Show overlay with solar lon/lat

overlay = ax.get_coords_overlay(SolarHGtoHPC())

lon = overlay[0]
lat = overlay[1]

lon.set_major_formatter('dd')

lon.set_ticklabel(size='x-small')
lat.set_ticklabel(size='x-small')

lon.set_axislabel('Solar Longitude', weight='bold')
lat.set_axislabel('Solar Latitude', weight='bold')

lon.set_ticks_position('tr')
lat.set_ticks_position('tr')


lon.set_ticks(spacing=5. * u.deg, color='white')
lat.set_ticks(spacing=5. * u.deg, color='white')

overlay.grid(color='white', alpha=0.3)

# Set title

ax.set_title("Some nearby star", y=1.12, size='large', weight='bold')

# Zoom in

ax.set_xlim(300., 1000.)
ax.set_ylim(2200., 2900.)

fig.savefig('sun_zoom.png', dpi=150, bbox_inches='tight')