notebook.community

Edit and run 





In [2]:



    


%pylab inline
# Matplotlib default settings
rcdef = plt.rcParams.copy()
pylab.rcParams['figure.figsize'] = 12, 10
pylab.rcParams['xtick.major.size'] = 8.0
pylab.rcParams['xtick.major.width'] = 1.5
pylab.rcParams['xtick.minor.size'] = 4.0
pylab.rcParams['xtick.minor.width'] = 1.5
pylab.rcParams['ytick.major.size'] = 8.0
pylab.rcParams['ytick.major.width'] = 1.5
pylab.rcParams['ytick.minor.size'] = 4.0
pylab.rcParams['ytick.minor.width'] = 1.5
rc('axes', linewidth=2)

import numpy as np
from astropy.io import fits 
from __future__ import division 
from astropy import units as u

import cubehelix  # Cubehelix color scheme from https://github.com/jradavenport/cubehelix
import copy



    


















    






Populating the interactive namespace from numpy and matplotlib

Parameters that are need to know:

  1. Make sure that correct WCS information is available in the header, or provide the pixel scale in both X and Y direction
  2. The parameters for background estimation in the Cold/Hot SEP run:
    • bSizeCold/fSizeCold: Sky background size;
    • bSizeHot / fSizeHot : Sky background filter size
  3. The parameters for object detections in the Cold/Hot SEP run (use Cold run as example, replace Cold with Hot for hot run):
    • nSigCold: thrCold = (nSigCold * bkgC.globalrms);
    • minDetPixCold: minimum number of pixels for a detection
    • convKerCold: 2-D array for the convolution kernel
    • deblendThrCold: deblend threshold TODO: Find ways to choose default values of these parameters.
  4. Location and Shape parameters for the galaxy (or galaxies) you want to fit:
    • galX, galY: The X/Y coordinates in pixel unit
    • galR: Size of the region dominated by the light from the galaxies
    • galQ, galPA: Axis ratio and size of the galaxy
  5. Increase the size of the mask by how many factors: mskGrowHot / mskGrowCold

Optional parameters

  1. contrast for determing the ZSCALE;
  2. options for the Cubehelix: start, rot, gamma
  3. skySigClipping: n-sigma clipping for mean sky value
  4. coaddExptime / coaddNcombine: Total exposure time and NCombine for "expTime corrected image" option

Read-in the Data





In [3]:



    


# Example of cutout image, z~0.4 BCG
#file = 'red_21572_9559_3,7_HSC-I.fits'  
file = 'red_21572_Icut.fits'  
data = fits.open(file)



    











In [4]:



    


#psfFile = 'hsc_coadd-9559-3,7-HSC-I_psf.fits'
psfFile = 'red_21572_Ipsf.fits'
psfData = fits.open(psfFile)
  • Right now, the cutout images from HSC is multi-HDU dataset
  • NOTICE: The image is not stored in the primary HDU
  • The 1/2/3 HDUs are Image/Mask/Variance




In [5]:



    


# Information and structure of the cutout image
data.info()



    


















    






Filename: red_21572_Icut.fits
No.    Name         Type      Cards   Dimensions   Format
0    PRIMARY     PrimaryHDU      20   ()              
1    IMAGE       ImageHDU        34   (801, 801)   float32   
2    MASK        ImageHDU        47   (801, 801)   int16   
3    VARIANCE    ImageHDU        34   (801, 801)   float32   
4    ARCHIVE_INDEX  BinTableHDU     40   305R x 7C    [1J, 1J, 1J, 1J, 1J, 64A, 64A]   
5    ARCHIVE_DATA  BinTableHDU     54   3R x 10C     [1J, 1J, 1J, 1J, 1J, 1K, 2J, 2J, 1J, 1D]   
6    ARCHIVE_DATA  BinTableHDU     68   21R x 9C     [2D, 2D, 4D, 72A, 72A, 1D, 72A, 72A, 72A]   
7    ARCHIVE_DATA  BinTableHDU     62   15R x 12C    [1J, 1J, 1J, 1J, 1J, 1K, 2J, 2J, 1J, 1K, 1J, 1D]   
8    ARCHIVE_DATA  BinTableHDU     59   15R x 8C     [1J, 1J, 1J, 1J, 1J, 1J, 2D, 1E]   
9    ARCHIVE_DATA  BinTableHDU     58   15R x 8C     [2J, 1J, 6D, 6D, 3J, 39366E, 2D, 2D]   
10   ARCHIVE_DATA  BinTableHDU     42   15R x 4C     [100D, 100D, 100D, 100D]   
11   ARCHIVE_DATA  BinTableHDU     35   192R x 2C    [64A, 1J]   
12   ARCHIVE_DATA  BinTableHDU    117   90R x 4C     [1J, 2J, 2J, 9D]   
13   ARCHIVE_DATA  BinTableHDU    117   90R x 4C     [1J, 2J, 2J, 1D]   
14   ARCHIVE_DATA  BinTableHDU     30   81R x 1C     [2D]   
15   ARCHIVE_DATA  BinTableHDU     44   12R x 5C     [1X, 2J, 2J, 1J, 1D]   
16   ARCHIVE_DATA  BinTableHDU     39   180R x 4C    [1J, 1J, 1J, 1D]   
17   ARCHIVE_DATA  BinTableHDU     28   1R x 4C      [1J, 1J, 2D, 32A]   
18   ARCHIVE_DATA  BinTableHDU     50   15R x 9C     [1J, 1J, 1J, 1J, 1J, 1K, 2J, 2J, 1D]   

















In [6]:



    


img = data[1].data  # Image array 
msk = data[2].data  # Mask plane
var = data[3].data  # Variance array 

sig = np.sqrt(var)  # Change the variance array into sigma image



    











