multiple_wcs.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.wcs import *

from wcsaxes.transforms import CurvedTransform

plt.rcParams['font.family'] = 'Arial'

# Read in 2MASS image in Equatorial coordinates
hdu = fits.open('20120222_000921_s4h2B.fts')[0]

fig = plt.figure(figsize=(12, 12))

for wcs_key in [None, 'A']:

    # Make colorscale figure
    if wcs_key is None:
        ax = WCSAxes(fig, [0.1, 0.1, 0.35, 0.35], wcs=WCS(hdu.header))    
    else:
        ax = WCSAxes(fig, [0.53, 0.1, 0.35, 0.35], wcs=WCS(hdu.header, key='A'))    

    fig.add_axes(ax)
    ax.imshow(hdu.data ** 0.5, origin='lower', cmap=plt.cm.gist_heat, vmax=1000)

    lon = ax.coords[0]
    lat = ax.coords[1]

    lon.set_ticks(color='white')
    lat.set_ticks(color='white')

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

    if wcs_key is None:
        lon.set_axislabel('Solar-x', weight='bold')
        lat.set_axislabel('Solar-y', weight='bold')
    else:
        lon.set_axislabel('Right Ascension', weight='bold')
        lat.set_axislabel('Declination', weight='bold')

    lon.set_ticks(spacing=10 * u.deg)
    lat.set_ticks(spacing=10 * u.deg)

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

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