stitch_cubes.py

import numpy as np
from astropy.io import fits
import reproject
from astropy import wcs

base_header = fits.getheader('G000_13CO21_SEDIGISM_DR1c_spatialdownsample_4x_downsampled_2kms.fits')
base_wcs = wcs.WCS(base_header)

new_header = base_wcs.to_header()

# 6 degrees across
new_header['NAXIS1'] = int(np.ceil(12.0 / np.abs(base_wcs.wcs.cdelt[0])))
new_header['CRPIX1'] = new_header['NAXIS1'] / 2
new_header['NAXIS2'] = base_header['NAXIS2']
new_header['NAXIS3'] = base_header['NAXIS3']
new_header['NAXIS'] = 3

out_data = np.zeros([new_header['NAXIS3'], new_header['NAXIS2'], new_header['NAXIS1'],])

for fn in (
    "G000_13CO21_SEDIGISM_DR1c_spatialdownsample_4x_downsampled_2kms.fits",
    "G001_13CO21_SEDIGISM_DR1c_spatialdownsample_4x_downsampled_2kms.fits",
    "G002_13CO21_SEDIGISM_DR1c_spatialdownsample_4x_downsampled_2kms.fits",
    "G003_13CO21_SEDIGISM_DR1c_spatialdownsample_4x_downsampled_2kms.fits",
    "G004_13CO21_SEDIGISM_DR1c_spatialdownsample_4x_downsampled_2kms.fits",
    "G005_13CO21_SEDIGISM_DR1c_spatialdownsample_4x_downsampled_2kms.fits",
    "G355_13CO21_SEDIGISM_DR1c_spatialdownsample_4x_downsampled_2kms.fits",
    "G356_13CO21_SEDIGISM_DR1c_spatialdownsample_4x_downsampled_2kms.fits",
    "G357_13CO21_SEDIGISM_DR1c_spatialdownsample_4x_downsampled_2kms.fits",
    "G358_13CO21_SEDIGISM_DR1c_spatialdownsample_4x_downsampled_2kms.fits",
    "G359_13CO21_SEDIGISM_DR1c_spatialdownsample_4x_downsampled_2kms.fits",):
    print("Resampling {0}".format(fn))
    hdu = fits.open(fn)
    new_array, weights = reproject.reproject_interp(hdu, new_header)

    out_data[weights.astype('bool')] += new_array[weights.astype('bool')]

outhdu = fits.PrimaryHDU(data=out_data, header=new_header)
outhdu.writeto('GalacticCenterMosiac_downsampled4x_2kms.fits', overwrite=True)