In [7]:



    


# Get the headers
imgHead = data[1].header # Correct WCS information should be included
mskHead = data[2].header
varHead = data[3].header



    











In [8]:



    


imgHead # Standard header information



    


















    
Out[8]:








XTENSION= 'IMAGE   '           / IMAGE extension                                
BITPIX  =                  -32 / number of bits per data pixel                  
NAXIS   =                    2 / number of data axes                            
NAXIS1  =                  801 / length of data axis 1                          
NAXIS2  =                  801 / length of data axis 2                          
PCOUNT  =                    0 / required keyword; must = 0                     
GCOUNT  =                    1 / required keyword; must = 1                     
EQUINOX =                2000. / Equinox of coordinates                         
RADESYS = 'ICRS    '           / Coordinate system for equinox                  
CRPIX1  =                3408. / WCS Coordinate reference pixel                 
CRPIX2  =              -12457. / WCS Coordinate reference pixel                 
CD1_1   = -4.66666666666667E-05 / WCS Coordinate scale matrix                   
CD1_2   =                   0. / WCS Coordinate scale matrix                    
CD2_1   =                   0. / WCS Coordinate scale matrix                    
CD2_2   = 4.66666666666667E-05 / WCS Coordinate scale matrix                    
CRVAL1  =     133.333333333333 / WCS Ref value (RA in decimal degrees)          
CRVAL2  =     0.74380165289257 / WCS Ref value (DEC in decimal degrees)         
CUNIT1  = 'deg     '                                                            
CUNIT2  = 'deg     '                                                            
CTYPE1  = 'RA---TAN'           / WCS Coordinate type                            
CTYPE2  = 'DEC--TAN'           / WCS Coordinate type                            
LTV1    =               -14592                                                  
LTV2    =               -30457                                                  
INHERIT =                    T                                                  
EXTTYPE = 'IMAGE   '                                                            
EXTNAME = 'IMAGE   '                                                            
CRVAL1A =                14592 / Column pixel of Reference Pixel                
CRVAL2A =                30457 / Row pixel of Reference Pixel                   
CRPIX1A =                    1 / Column Pixel Coordinate of Reference           
CRPIX2A =                    1 / Row Pixel Coordinate of Reference              
CTYPE1A = 'LINEAR  '           / Type of projection                             
CTYPE2A = 'LINEAR  '           / Type of projection                             
CUNIT1A = 'PIXEL   '           / Column unit                                    
CUNIT2A = 'PIXEL   '           / Row unit                                       

















In [63]:



    


# Get the pixel scale of the image 
try: 
    pixScaleX = (np.abs(imgHead['CD1_1']) * 3600.0)
    pixScaleY = (np.abs(imgHead['CD2_2']) * 3600.0)
except KeyError:
    pixScaleX = pixScaleY = 0.168
else: 
    pixScaleX = pixScaleY = 0.168

print "The pixel scale in X/Y directions are %7.4f / %7.4f arcsecs" % (pixScaleX, pixScaleY)

# Get the image size 
imgSizeX = imgHead['NAXIS1']
imgSizeY = imgHead['NAXIS2']

print "The image size in X/Y directions are %d / %d pixels" % (imgSizeX, imgSizeY)
print "                            %10.2f / %10.2f arcsecs" % (imgSizeX * pixScaleX, imgSizeY * pixScaleY)



    


















    






The pixel scale in X/Y directions are  0.1680 /  0.1680 arcsecs
The image size in X/Y directions are 801 / 801 pixels
                                134.57 /     134.57 arcsecs

Read in the PSF Image

  • TODO: Right now the PSF image is not extracted using the exact coordinate of the galaxy




In [10]:



    


psfData.info()
imgPSF = psfData[0].data



    


















    






Filename: red_21572_Ipsf.fits
No.    Name         Type      Cards   Dimensions   Format
0    PRIMARY     PrimaryHDU      16   (41, 41)     float64   

















In [11]:



    


def zscale(img, contrast=0.25, samples=500):

    # Image scaling function form http://hsca.ipmu.jp/hscsphinx/scripts/psfMosaic.html
    ravel = img.ravel()
    if len(ravel) > samples:
        imsort = np.sort(np.random.choice(ravel, size=samples))
    else:
        imsort = np.sort(ravel)

    n = len(imsort)
    idx = np.arange(n)

    med = imsort[n/2]
    w = 0.25
    i_lo, i_hi = int((0.5-w)*n), int((0.5+w)*n)
    p = np.polyfit(idx[i_lo:i_hi], imsort[i_lo:i_hi], 1)
    slope, intercept = p

    z1 = med - (slope/contrast)*(n/2-n*w)
    z2 = med + (slope/contrast)*(n/2-n*w)

    return z1, z2

Show the Input Image





In [12]:



    


fig = plt.figure(figsize=(10, 10))
fig.subplots_adjust(hspace=0.1, wspace=0.1,
                    top=0.95, right=0.95)
ax = gca()
fontsize = 14
ax.minorticks_on()

for tick in ax.xaxis.get_major_ticks():
    tick.label1.set_fontsize(fontsize)
for tick in ax.yaxis.get_major_ticks():
    tick.label1.set_fontsize(fontsize)
ax.set_title('i-band Image', fontsize=25, fontweight='bold')
ax.title.set_position((0.5,1.01))

imin, imax = zscale(img, contrast=0.05, samples=500)

ax.imshow(np.arcsinh(img), interpolation="none", 
           vmin=imin, vmax=imax,
           cmap=cubehelix.cmap(start=0.5, rot=-0.8, minSat=1.2, maxSat=1.2, 
                               minLight=0., maxLight=1., gamma=1.0))



    


















    
Out[12]:








<matplotlib.image.AxesImage at 0x10479be50>