build_wisemosaic.py

from astropy.io import fits
from astroquery.irsa import Irsa
from astropy.wcs import wcs
from astropy.coordinates import SkyCoord
import astropy.units as u
import urllib
import os
import montage_wrapper as montage
import numpy as np

filename = 'wise22_cygnusx_mosaic.fits'
outfile = 'allwise_cygnusx.fits'
w = wcs.WCS(filename)
corners = w.calc_footprint(fits.getheader(filename))
cosdec = np.cos(w.wcs.crval[1]*np.pi/180)
#Create a buffer around the image containing the target area.
offsets = np.array([[+1/cosdec,-1],
                    [+1/cosdec,+1],
                    [-1/cosdec,+1],
                    [-1/cosdec,-1]])*2.0
corners = corners+offsets
corners = [SkyCoord(ra=cc[0],dec=cc[1],unit='deg',frame='fk5') for cc in corners]
table = Irsa.query_region(catalog='wise_allwise_p3am_cdd',spatial='Polygon',polygon = corners,width=2*u.deg)

b4 = table['band'==4]

os.mkdir('tmp')

for row in table:
    if row['band']==4:
        coadd = row['coadd_id']    
        uri = 'http://irsa.ipac.caltech.edu/ibe/data/wise/allwise/p3am_cdd/'
        DownloadName = row['coadd_id']+'-w4-int-3.fits.gz'
        uri = uri+'/'+coadd[0:2]+'/'+coadd[0:4]+'/'+coadd+'/'+DownloadName
        print(uri)
        TheFile = urllib.URLopener()
        TheFile.retrieve(uri,'tmp/'+DownloadName)

os.remove('tmphdr.txt')
header = fits.getheader(filename)
header.totextfile('tmphdr.txt')
montage.mosaic('tmp','wise_mosaic',
               header = 'tmphdr.txt',
               background_match=True)
hdu = fits.open('wise_mosaic/mosaic.fits')
data = hdu[0].data
hdr = hdu[0].header
# Use original pixel area.
PixelAreaArcsec2 = np.abs(np.linalg.det(
    wcs.WCS('tmp/'+DownloadName).pixel_scale_matrix*3600))
PixelAreaSr = PixelAreaArcsec2/4.254e10
data = data * 5.326e-11 / PixelAreaSr
hdr['BUNIT'] = 'MJy/sr'
hduout = fits.PrimaryHDU(data,header=hdr)
hduout.writeto(outfile,clobber=True)