fusemusealma.py

def fusemusealma(muse_moment_fits, alma_moment_fits):
    '''
    Open MUSE and ALMA HDUs, pick up their WCSs, 
    drop needless ALMA WCS axes, reproject 
    ALMA to MUSE, then output new HDUs.  
    '''
    museHDU = fits.open(get_pkg_data_filename(muse_moment_fits))
    almaHDU = fits.open(get_pkg_data_filename(alma_moment_fits))

    w_alma = WCS(almaHDU[0])
    w_muse = WCS(museHDU[0])

    alma_3wcs = w_alma.dropaxis(3)
    alma_2wcs = alma_3wcs.dropaxis(2)

    newalmaHeader = alma_2wcs.to_header()
    newalmaHDU = fits.PrimaryHDU(data=almaHDU[0].data, header=newalmaHeader)
    
    registered_alma_data, registered_alma_footprint = reproject_interp(newalmaHDU, museHDU[0].header)
    registered_alma_data[registered_alma_data==0] = np.nan
    
    ax1 = plt.subplot(1,2,1, projection=WCS(museHDU[0].header))
    ax1.imshow(museHDU[0].data, origin='lower')
    ax1.coords['ra'].set_axislabel('Right Ascension')
    ax1.coords['dec'].set_axislabel('Declination')
    ax1.set_title('MUSE Original')

    ax2 = plt.subplot(1,2,2, projection=WCS(newalmaHDU.header))
    ax2.imshow(registered_alma_data, origin='lower')
    ax2.coords['ra'].set_axislabel('Right Ascension')
    ax2.coords['dec'].set_axislabel('Declination')
    ax2.set_title('Registered ALMA Image')
    
    ax1.set_xlim(100, 220)
    ax1.set_ylim(100,300)
    ax2.set_xlim(100, 220)
    ax2.set_ylim(100,300)
    
    plt.show()
    print("Successfully reprojected {} to {}".format(alma_moment_fits, muse_moment_fits))
    
    return museHDU, almaHDU