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>










    
























In [13]:



    


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 "Background"', fontsize=25, fontweight='bold')
ax.title.set_position((0.5,1.01))

imgMsk = copy.deepcopy(img) 
imgMsk[msk > 0] = np.nan 
    
mmin, mmax = zscale(imgMsk, contrast=0.6, samples=500)
cmap = cubehelix.cmap(start=0.5, rot=-0.8, minSat=1.2, maxSat=1.2, 
                      minLight=0., maxLight=1., gamma=1.0)
cmap.set_bad('k',1.)

ax.imshow(np.arcsinh(imgMsk), interpolation="none", 
          vmin=mmin, vmax=mmax, cmap=cmap)



    


















    
Out[13]:








<matplotlib.image.AxesImage at 0x1080b4bd0>

Show the Mask Plane





In [14]:



    


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 Mask Plane', fontsize=25, fontweight='bold')
ax.title.set_position((0.5,1.01))

ax.yaxis.set_major_formatter(plt.NullFormatter())

ax.imshow(msk, cmap=cubehelix.cmap(reverse=True))



    


















    
Out[14]:








<matplotlib.image.AxesImage at 0x107b523d0>

Show the Sigma Image





In [15]:



    


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 Sigma Image', fontsize=25, fontweight='bold')
ax.title.set_position((0.5,1.01))

ax.yaxis.set_major_formatter(plt.NullFormatter())

smin, smax = zscale(sig, contrast=0.03, samples=500)
ax.imshow(np.arcsinh(sig), interpolation="none", 
          vmin=smin, vmax=smax, 
          cmap=cubehelix.cmap(start=0.5, rot=-0.8, reverse=True, minSat=1.2, maxSat=1.2, 
                               minLight=0., maxLight=1., gamma=0.5))



    


















    
Out[15]:








<matplotlib.image.AxesImage at 0x1029f3690>

Show the PSF images





In [16]:



    


fig, axes = plt.subplots(1, 2, sharex=True, sharey=True, figsize=(10,5))
plt.subplots_adjust(left=0.03, bottom=0.03, right=0.995, top=0.99, 
                    wspace=0.01, hspace=0.01)

# Highlight the Inner part of PSF
pmin, pmax = zscale(imgPSF, contrast=0.05, samples=500)
axes[0].imshow(np.arcsinh(imgPSF), interpolation='none', 
               cmap=cubehelix.cmap(reverse=True, minSat=1.2, maxSat=1.2, 
                                   minLight=0., maxLight=1., gamma=1.0))
axes[0].set_title('Inner PSF', fontsize=26)
axes[0].title.set_position((0.5,1.01))

# Highlight the Outer part of PSF
pmin, pmax = zscale(imgPSF, contrast=0.1, samples=500)
axes[1].imshow(np.arcsinh(imgPSF), interpolation='none', 
               vmin=pmin, vmax=pmax, 
               cmap=cubehelix.cmap(reverse=True, minSat=1.2, maxSat=1.2, 
                                   minLight=0., maxLight=1., gamma=0.5))
axes[1].set_title('Outer PSF', fontsize=26)
axes[1].title.set_position((0.5,1.01))

Show Image Using APLpy Package





In [17]:



    


# Try ApLPY 
import aplpy 
fig = aplpy.FITSFigure(file, hdu=1)
fig.show_grayscale(interpolation='none')
fig.add_label(133.4737, 1.3449933, 'BCG', color='w')
fig.add_label(0.1, 0.9, '(a)', relative=True, color='w', size=25)
fig.show_ellipses(133.4737, 1.3449933, 0.02, 0.03, color='y')

fig.add_grid()

fig.add_scalebar(0.05)
fig.scalebar.set_length(10 * u.arcsecond)
fig.scalebar.set_label('5 asec')
fig.scalebar.set_corner('top right')
fig.scalebar.set_alpha(0.9)
fig.scalebar.set_color('white')



    


















    






WARNING: Cannot determine equinox. Assuming J2000. [aplpy.wcs_util]
WARNING:astropy:Cannot determine equinox. Assuming J2000.
WARNING: Cannot determine equinox. Assuming J2000. [aplpy.wcs_util]
WARNING:astropy:Cannot determine equinox. Assuming J2000.
/Users/songhuang/anaconda/lib/python2.7/site-packages/aplpy/labels.py:432: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal
  if self.coord == x or self.axis.apl_tick_positions_world[ipos] > 0:
INFO:astropy:Auto-setting vmin to -8.522e-01
INFO:astropy:Auto-setting vmax to  5.392e+00










    






INFO: Auto-setting vmin to -8.522e-01 [aplpy.core]
INFO: Auto-setting vmax to  5.392e+00 [aplpy.core]

Histogram for Pixels without Any Mask





In [64]:



    


from astropy.stats import sigma_clip

indNoMask = np.where(msk == 0)
skySigClipping = 3.0

pixNoMask = sigma_clip(img[indNoMask].flatten(), skySigClipping, 3)
pixAll    = img.flatten()



    











In [19]:



    


from astroML.plotting import hist
import scipy 

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

ax.set_xlim(-0.7, 1.0)

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)

counts1, bins2, patches3 = hist(pixNoMask, bins='knuth', ax=ax, alpha=0.7, 
                                color='r', histtype='stepfilled', normed=True)
counts1, bins2, patches3 = hist(pixAll,  bins='knuth', ax=ax, alpha=0.9, 
                                color='k', histtype='step', normed=True, linewidth=2)

ax.axvline(0.0, linestyle='-', color='k', linewidth=1.5)
ax.axvline(np.nanmedian(pixNoMask), linestyle='--', color='b', linewidth=1.5)

ax.set_xlabel('Pixel Value', fontsize=20)

ax.text(0.3, 3.2, "Min   : %8.5f" % np.nanmin(pixNoMask),  fontsize=22)
ax.text(0.3, 2.9, "Max   : %8.5f" % np.nanmax(pixNoMask),  fontsize=22)
ax.text(0.3, 2.6, "Mean  : %8.5f" % np.nanmean(pixNoMask), fontsize=22)
ax.text(0.3, 2.3, "StdDev: %8.5f" % np.nanstd(pixNoMask),  fontsize=22)
ax.text(0.3, 2.0, "Median: %8.5f" % np.nanmedian(pixNoMask),  fontsize=22)
ax.text(0.3, 1.7, "Skew  : %8.5f" % scipy.stats.skew(pixNoMask),  fontsize=22)



    


















    
Out[19]:








<matplotlib.text.Text at 0x112677290>

Making Object Mask for the Central Galaxy

Use the SEP.background to re-estimate the background of the image





In [20]:



    


import sep

bSizeCold = 80
fSizeCold = 20

bSizeHot  =  8
fSizeHot  =  4

imgC = copy.deepcopy(img)
imgC = imgC.byteswap(True).newbyteorder()

imgH = copy.deepcopy(img)
imgH = imgH.byteswap(True).newbyteorder()

bkgC = sep.Background(imgC, None, bSizeCold, bSizeCold, 
                     fSizeCold, fSizeCold)
bkgH = sep.Background(imgH, None, bSizeHot, bSizeHot, 
                     fSizeHot, fSizeHot)

globalBkgC = bkgC.globalback
globalRmsC = bkgC.globalrms

globalBkgH = bkgH.globalback
globalRmsH = bkgH.globalrms

Show the new "Background"





In [21]:



    


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('SEP Background (Cold)', fontsize=25, fontweight='bold')
ax.title.set_position((0.5,1.01))

cmap = cubehelix.cmap(start=0.5, rot=-0.8, minSat=1.2, maxSat=1.2, 
                      minLight=0., maxLight=1., gamma=1.0)
cmap.set_bad('k',1.)

ax.imshow(np.arcsinh(bkgC.back()), interpolation="none", 
          vmin=mmin, vmax=mmax, cmap=cmap)
ax.text(0.6, 0.92, "Global BKG : %7.4f" % globalBkgC, 
        horizontalalignment='center',
        verticalalignment='center',
        transform = ax.transAxes, fontsize=20)
ax.text(0.6, 0.87, "Global RMS : %7.4f" % globalRmsC, 
        horizontalalignment='center',
        verticalalignment='center',
        transform = ax.transAxes, fontsize=20)



    


















    
Out[21]:








<matplotlib.text.Text at 0x10bd97a50>










    
























In [22]:



    


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('SEP Background (Hot)', fontsize=25, fontweight='bold')
ax.title.set_position((0.5,1.01))

cmap = cubehelix.cmap(start=0.5, rot=-0.8, minSat=1.2, maxSat=1.2, 
                      minLight=0., maxLight=1., gamma=1.0)
cmap.set_bad('k',1.)

ax.imshow(np.arcsinh(bkgH.back()), interpolation="none", 
          vmin=mmin, vmax=mmax, cmap=cmap)
ax.text(0.6, 0.92, "Global BKG : %7.4f" % globalBkgH, 
        horizontalalignment='center',
        verticalalignment='center',
        transform = ax.transAxes, fontsize=20)
ax.text(0.6, 0.87, "Global RMS : %7.4f" % globalRmsH, 
        horizontalalignment='center',
        verticalalignment='center',
        transform = ax.transAxes, fontsize=20)



    


















    
Out[22]:








<matplotlib.text.Text at 0x112bd9490>

Subtract that "Background" from the Original Image





In [23]:



    


bkgC.subfrom(imgC) 
bkgH.subfrom(imgH)



    











In [28]:



    


def getEll2Plot(objects, radius=None): 
    
    from matplotlib.patches import Ellipse
    
    x  = objects['x'].copy()
    y  = objects['y'].copy()
    pa = objects['theta'].copy() # in unit of radian

    if radius is not None: 
        a = radius.copy()
        b = radius.copy()*(objects['b'].copy()/objects['a'].copy())
    else: 
        a  = objects['a'].copy()
        b  = objects['b'].copy()
    
    ells = [Ellipse(xy=np.array([x[i], y[i]]), 
                    width=np.array(2.0 * b[i]), 
                    height=np.array(2.0 * a[i]), 
                    angle=np.array(pa[i]*180.0/np.pi + 90.0)) 
            for i in range(x.shape[0])]
    
    return ells

Object Detection

The Low-Threshold Cold Run





In [25]:



    


nSigCold       = 1.5
thrCold        = (nSigCold * bkgC.globalrms)
minDetPixCold  = 5
# Tophat_3.0_3x3
convKerCold    = np.asarray([[0.560000, 0.980000, 0.560000], 
                             [0.980000, 1.000000, 0.980000],
                             [0.560000, 0.980000, 0.560000]])
deblendThrCold = 24
deblendConCold = 0.005

objCold = sep.extract(imgC, thrCold, minarea=minDetPixCold, 
                      conv=convKerCold, deblend_nthresh=deblendThrCold, 
                      deblend_cont=deblendConCold)



    











In [26]:



    


xCold  = objCold['x'].copy()
yCold  = objCold['y'].copy()
aCold  = objCold['a'].copy()
bCold  = objCold['b'].copy()
paCold = objCold['theta'].copy()
fluxCold = objCold['cflux'].copy()

numCold = xCold.shape[0]
print "%d objects are detected in the cold run" % numCold

# TODO: For some reasons, the Kron radius does not work?
radKronCold, flagKronCold = sep.kron_radius(imgC, xCold, yCold, 
                                            aCold, bCold, paCold, 
                                            np.empty(numCold).fill(6), 
                                            mask=None, maskthresh=0.0)



    


















    






276 objects are detected in the cold run

















In [29]:



    


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('Cold-mode Detections', fontsize=25, fontweight='bold')
ax.title.set_position((0.5,1.01))

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

ax.imshow(np.arcsinh(imgC), 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))
ellsCold = getEll2Plot(objCold, radius=(aCold*4.0))

for e in ellsCold:
    ax.add_artist(e)
    e.set_clip_box(ax.bbox)
    e.set_alpha(0.8)
    e.set_edgecolor('b')
    e.set_facecolor('none')
    e.set_linewidth(1.5)

The High-Threshold Hot Run





In [30]:



    


nSigHot       = 3.0
thrHot        = (nSigHot * bkgH.globalrms)
minDetPixHot  = 4
# Tophat_3.0_3x3
convKerHot    = np.asarray([[0.560000, 0.980000, 0.560000], 
                            [0.980000, 1.000000, 0.980000],
                            [0.560000, 0.980000, 0.560000]])
deblendThrHot = 24
deblendConHot = 0.005

objHot = sep.extract(imgH, thrHot, minarea=minDetPixHot, 
                     conv=convKerHot, deblend_nthresh=deblendThrHot, 
                     deblend_cont=deblendConHot)



    











In [31]:



    


xHot   = objHot['x'].copy()
yHot   = objHot['y'].copy()
aHot   = objHot['a'].copy()
bHot   = objHot['b'].copy()
paHot  = objHot['theta'].copy()
fluxHot = objHot['cflux'].copy()

numHot = xHot.shape[0]
print "%d objects are detected in the hot run" % numHot



    


















    






182 objects are detected in the hot run

















In [32]:



    


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('Hot-mode Detections', fontsize=25, fontweight='bold')
ax.title.set_position((0.5,1.01))

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

ax.imshow(np.arcsinh(imgH), 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))
ellsHot = getEll2Plot(objHot, radius=(aHot*5.0))

for e in ellsHot:
    ax.add_artist(e)
    e.set_clip_box(ax.bbox)
    e.set_alpha(0.8)
    e.set_edgecolor('r')
    e.set_facecolor('none')
    e.set_linewidth(1.5)

Making Object Mask Image





In [33]:



    


objMskCold = np.zeros(shape(imgC), dtype = uint8)
sep.mask_ellipse(objMskCold, xCold, yCold, aCold, bCold, paCold, r=8.0)

plt.imshow(objMskCold, cmap='gray')



    


















    
Out[33]:








<matplotlib.image.AxesImage at 0x112f34cd0>










    
























In [34]:



    


objMskHot = np.zeros(shape(imgH), dtype = uint8)
sep.mask_ellipse(objMskHot, xHot, yHot, aHot, bHot, paHot, r=4.0)

plt.imshow(objMskHot, cmap='gray')



    


















    
Out[34]:








<matplotlib.image.AxesImage at 0x11164c9d0>

Remove the Central Galaxy from the Mask Array





In [35]:



    


galX  = 400.5
galY  = 400.5 
galR  = 80.0    # Radius of the galaxy 
galQ  = 0.9     # Axis Ratio of the galaxy 
galPA = 90.0    # Position angle of the galaxy



    











In [36]:



    


objColdNew = copy.deepcopy(objCold) 
objHotNew  = copy.deepcopy(objHot) 

galCenDistCold = np.sqrt((xCold - galX)**2 + (yCold - galY)**2)  # Distance between detected objects 
galCenDistHot  = np.sqrt((xHot  - galX)**2 + (yHot  - galY)**2)  #   and the target

objColdNew = np.delete(objColdNew, np.where(galCenDistCold <= galR))  # Clear the central region 
#objColdNew = np.delete(objColdNew, np.argmin(galCenDistCold))  # Remove the one with the smallest distance 

objHotNew  = np.delete(objHotNew,  np.argmin(galCenDistHot))   # Remove the one with the smallest distance

print "Cold Mode: %d --> %d" % (objCold.shape[0], objColdNew.shape[0])
print "Hot  Mode: %d --> %d" % (objHot.shape[0],  objHotNew.shape[0])



    


















    






Cold Mode: 276 --> 264
Hot  Mode: 182 --> 181

















In [38]:



    


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('Cold-mode / No Central Galaxy', fontsize=25, fontweight='bold')
ax.title.set_position((0.5,1.01))

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

ax.imshow(np.arcsinh(imgC), 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))
ellsCold = getEll2Plot(objColdNew, radius=(objColdNew['a']*5.0))

for e in ellsCold:
    ax.add_artist(e)
    e.set_clip_box(ax.bbox)
    e.set_alpha(0.8)
    e.set_edgecolor('b')
    e.set_facecolor('none')
    e.set_linewidth(1.5)
  • Bright object around central galaxy can still be detected in the cold mode
  • To make the size of the mask large enough for most objects, the masks for objects close to central galaxy might get too large in the colde mode
  • TODO: Need a better scheme to increase the size of mask
  • Work Around: Remove all masks within certain radius of the central galaxies in the colde mode




In [42]:



    


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('Hot-mode / No Central Galaxy', fontsize=25, fontweight='bold')
ax.title.set_position((0.5,1.01))

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

ax.imshow(np.arcsinh(imgH), 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))
ellsHot = getEll2Plot(objHotNew, radius=(objHotNew['a']*4.0))

for e in ellsHot:
    ax.add_artist(e)
    e.set_clip_box(ax.bbox)
    e.set_alpha(0.8)
    e.set_edgecolor('r')
    e.set_facecolor('none')
    e.set_linewidth(1.5)

Combine the Masks





In [43]:



    


mskGrowCold = 8.0  # 

objMskColdNew = np.zeros(shape(imgC), dtype = uint8)
xColdNew  = objColdNew['x'].copy()
yColdNew  = objColdNew['y'].copy()
aColdNew  = objColdNew['a'].copy()
bColdNew  = objColdNew['b'].copy()
paColdNew = objColdNew['theta'].copy()

sep.mask_ellipse(objMskColdNew, xColdNew, yColdNew, 
                 aColdNew, bColdNew, paColdNew, r=mskGrowCold)

plt.imshow(objMskColdNew, cmap='gray')



    


















    
Out[43]:








<matplotlib.image.AxesImage at 0x116f1d210>










    
























In [44]:



    


mskGrowHot = 4.2  # 

objMskHotNew = np.zeros(shape(imgH), dtype = uint8)
xHotNew  = objHotNew['x'].copy()
yHotNew  = objHotNew['y'].copy()
aHotNew  = objHotNew['a'].copy()
bHotNew  = objHotNew['b'].copy()
paHotNew = objHotNew['theta'].copy()

sep.mask_ellipse(objMskHotNew, xHotNew, yHotNew, 
                 aHotNew, bHotNew, paHotNew, r=mskGrowHot)

plt.imshow(objMskHotNew, cmap='gray')



    


















    
Out[44]:








<matplotlib.image.AxesImage at 0x110b378d0>










    
























In [45]:



    


objMskComb = (objMskColdNew | objMskHotNew)   # Combine the two mask 

plt.imshow(objMskComb, cmap='gray')



    


















    
Out[45]:








<matplotlib.image.AxesImage at 0x10e6543d0>










    
























In [46]:



    


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 Masked Galaxy Image', fontsize=25, fontweight='bold')
ax.title.set_position((0.5,1.01))

galMsk = copy.deepcopy(img)
mmin, mmax = zscale(galMsk, contrast=0.6, samples=500)

galMsk[objMskComb > 0] = np.nan 
    
cmap = cubehelix.cmap(start=0.5, rot=-0.8, minSat=1.2, maxSat=1.2, 
                      minLight=0., maxLight=1., gamma=1.0)
cmap.set_bad('k',1.)

ax.imshow(np.arcsinh(galMsk), interpolation="none", 
          vmin=mmin, vmax=mmax, cmap=cmap)



    


















    
Out[46]:








<matplotlib.image.AxesImage at 0x110b90490>

Save the Input Image, Object Mask, and Sigma Images to Individual Files





In [54]:



    


import os
filePrefix = os.path.splitext(file)[0]

imgFile = filePrefix + '_ori.fits' 
sigFile = filePrefix + '_sig.fits' 
mskFile = filePrefix + '_msk.fits' 

# Just in case, save a "exposure time corrected" version of the image 
corFile = filePrefix + '_cor.fits' 
coaddExptime  = 900.0   # Assumed a 900s total exposure time 
coaddNcombine = 3       # NCOMBINE = 3

Update Keywords in the Headers





In [48]:



    


priHead = copy.deepcopy(data[0].header)

# Get the photometric zeropoint 
photoZP = 2.5 * np.log10(priHead['FLUXMAG0'])
print "The photometric zeropoint is %8.4f mag" % (photoZP)

priHead.remove('NAXIS')
priHead.remove('BITPIX')
priHead.remove('EXPTIME')
priHead



    


















    






The photometric zeropoint is  27.0000 mag










    
Out[48]:








SIMPLE  =                    T / file does conform to FITS standard             
EXTEND  =                    T / FITS dataset may contain extensions            
COMMENT   FITS (Flexible Image Transport System) format is defined in 'Astronomy
COMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H 
HIERARCH HSCPIPE_VERSION = '3.4.1-3-gc2ac19a'                                   
HOST    = 'analysis16.hsc.ipmu.jp'                                              
USER    = 'yasuda  '                                                            
HIERARCH ROOT_PATH = '/lustre/Subaru/SSP/rerun/yasuda/SSP3.4.1_20141224'        
AR_HDU  =                    5 / HDU containing the archive used to store ancill
HIERARCH COADD_INPUTS_ID =   1 / archive ID for coadd inputs catalogs           
HIERARCH AP_CORR_MAP_ID =  245 / archive ID for aperture correction map         
PSF_ID  =                  259 / archive ID for the Exposure's main Psf         
WCS_ID  =                  261 / archive ID for the Exposure's main Wcs         
FILTER  = 'i       '                                                            
TIME-MID= '1969-12-31T23:59:51.999918210Z'                                      
FLUXMAG0=     63095734448.0194                                                  
HIERARCH FLUXMAG0ERR =      0.                                                  

















In [49]:



    


imgHnew = copy.deepcopy(imgHead)
varHnew = copy.deepcopy(varHead)
mskHnew = copy.deepcopy(mskHead)

imgHnew.remove('EXTTYPE')
imgHnew.remove('EXTNAME')
imgHnew.remove('XTENSION')
imgHnew.set('PSCALEX', pixScaleX)
imgHnew.set('PSCALEY', pixScaleX)
imgHnew.set('GAIN', 3.0)        # This is just typical value for single HSC CCD 
imgHnew.set('EXPTIME', 1.0)     # The coadd image is normalized to 1.0s 
imgHnew.set('PHOTOZP', photoZP) # Photometric zeropoint 

varHnew.remove('EXTTYPE')
varHnew.remove('EXTNAME')
varHnew.remove('XTENSION')
varHnew.set('PSCALEX', pixScaleX)
varHnew.set('PSCALEY', pixScaleX)
varHnew.set('GAIN', 3.0)        # This is just typical value for single HSC CCD 
varHnew.set('EXPTIME', 1.0)     # The coadd image is normalized to 1.0s 
varHnew.set('PHOTOZP', photoZP) # Photometric zeropoint 

mskHnew.remove('EXTTYPE')
mskHnew.remove('EXTNAME')
mskHnew.remove('XTENSION')
mskHnew.set('PSCALEX', pixScaleX)
mskHnew.set('PSCALEY', pixScaleX)



    











In [50]:



    


priHead + imgHnew



    


















    
Out[50]:








