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
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)
In [5]:
# Information and structure of the cutout image
data.info()
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]:
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)
In [10]:
psfData.info()
imgPSF = psfData[0].data
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
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]:
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]:
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]:
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]:
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))
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')
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]:
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
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]:
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]:
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
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)
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)
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
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)
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]:
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]:
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])
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)
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)
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]:
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]:
In [45]:
objMskComb = (objMskColdNew | objMskHotNew) # Combine the two mask
plt.imshow(objMskComb, cmap='gray')
Out[45]:
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]:
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
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
Out[48]:
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]:
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)
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)
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]:
In [ ]: