stitch_fields_bigger_DRc.py

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

base_header = fits.getheader('moments/G0_13CO21_downsampled_mom0.fits')
base_wcs = wcs.WCS(base_header).celestial

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['NAXIS'] = 2

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

for fn in ('moments/G001_13CO21_SEDIGISM_DR1c_mom0.fits',
           'moments/G000_13CO21_SEDIGISM_DR1c_mom0.fits',
           'moments/G002_13CO21_SEDIGISM_DR1c_mom0.fits',
           'moments/G003_13CO21_SEDIGISM_DR1c_mom0.fits',
           'moments/G004_13CO21_SEDIGISM_DR1c_mom0.fits',
           'moments/G005_13CO21_SEDIGISM_DR1c_mom0.fits',
           'moments/G359_13CO21_SEDIGISM_DR1c_mom0.fits',
           'moments/G358_13CO21_SEDIGISM_DR1c_mom0.fits',
           'moments/G356_13CO21_SEDIGISM_DR1c_mom0.fits',
           'moments/G355_13CO21_SEDIGISM_DR1c_mom0.fits',
           'moments/G357_13CO21_SEDIGISM_DR1c_mom0.fits',):
    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('moments/G0_larger_13CO21_mom0_mosaic_DR1c.fits', overwrite=True)


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

for fn in ('max/G001_13CO21_SEDIGISM_DR1c_max_K.fits',
           'max/G000_13CO21_SEDIGISM_DR1c_max_K.fits',
           'max/G002_13CO21_SEDIGISM_DR1c_max_K.fits',
           'max/G003_13CO21_SEDIGISM_DR1c_max_K.fits',
           'max/G004_13CO21_SEDIGISM_DR1c_max_K.fits',
           'max/G005_13CO21_SEDIGISM_DR1c_max_K.fits',
           'max/G359_13CO21_SEDIGISM_DR1c_max_K.fits',
           'max/G358_13CO21_SEDIGISM_DR1c_max_K.fits',
           'max/G356_13CO21_SEDIGISM_DR1c_max_K.fits',
           'max/G355_13CO21_SEDIGISM_DR1c_max_K.fits',
           'max/G357_13CO21_SEDIGISM_DR1c_max_K.fits',):
    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('max/G0_larger_13CO21_max_K_mosaic_DR1c.fits', overwrite=True)


for fn in ('moments/G001_13CO21_SEDIGISM_DR1c_mom1.fits',
           'moments/G000_13CO21_SEDIGISM_DR1c_mom1.fits',
           'moments/G002_13CO21_SEDIGISM_DR1c_mom1.fits',
           'moments/G003_13CO21_SEDIGISM_DR1c_mom1.fits',
           'moments/G004_13CO21_SEDIGISM_DR1c_mom1.fits',
           'moments/G005_13CO21_SEDIGISM_DR1c_mom1.fits',
           'moments/G359_13CO21_SEDIGISM_DR1c_mom1.fits',
           'moments/G358_13CO21_SEDIGISM_DR1c_mom1.fits',
           'moments/G356_13CO21_SEDIGISM_DR1c_mom1.fits',
           'moments/G355_13CO21_SEDIGISM_DR1c_mom1.fits',
           'moments/G357_13CO21_SEDIGISM_DR1c_mom1.fits',):
    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('moments/G0_larger_13CO21_mom1_mosaic_DR1c.fits', overwrite=True)