In [1]:
%matplotlib inline
from __future__ import (division,
print_function)
import os
import numpy as np
from astropy.io import fits
from astropy import wcs
from astropy import units as u
from astropy.nddata import Cutout2D
from astropy.coordinates import SkyCoord
from astropy import coordinates as coords
import matplotlib.pyplot as plt
In [2]:
baseDir = "/Volumes/astro5/project/perseus/work/"
dataDir = os.path.join(baseDir, 'data')
gImg = os.path.join(dataDir, 'PERSEUS30sec_HSC-G_img.fits')
iImg = os.path.join(dataDir, 'PERSEUS30sec_HSC-I_img.fits')
gVar = os.path.join(dataDir, 'PERSEUS30sec_HSC-G_variance.fits')
iVar = os.path.join(dataDir, 'PERSEUS30sec_HSC-I_variance.fits')
In [3]:
# HSC-I band images
iHeader = fits.open(iImg)[0].header
iData = fits.open(iImg)[0].data
iVar = fits.open(iVar)[0].data
# HSC-G band images
gHeader = fits.open(gImg)[0].header
gData = fits.open(gImg)[0].data
gVar = fits.open(gVar)[0].data
# Read WCS information from the header
iWcs = wcs.WCS(iHeader)
gWcs = wcs.WCS(gHeader)
In [4]:
def img_cutout(img, wcs, ra, dec,
imgsize=4.0, pix=0.168,
prefix='perseus_cutout'):
"""
Generate cutout image.
"""
imgsize = (4.0 * 60.0) / pix
cenX, cenY = wcs.wcs_world2pix(ra, dec, 0)
cenPos = (int(cenX), int(cenY))
shape = (imgsize, imgsize)
# Generate cutout
cutout = Cutout2D(img, cenPos, shape, wcs=wcs)
# Update the header
header = cutout.wcs.to_header()
# Build a HDU
hdu = fits.PrimaryHDU(header=header)
hdu.data = cutout.data
# Save FITS image
fitsFile = prefix + '.fits'
hdu.writeto(fitsFile, overwrite=True)
return cutout
In [6]:
udg08_radec = ["3:18:07.782 +41:27:05.106"]
udg10_radec = ["3:18:56.291 +41:15:01.878"]
# R16 3:18:36.384 +41:11:31.075
# R24 3:18:35.750 +41:48:30.462
udg08_coord = SkyCoord(udg08_radec, unit=(u.hourangle, u.deg))
udg10_coord = SkyCoord(udg10_radec, unit=(u.hourangle, u.deg))
In [7]:
udg08_img = img_cutout(iData, iWcs,
udg08_coord.ra.value[0],
udg08_coord.dec.value[0],
imgsize=4.0,
prefix='perseus_udg_08_i_img')
udg08_var = img_cutout(iVar, iWcs,
udg08_coord.ra.value[0],
udg08_coord.dec.value[0],
imgsize=4.0,
prefix='perseus_udg_08_i_var')
In [33]:
plt.imshow(np.arcsinh(udg08_img.data), origin='lower')
Out[33]:
In [8]:
udg10_img = img_cutout(iData, iWcs,
udg10_coord.ra.value[0],
udg10_coord.dec.value[0],
imgsize=4.0,
prefix='perseus_udg_10_i_img')
udg10_var = img_cutout(iVar, iWcs,
udg10_coord.ra.value[0],
udg10_coord.dec.value[0],
imgsize=4.0,
prefix='perseus_udg_10_i_var')
In [36]:
plt.imshow(np.arcsinh(udg10_img.data), origin='lower')
Out[36]:
In [9]:
udg08_img = img_cutout(gData, gWcs,
udg08_coord.ra.value[0],
udg08_coord.dec.value[0],
imgsize=4.0,
prefix='perseus_udg_08_g_img')
udg08_var = img_cutout(gVar, gWcs,
udg08_coord.ra.value[0],
udg08_coord.dec.value[0],
imgsize=4.0,
prefix='perseus_udg_08_g_var')
In [10]:
plt.imshow(np.arcsinh(udg10_img.data), origin='lower')
Out[10]: