Experiments in helioviewer JP2 image reprojection

gistfile1.txt

import warnings
#warnings.filterwarnings('ignore')
import copy

import numpy as np

import matplotlib.pyplot as plt
import astropy.units as u
from astropy.coordinates import SkyCoord

import sunpy.map
from sunpy.data.sample import AIA_171_ROLL_IMAGE
from sunpy.coordinates.ephemeris import get_body_heliographic_stonyhurst
from sunpy.map.header_helper import make_fitswcs_header
from sunpy.coordinates import frames

from reproject import reproject_exact, reproject_interp

# Locations of some test images
aia_file = '/Users/ireland/Data/jp2/20110820/2011_08_20__16_08_00_34__SDO_AIA_AIA_171.jp2'
stereo_file =  '/Users/ireland/Data/jp2/20110820/2011_08_20__16_14_45_680__STEREO-B_SECCHI_EUVI_171.jp2'

# Load in the first map
first_map = sunpy.map.Map(stereo_file)

# Put in some data to assess how the transformation is working
first_map.data[1024-50:1024+50, 1600-30:1600+30] = 255

# Make the new WCS
# Number of pixels in the x and y direction
nx = int(2048/4)
ny = int(2048/4)

# Calculate the pixel scale so that it spans the physical range of the map
lx = first_map.dimensions.x * first_map.scale.axis1
ly = first_map.dimensions.y * first_map.scale.axis2
cdelt1 = (lx/nx).to(u.arcsec).value
cdelt2 = (ly/ny).to(u.arcsec).value
scale = [cdelt1, cdelt2]*u.arcsec

# Define the observer
obstime = first_map.date
observer = get_body_heliographic_stonyhurst('earth', obstime)
observer = SkyCoord(-45*u.deg, -30*u.deg, 0.5*u.AU, obstime=obstime, frame=frames.HeliographicStonyhurst)
new_observer = SkyCoord(0*u.arcsec, 0*u.arcsec, obstime=obstime, observer=observer, frame=frames.Helioprojective)

# Define the new WCS via a new header
new_data = np.zeros((ny, nx))
new_header = make_fitswcs_header(new_data, new_observer, scale=scale/u.pix, instrument=first_map.instrument,
                                     observatory=first_map.observatory,
                                     wavelength=first_map.wavelength)
new_map = sunpy.map.Map((new_data, new_header))
new_wcs = new_map.wcs.deepcopy()

# Reproject and create the output map
output, footprint = reproject_interp((first_map.data, first_map.wcs), new_wcs, new_data.shape)
outmap = sunpy.map.Map((output, new_header))

# Make the plot
fig = plt.figure(figsize=(10,5))

ax1 = fig.add_subplot(1,2,1, projection=first_map)
first_map.plot(axes=ax1)
first_map.draw_grid()

ax2 = fig.add_subplot(1,2,2, projection=outmap)
outmap.plot(axes=ax2)
outmap.draw_grid(color='r')
loc = 'lon={:n}, lat={:n}, r={:n}AU'.format(observer.lon.to(u.deg).value, observer.lat.to(u.deg).value, observer.radius.to(u.AU).value)
ax2.set_title('First map data visible from new observer location\n{:s}\n(no solar differential rotation)'.format(loc))
ax2.grid('on', linestyle=":", color='c')