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


Populating the interactive namespace from numpy and matplotlib

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))

Old results from 2016-01

  • redbcg_old.fits
  • nonbcg_old.fits

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')


WARNING: UnitsWarning: 'dex' did not parse as fits unit: At col 0, Unit 'dex' not supported by the FITS standard.  [astropy.units.core]

Modify the old stellar mass estimates


In [ ]: