Reproject JP2 files

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 reproject import reproject_interp


def limb_coords(smap):
    r = smap.rsun_obs - 1 * u.arcsec  # remove one arcsec so it's on disk.
    # Adjust the following range if you only want to plot on STEREO_A
    th = np.linspace(-180 * u.deg, 0 * u.deg)
    x = r * np.sin(th)
    y = r * np.cos(th)
    return SkyCoord(x, y, frame=smap.coordinate_frame)



#map_aia = sunpy.map.Map(AIA_171_ROLL_IMAGE)
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(aia_file)
first_wcs =  first_map.wcs.deepcopy()
first_limb = limb_coords(first_map)

# Load in the second map
second_map = sunpy.map.Map(stereo_file)
second_wcs = second_map.wcs.deepcopy()
second_limb = limb_coords(second_map)

# We want to change the pixel scale and the image size, wcs expects cdelt in deg.
#venus_wcs = map_aia.wcs.deepcopy()
#venus_wcs.wcs.cdelt = ([1.5, 1.5]*u.arcsec).to(u.deg)
#venus_wcs.wcs.crpix = [1024, 1024]
#venus = get_body_heliographic_stonyhurst('venus', map_aia.date)
#venus_wcs.heliographic_observer = venus

# Reproject the 
output, footprint = reproject_interp((first_map.data, first_map.wcs), second_wcs, (2048, 2048), order='nearest-neighbor')

out_header = copy.deepcopy(first_map.meta)
out_header.update(second_wcs.to_header())

outmap = sunpy.map.Map((output, out_header))
outmap._remove_existing_observer_location()
outmap.meta['hgln_obs'] = second_map.observer_coordinate.lon.to_value(u.deg)
outmap.meta['hglt_obs'] = second_map.observer_coordinate.lat.to_value(u.deg)
outmap.meta['dsun_obs'] = second_map.observer_coordinate.radius.to_value(u.m)

fig = plt.figure(figsize=(15,5))

ax1 = fig.add_subplot(1,3,1, projection=first_map)
first_map.plot(axes=ax1)
first_map.draw_grid()
ax1.plot_coord(second_limb, color='r')

ax2 = fig.add_subplot(1,3,2, projection=second_map)
second_map.plot(axes=ax2)
second_map.draw_grid(color='k')
ax2.plot_coord(first_limb, color='r')

ax3 = fig.add_subplot(1,3,3, projection=outmap)
outmap.plot(axes=ax3)
outmap.draw_grid(color='r')
ax3.set_title('First map data visible from second map observer location\n(no solar differential rotation)')
ax3.grid('on', linestyle=":", color='c')
#ax3.plot_coord(first_limb, color='o')