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