SunPy test of the astropy reproject. This requires this branch of SunPy: https://github.com/VaticanCameos/sunpy/tree/coordinates-framework and this branch of astropy: https://github.com/astrofrog/astropy/tree/image-reproject

reproject.py

# -*- coding: utf-8 -*-
"""
Created on Sun Jul 13 19:35:54 2014

@author: stuart
"""
import copy

import numpy as np
import matplotlib.pyplot as plt

from astropy.wcs import WCS, WCSSUB_CELESTIAL
from astropy.wcs.utils import wcs_to_celestial_frame, WCS_FRAME_MAPPINGS

from astropy.image.reproject import reproject_image_2d

import sunpy.map
import sunpy.coordinates


def solar_wcs_to_frame(wcs):
    wcs = wcs.sub([WCSSUB_CELESTIAL])

    xcoord = wcs.wcs.ctype[0][:4]
    ycoord = wcs.wcs.ctype[1][:4]
    
    if xcoord == b'HPLN' and ycoord == b'HPLT':
        frame = sunpy.coordinates.HelioProjective()
    
    else:
        frame = None
    
    return frame

WCS_FRAME_MAPPINGS.append(solar_wcs_to_frame)

def get_wcs(amap):
    w2 = WCS(naxis=2)
    w2.wcs.crpix = [amap.reference_pixel['x'], amap.reference_pixel['y']]
    w2.wcs.cdelt = [amap.scale['x'], amap.scale['y']]
    w2.wcs.crval = [amap.reference_coordinate['x'], amap.reference_coordinate['y']]
    w2.wcs.ctype = [amap.coordinate_system['x'], amap.coordinate_system['y']]
    w2.wcs.pc = amap.rotation_matrix
    w2.wcs.cunit = [amap.units['x'], amap.units['y']]   

    return w2

out_ctype = ['HPLN-CAR', 'HPLT-CAR']

aiamap = sunpy.map.Map("~/SyncBox/Programming/SunPy/aia.lev1.193A_2012-06-10T13:00:07.84Z.image_lev1.fits")

in_wcs = get_wcs(aiamap)
out_wcs = copy.deepcopy(in_wcs)
out_wcs.wcs.ctype = out_ctype
out_wcs.wcs.pc = np.identity(2)
out_wcs.wcs.crpix = [512, 1024]

arr_new = reproject_image_2d(aiamap.data, in_wcs, out_wcs,
                             (4096,4096), mode='interpolation', order=0)

plt.imshow(arr_new)
plt.show()