SIMPLE  =                    T / file does conform to FITS standard             
EXTEND  =                    T / FITS dataset may contain extensions            
COMMENT   FITS (Flexible Image Transport System) format is defined in 'Astronomy
COMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H 
HIERARCH HSCPIPE_VERSION = '3.4.1-3-gc2ac19a'                                   
HOST    = 'analysis16.hsc.ipmu.jp'                                              
USER    = 'yasuda  '                                                            
HIERARCH ROOT_PATH = '/lustre/Subaru/SSP/rerun/yasuda/SSP3.4.1_20141224'        
AR_HDU  =                    5 / HDU containing the archive used to store ancill
HIERARCH COADD_INPUTS_ID =   1 / archive ID for coadd inputs catalogs           
HIERARCH AP_CORR_MAP_ID =  245 / archive ID for aperture correction map         
PSF_ID  =                  259 / archive ID for the Exposure's main Psf         
WCS_ID  =                  261 / archive ID for the Exposure's main Wcs         
FILTER  = 'i       '                                                            
TIME-MID= '1969-12-31T23:59:51.999918210Z'                                      
FLUXMAG0=     63095734448.0194                                                  
HIERARCH FLUXMAG0ERR =      0.                                                  
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                                                  
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                                       
PSCALEX =   0.1680000000000001                                                  
PSCALEY =   0.1680000000000001                                                  
GAIN    =                  3.0                                                  
EXPTIME =                  1.0                                                  
PHOTOZP =                 27.0                                                  

















In [51]:



    


imgHdu = fits.PrimaryHDU(img)
sigHdu = fits.PrimaryHDU(sig) 
mskHdu = fits.PrimaryHDU(objMskComb)



    











In [52]:



    


imgHdu.header = (priHead + imgHnew)
sigHdu.header = (priHead + varHnew)
mskHdu.header = (priHead + mskHnew)



    











In [53]:



    


imgHdu.writeto(imgFile, clobber=True)
sigHdu.writeto(sigFile, clobber=True)
mskHdu.writeto(mskFile, clobber=True)

To test GALFIT, save a "Exposure-time corrected" version of the image

  • Use reasonable EXPTIME and NCOMBINE keywords, and change the pixel unit from counts/sec to counts. 




In [56]:



    


# 
corHnew = copy.deepcopy(imgHead)
# 
imgCor = img * coaddExptime

corHnew.remove('EXTTYPE')
corHnew.remove('EXTNAME')
corHnew.remove('XTENSION')
corHnew.set('EXPTIME',  coaddExptime) 
corHnew.set('NCOMBINE', coaddNcombine)     
corHnew.set('PSCALEX', pixScaleX)
corHnew.set('PSCALEY', pixScaleX)
corHnew.set('GAIN', 3.0)        # This is just typical value for single HSC CCD 
corHnew.set('PHOTOZP', photoZP) # Photometric zeropoint 

corHdu = fits.PrimaryHDU(imgCor)
corHdu.header = (priHead + corHnew)
corHdu.writeto(corFile, clobber=True)

Automatically decide which galaxy to fit, and which to mask

  • Use level of "deblendedness" to decide?

Exam the objects detected by SEP





In [69]:



    


fig = plt.figure(figsize=(8, 8))
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_xlabel('log(Flux)',      fontsize=23)
ax.set_ylabel('Major Axis Radius (pixel)', fontsize=23)

ax.plot(np.log10(objCold['cflux']), objCold['a'], 'o', 
        color='b', label='Cold')

ax.plot(np.log10(objHot['cflux']),  objHot['a'], 'o', 
        color='r', label='Hot')

# Put legend on
ax.legend(loc=[0.06, 0.78], fontsize=20)



    


















    
Out[69]:








<matplotlib.legend.Legend at 0x118177f10>










    
























In [ ]:

Content source: dr-guangtou/hs_galphot

Similar notebooks:

notebook.community | gallery | about