In [1]:
%pylab inline
from __future__ import (division, print_function)
import os
import sys
import copy
import fnmatch
import warnings
# Numpy & Scipy
import scipy
import numpy as numpy
# Astropy related
from astropy.io import fits
from astropy import wcs
from astropy import units as u
from astropy.table import Table, Column, vstack
from astropy.stats import sigma_clip
# Matplotlib
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
from matplotlib.ticker import NullFormatter
from matplotlib.ticker import MaxNLocator
# from astroML.plotting import hist
# plt.ioff()
# ColorMap
from palettable.colorbrewer.sequential import PuBu_5, OrRd_6
cmap1 = PuBu_5.mpl_colormap
cmap2 = OrRd_6.mpl_colormap
# 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
mpl.rcParams['legend.numpoints'] = 1
rc('axes', linewidth=2)
# Shapely related imports
from shapely.geometry import Polygon, LineString, Point
from shapely import wkb
from shapely.ops import cascaded_union
from shapely.prepared import prep
from descartes import PolygonPatch
In [2]:
# SDSS pivot wavelength
sdss_u_pivot = 3551.0
sdss_g_pivot = 4686.0
sdss_r_pivot = 6165.0
sdss_i_pivot = 7481.0
sdss_z_pivot = 8931.0
# GALEX pivot wavelength
galex_fuv_pivot = 1535.0
galex_nuv_pivot = 2301.0
# WISE pivot wavelength
wise_w1_pivot = 34000.0
wise_w2_pivot = 46000.0
# HSC pivot wavelength
hsc_g_pivot = 4782.2
hsc_r_pivot = 6101.7
hsc_i_pivot = 7648.0
hsc_z_pivot = 8883.0
hsc_y_pivot = 9750.8
hscFiltWave = np.asarray([hsc_g_pivot, hsc_r_pivot, hsc_i_pivot, hsc_z_pivot, hsc_y_pivot])
In [3]:
def addMaggies(inputCat, magType='mag_cmodel',
snr=None,
filters=['g', 'r', 'i', 'z', 'y'],
maggiesCol='cmodel_maggies',
ivarsCol='cmodel_ivars',
saveNew=True, outputCat=None,
sortCol=None):
"""Convert the magnitude and error into Maggies and IvarMaggies."""
if not os.path.isfile(inputCat):
raise Exception("Can not find input catalog: %s" % s)
data = Table.read(inputCat, format='fits')
maggies = np.dstack((map(lambda f: hscMag2Flux(data[f + magType] -
data['a_' + f], unit='maggy'), filters)))[0]
if snr is None:
ivars = np.dstack((map(lambda f: hscMagerr2Ivar(
hscMag2Flux(data[f + magType] - data['a_' + f], unit='maggy'),
data[f + magType + '_err']), filters)))[0]
else:
ivars = np.dstack((map(lambda f: hscFluxSNR2Ivar(
hscMag2Flux(data[f + magType] - data['a_' + f], unit='maggy'),
snr), filters)))[0]
data.add_column(Column(name=maggiesCol, data=maggies))
data.add_column(Column(name=ivarsCol, data=ivars))
if sortCol is not None:
data.sort(sortCol)
if saveNew:
if outputCat is None:
newCat = inputCat.replace('.fits', '_' + magType + '_maggies.fits')
else:
newCat = outputCat
data.write(newCat, format='fits', overwrite=True)
return data
In [4]:
def hscFlux2AB(flux, zero=27.0):
"""
Convert HSC flux in unit of ADU to AB magnitude.
So far, constant zeropoint is applied to the calibration
"""
try:
mag = -2.5 * np.log10(flux) + zero
except NameError:
import numpy as np
mag = -2.5 * np.log10(flux) + zero
return mag
def hscMag2Flux(mag, unit='maggy'):
"""
Convert HSC AB magnitude into physical flux.
Three units can be used here:
unit='maggy/nanomaggy/jy'
"""
flux = 10.0 ** (-0.4 * mag)
if unit.lower().strip() == 'jy':
return (flux * 3631.0)
elif unit.lower().strip() == 'maggy':
return flux
elif unit.lower().strip() == 'nanomaggy':
return (flux * 1.0E-9)
else:
raise Exception("## Wrong unit, should be jy/maggy/nanomaggy")
def hscMaggy2AB(flux):
"""
Convert flux in unit of Maggies into AB magnitude
"""
return (np.log10(flux) / -0.4)
def hscMaggyErr2ABErr(flux, fluxErr, ivar=False):
"""
Convert (flux, fluxErr) into AB magnitude error
"""
if ivar:
fluxErr = np.sqrt(1.0 / fluxErr)
return (2.5 * np.log10((flux + fluxErr) / flux))
def hscMagerr2Ivar(flux, magErr):
"""
Get the inverse variance of flux estimates from Flux and magErr
"""
fluxErr = flux * ((10.0 ** (magErr/2.5)) - 1.0)
return (1.0 / (fluxErr ** 2.0))
def hscMagerr2Fluxerr(flux, magErr):
"""
Get the inverse variance of flux estimates from Flux and magErr
"""
fluxErr = flux * ((10.0 ** (magErr/2.5)) - 1.0)
return fluxErr
def hscFluxSNR2Ivar(flux, snr):
"""
Estimate inverse variance of flux error using HSC flux and snr
"""
fluxErr = flux / snr
return (1.0 / (fluxErr ** 2.0))
In [5]:
# Working directory
dataDir = '/Users/songhuang/Downloads/dr15b'
galDir = os.path.join(dataDir, 'wide_galaxy')
galWide = 'dr1_gal21_cmodel_i.fits'
mosaicDir = os.path.join(dataDir, 'basic/mosaic')
mosaicPre = 's15b_wide_i_mosaic_REG.fits'
polyDir = os.path.join(dataDir, 'basic/polymask/wide')
acpFormat = 'dr1_wide_HSC-FILTER_patch_REG.wkb'
rejFormat = 'dr1_wide_HSC-FILTER_wkbBig_REG.wkb'
# Wide fields used:
## GAMA09; GAMA15; WIDE12; XMM-LSS; HECTOMAP; VVDS
fields = ['g09', 'g15', 'w12', 'xmm', 'hec', 'vvd']
filters = ['G', 'R', 'I', 'Z', 'Y']
# Spec-z catalog
speczCat = os.path.join(dataDir, 'basic/specz/dr1_specz.fits')
# SDSS Master
sdssMaster = os.path.join(dataDir, 'sdss', 'sdss_dr12_i20.5_master.fits')
# GAMA Master
gamaMaster = os.path.join(dataDir, 'gama', 'gama_dr2_master.fits')
# redMaPPer Master
redbcgMaster = os.path.join(dataDir, 'redmapper', 'redmapper_dr8_bcg_master.fits')
redmemMaster = os.path.join(dataDir, 'redmapper', 'redmapper_dr8_mem_master.fits')
In [6]:
redbcg_old = 'redbcg_old.fits'
nonbcg_old = 'nonbcg_old.fits'
redbcg_old_b = 'redbcg_old_s15b.fits'
nonbcg_old_b = 'nonbcg_old_s15b.fits'
filter1 = ['g', 'r', 'i', 'z', 'y']
filter2 = ['g', 'r', 'i', 'z']
snr1 = 100
snr2 = 50
In [7]:
redbcg_old = os.path.join(dataDir, 'dr15a', redbcg_old)
nonbcg_old = os.path.join(dataDir, 'dr15a', nonbcg_old)
# Old results, Old cModel
## SNR=100
redOld1 = addMaggies(redbcg_old, sortCol='z_use', filters=filter1,
snr=snr1, outputCat='redbcg_old_snr100_5band.fits')
redOld2 = addMaggies(redbcg_old, sortCol='z_use', filters=filter2,
snr=snr1, outputCat='redbcg_old_snr100_4band.fits')
nonOld1 = addMaggies(nonbcg_old, sortCol='z_use', filters=filter1,
snr=snr1, outputCat='nonbcg_old_snr100_5band.fits')
nonOld2 = addMaggies(nonbcg_old, sortCol='z_use', filters=filter2,
snr=snr1, outputCat='nonbcg_old_snr100_4band.fits')
## SNR=50
redOld3 = addMaggies(redbcg_old, sortCol='z_use', filters=filter1,
snr=snr2, outputCat='redbcg_old_snr50_5band.fits')
redOld4 = addMaggies(redbcg_old, sortCol='z_use', filters=filter2,
snr=snr2, outputCat='redbcg_old_snr50_4band.fits')
nonOld3 = addMaggies(nonbcg_old, sortCol='z_use', filters=filter1,
snr=snr2, outputCat='nonbcg_old_snr50_5band.fits')
nonOld4 = addMaggies(nonbcg_old, sortCol='z_use', filters=filter2,
snr=snr2, outputCat='nonbcg_old_snr50_4band.fits')
# Old results, New cModel
redbcg_old_b = os.path.join(dataDir, 'dr15a', redbcg_old_b)
nonbcg_old_b = os.path.join(dataDir, 'dr15a', nonbcg_old_b)
## SNR=100
redOld5 = addMaggies(redbcg_old_b, sortCol='z_use', filters=filter1,
snr=snr1, outputCat='redbcg_old_b_snr100_5band.fits')
redOld6 = addMaggies(redbcg_old_b, sortCol='z_use', filters=filter2,
snr=snr1, outputCat='redbcg_old_b_snr100_4band.fits')
nonOld5 = addMaggies(nonbcg_old_b, sortCol='z_use', filters=filter1,
snr=snr1, outputCat='nonbcg_old_b_snr100_5band.fits')
nonOld6 = addMaggies(nonbcg_old_b, sortCol='z_use', filters=filter2,
snr=snr1, outputCat='nonbcg_old_b_snr100_4band.fits')
## SNR=50
redOld7 = addMaggies(redbcg_old_b, sortCol='z_use', filters=filter1,
snr=snr2, outputCat='redbcg_old_b_snr50_5band.fits')
redOld8 = addMaggies(redbcg_old_b, sortCol='z_use', filters=filter2,
snr=snr2, outputCat='redbcg_old_b_snr50_4band.fits')
nonOld7 = addMaggies(nonbcg_old_b, sortCol='z_use', filters=filter1,
snr=snr2, outputCat='nonbcg_old_b_snr50_5band.fits')
nonOld8 = addMaggies(nonbcg_old_b, sortCol='z_use', filters=filter2,
snr=snr2, outputCat='nonbcg_old_b_snr50_4band.fits')
In [ ]: