

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

Read in the data



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)

Generate cutouts



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]:





<matplotlib.image.AxesImage at 0x6f70ceb10>



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]:





<matplotlib.image.AxesImage at 0x6f78cbc90>

g-band



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]:





<matplotlib.image.AxesImage at 0x6f5ff7890>