notebook.community

Edit and run 





In [1]:



    


%matplotlib inline
%reload_ext autoreload
from __future__ import (division, 
                        print_function)

import os
import sys
import copy
import glob
import fnmatch
import warnings
import collections

import numpy as np
import scipy
try:
    from scipy.stats import scoreatpercentile
except:
    scoreatpercentile = False
from scipy.interpolate import interp1d
from scipy.ndimage import gaussian_filter
import cPickle as pickle

# Astropy
from astropy.io import fits
from astropy    import units as u
from astropy.stats import sigma_clip
from astropy.table import Table, Column, join
from astropy.utils.console import ProgressBar
from astropy.convolution import convolve, Box1DKernel

# AstroML
from astroML.plotting import hist
from astroML.density_estimation import KNeighborsDensity
try:
    from sklearn.neighbors import KernelDensity
    use_sklearn_KDE = True
except:
    import warnings
    warnings.warn("KDE will be removed in astroML version 0.3.  Please "
                  "upgrade to scikit-learn 0.14+ and use "
                  "sklearn.neighbors.KernelDensity.", DeprecationWarning)
    from astroML.density_estimation import KDE
    use_sklearn_KDE = False
from sklearn.neighbors import KDTree
from sklearn.neighbors import BallTree

# Matplotlib related
import matplotlib as mpl
mpl.style.use('classic')
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib.patches import Ellipse
from matplotlib.ticker import NullFormatter, MaxNLocator, FormatStrFormatter
from matplotlib.collections import PatchCollection

# Personal
import hscUtils as hUtil
#import galSBP
#import coaddCutoutGalfitSimple as gSimple 

# Cosmology
import cosmology
c=cosmology.Cosmo(H0=70.0, omega_m=0.3, omega_l=0.7, flat=1)

import emcee
import corner

# Color map 
from palettable.colorbrewer.sequential import Greys_9, OrRd_9, Blues_9, Purples_9, YlGn_9
BLK = Greys_9.mpl_colormap
ORG = OrRd_9.mpl_colormap
BLU = Blues_9.mpl_colormap
GRN = YlGn_9.mpl_colormap
PUR = Purples_9.mpl_colormap

# Personal tools
from hscUtils import songPlotSetup, removeIsNullCol
from hscUtils import confidence_interval, ma_confidence_interval_1d, confidence_interval_1d

## Constants
# 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

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

"""
Absolute magnitude of the Sun in HSC filters
Right now, just use the DES filters
"""
SUN_G = 5.08
SUN_R = 4.62
SUN_I = 4.52
SUN_Z = 4.52
SUN_Y = 4.51

# Solar stellar metallicity 
Z_SUN = 0.02

# definitions for the axes
left, width    = 0.15, 0.64
right          = left + width 
bottom, height = 0.13, 0.86
bottom_h = left_h = left + width + 0.02
recScat = [0.16, 0.11, 0.59, 0.88]
recHist = [0.75, 0.11, 0.24, 0.88]
SBP1 = [0.124, 0.085, 0.865, 0.33]
SBP2 = [0.124, 0.41, 0.865, 0.55]
EC1 = [0.19, 0.14, 0.65, 0.65]
EC2 = [0.19, 0.79, 0.65, 0.18]
EC3 = [0.84, 0.14, 0.157, 0.65]
REC = [0.12, 0.11, 0.87, 0.87]
COG1 = [0.143, 0.10, 0.850, 0.43]
COG2 = [0.143, 0.53, 0.850, 0.43]

# Universal RSMA array
RSMA_COMMON = np.arange(0.4, 4.2, 0.01)
RR50_COMMON = np.arange(0.0, 9.0, 0.01)
EMPTY = (RSMA_COMMON * np.nan)

# Color 
BLUE0 = "#92c5de"
BLUE1 = "#0571b0"
RED0 = "#f4a582"
RED1 = "#ca0020"
PURPLE0 = '#af8dc3'
PURPLE1 = '#762a83'
BROWN0 = '#bf812d'
BROWN1 = '#543005'
GREEN0 = '#7fbf7b'
GREEN1 = '#1b7837'

# Color maps for different filters
from palettable.colorbrewer.sequential import Blues_9, Greens_9, Reds_9, Purples_9, Greys_9
GCMAP = Blues_9.mpl_colormap
RCMAP = Greens_9.mpl_colormap
ICMAP = Reds_9.mpl_colormap
ZCMAP = Purples_9.mpl_colormap
YCMAP = Greys_9.mpl_colormap

# 3-sigma
SIGMA1 = 0.3173
SIGMA2 = 0.0455
SIGMA3 = 0.0027

#import seaborn as sns
#sns.set(color_codes=False)

Location and Files





In [2]:



    


# Location of the data
homeDir = os.getenv('HOME')
synpipeDir = os.path.join(homeDir, 'astro4/synpipe/')
sampleDir = os.path.join(synpipeDir, 'sample/cosmos')

# Samples
cosPhoFile = os.path.join(sampleDir, 'cosmos_mag25.2_shuang_photo.fits')
cosHscFile = os.path.join(sampleDir, 'cosmos_mag25.2_shuang_photo_hscmatch.fits')

cosPho = Table.read(cosPhoFile, format='fits')
cosHsc = Table.read(cosHscFile, format='fits')

# Outputs 
figDir = os.path.join(synpipeDir, 'figure')
outDir = os.path.join(synpipeDir, 'outputs')
inDir = os.path.join(synpipeDir, 'inputs')

starDir = os.path.join(outDir, 'star')
galaxyDir = os.path.join(outDir, 'galaxy')

Test outputs for stars





In [3]:



    


glob.glob(inDir + '/*gal*')



    


















    
Out[3]:








['/Users/song/astro4/synpipe/inputs/galaxy_all_9699_HSC-Y.fits',
 '/Users/song/astro4/synpipe/inputs/galaxy_all_8764_HSC-G.fits',
 '/Users/song/astro4/synpipe/inputs/galaxy_all_8764_HSC-I.fits',
 '/Users/song/astro4/synpipe/inputs/galaxy_all_8764_HSC-R.fits',
 '/Users/song/astro4/synpipe/inputs/galaxy_all_8764_HSC-Y.fits',
 '/Users/song/astro4/synpipe/inputs/galaxy_all_8764_HSC-Z.fits',
 '/Users/song/astro4/synpipe/inputs/galaxy_all_9699_HSC-G.fits',
 '/Users/song/astro4/synpipe/inputs/galaxy_all_9699_HSC-I.fits',
 '/Users/song/astro4/synpipe/inputs/galaxy_all_9699_HSC-R.fits',
 '/Users/song/astro4/synpipe/inputs/galaxy_all_9699_HSC-Z.fits']

















In [4]:



    


glob.glob(galaxyDir + '/*syn*fits')



    


















    
Out[4]:








['/Users/song/astro4/synpipe/outputs/galaxy/gal_all1_HSC-G_syn.fits',
 '/Users/song/astro4/synpipe/outputs/galaxy/gal_all1_HSC-G_syn_multipix.fits',
 '/Users/song/astro4/synpipe/outputs/galaxy/gal_all1_HSC-Z_syn.fits',
 '/Users/song/astro4/synpipe/outputs/galaxy/gal_all1_HSC-Z_syn_multipix.fits',
 '/Users/song/astro4/synpipe/outputs/galaxy/gal_all2_HSC-G_syn.fits',
 '/Users/song/astro4/synpipe/outputs/galaxy/gal_all2_HSC-R_syn.fits',
 '/Users/song/astro4/synpipe/outputs/galaxy/gal_all2_HSC-Y_syn.fits',
 '/Users/song/astro4/synpipe/outputs/galaxy/gal_all2_HSC-Z_syn.fits',
 '/Users/song/astro4/synpipe/outputs/galaxy/gal_all1_HSC-R_syn.fits',
 '/Users/song/astro4/synpipe/outputs/galaxy/gal_all1_HSC-R_syn_multipix.fits',
 '/Users/song/astro4/synpipe/outputs/galaxy/gal_all2_HSC-I_syn.fits',
 '/Users/song/astro4/synpipe/outputs/galaxy/gal_all1_HSC-I_syn.fits',
 '/Users/song/astro4/synpipe/outputs/galaxy/gal_all1_HSC-I_syn_multipix.fits',
 '/Users/song/astro4/synpipe/outputs/galaxy/gal_all1_HSC-Y_syn.fits',
 '/Users/song/astro4/synpipe/outputs/galaxy/gal_all1_HSC-Y_syn_multipix.fits']

Input catalogs





In [13]:



    


inputG1 = Table.read(os.path.join(inDir, 'galaxy_all_8764_HSC-G.fits'), format='fits')
inputG2 = Table.read(os.path.join(inDir, 'galaxy_all_9699_HSC-G.fits'), format='fits')

inputR1 = Table.read(os.path.join(inDir, 'galaxy_all_8764_HSC-R.fits'), format='fits')
inputR2 = Table.read(os.path.join(inDir, 'galaxy_all_9699_HSC-R.fits'), format='fits')

inputI1 = Table.read(os.path.join(inDir, 'galaxy_all_8764_HSC-I.fits'), format='fits')
inputI2 = Table.read(os.path.join(inDir, 'galaxy_all_9699_HSC-I.fits'), format='fits')

inputZ1 = Table.read(os.path.join(inDir, 'galaxy_all_8764_HSC-Z.fits'), format='fits')
inputZ2 = Table.read(os.path.join(inDir, 'galaxy_all_9699_HSC-Z.fits'), format='fits')

inputY1 = Table.read(os.path.join(inDir, 'galaxy_all_8764_HSC-Y.fits'), format='fits')
inputY2 = Table.read(os.path.join(inDir, 'galaxy_all_9699_HSC-Y.fits'), format='fits')



    











In [5]:



    


print(inputG1.colnames)



    


















    






['RA', 'Dec', 'modelID', 'sersic_n', 'reff', 'b_a', 'theta', 'mag_I', 'mag_G', 'mag_R', 'mag_Z', 'mag_Y', 'ID', 'mag']

Matches with the Synthetic-added images





In [3]:



    


galG1 = os.path.join(galaxyDir, 'gal_all1_HSC-G_syn.fits') 
galR1 = os.path.join(galaxyDir, 'gal_all1_HSC-R_syn.fits') 
galI1 = os.path.join(galaxyDir, 'gal_all1_HSC-I_syn.fits') 
galZ1 = os.path.join(galaxyDir, 'gal_all1_HSC-Z_syn.fits') 
galY1 = os.path.join(galaxyDir, 'gal_all1_HSC-Y_syn.fits') 

galG1m = os.path.join(galaxyDir, 'gal_all1_HSC-G_syn_multipix.fits') 
galR1m = os.path.join(galaxyDir, 'gal_all1_HSC-R_syn_multipix.fits') 
galI1m = os.path.join(galaxyDir, 'gal_all1_HSC-I_syn_multipix.fits') 
galZ1m = os.path.join(galaxyDir, 'gal_all1_HSC-Z_syn_multipix.fits') 
galY1m = os.path.join(galaxyDir, 'gal_all1_HSC-Y_syn_multipix.fits') 

galG2 = os.path.join(galaxyDir, 'gal_all2_HSC-G_syn.fits') 
galR2 = os.path.join(galaxyDir, 'gal_all2_HSC-R_syn.fits') 
galI2 = os.path.join(galaxyDir, 'gal_all2_HSC-I_syn.fits') 
galZ2 = os.path.join(galaxyDir, 'gal_all2_HSC-Z_syn.fits') 
galY2 = os.path.join(galaxyDir, 'gal_all2_HSC-Y_syn.fits')



    











In [4]:



    


catG1 = Table.read(galG1, format='fits')
catR1 = Table.read(galR1, format='fits')
catI1 = Table.read(galI1, format='fits')
catZ1 = Table.read(galZ1, format='fits')
catY1 = Table.read(galY1, format='fits')

catG1m = Table.read(galG1m, format='fits')
catR1m = Table.read(galR1m, format='fits')
catI1m = Table.read(galI1m, format='fits')
catZ1m = Table.read(galZ1m, format='fits')
catY1m = Table.read(galY1m, format='fits')

catG2 = Table.read(galG2, format='fits')
catR2 = Table.read(galR2, format='fits')
catI2 = Table.read(galI2, format='fits')
catZ2 = Table.read(galZ2, format='fits')
catY2 = Table.read(galY2, format='fits')

Matches with the Original images





In [6]:



    


oriG1 = os.path.join(galaxyDir, 'gal_all1_HSC-G_ori.fits') 
oriR1 = os.path.join(galaxyDir, 'gal_all1_HSC-R_ori.fits') 
oriI1 = os.path.join(galaxyDir, 'gal_all1_HSC-I_ori.fits') 
oriZ1 = os.path.join(galaxyDir, 'gal_all1_HSC-Z_ori.fits') 
oriY1 = os.path.join(galaxyDir, 'gal_all1_HSC-Y_ori.fits') 

oriG2 = os.path.join(galaxyDir, 'gal_all2_HSC-G_ori.fits') 
oriR2 = os.path.join(galaxyDir, 'gal_all2_HSC-R_ori.fits') 
oriI2 = os.path.join(galaxyDir, 'gal_all2_HSC-I_ori.fits') 
oriZ2 = os.path.join(galaxyDir, 'gal_all2_HSC-Z_ori.fits') 
oriY2 = os.path.join(galaxyDir, 'gal_all2_HSC-Y_ori.fits')



    











In [7]:



    


oriG1 = Table.read(oriG1, format='fits')
oriR1 = Table.read(oriR1, format='fits')
oriI1 = Table.read(oriI1, format='fits')
oriZ1 = Table.read(oriZ1, format='fits')
oriY1 = Table.read(oriY1, format='fits')

oriG2 = Table.read(oriG2, format='fits')
oriR2 = Table.read(oriR2, format='fits')
oriI2 = Table.read(oriI2, format='fits')
oriZ2 = Table.read(oriZ2, format='fits')
oriY2 = Table.read(oriY2, format='fits')

Useful columns





In [12]:



    


print("# Number of columns : %d" % len(catG1.colnames))
print("\n# Columns related to flux")
print([col for col in catG1.colnames if 'flux' in col])
print("\n# Columns related to mag")
print([col for col in catG1.colnames if 'mag' in col])
print("\n# Columns related to coord")
print([col for col in catG1.colnames if 'coord' in col])
print("\n# Columns related to fake")
print([col for col in catG1.colnames if 'fake' in col.lower()])
print("\n# Columns related to id")
print([col for col in catG1.colnames if 'id' in col.lower()])
print("\n# Columns related to blend")
print([col for col in catG1.colnames if 'blend' in col.lower()])
print("\n# Columns related to flag")
print([col for col in catG1.colnames if 'flag' in col.lower()])



    


















    






# Number of columns : 360

# Columns related to flux
['blendedness.abs.flux', 'blendedness.abs.flux.child', 'blendedness.abs.flux.parent', 'blendedness.raw.flux', 'blendedness.raw.flux.child', 'blendedness.raw.flux.parent', 'cmodel.dev.flux', 'cmodel.dev.flux.apcorr', 'cmodel.dev.flux.apcorr.err', 'cmodel.dev.flux.err', 'cmodel.dev.flux.flags', 'cmodel.dev.flux.flags.apcorr', 'cmodel.dev.flux.inner', 'cmodel.exp.flux', 'cmodel.exp.flux.apcorr', 'cmodel.exp.flux.apcorr.err', 'cmodel.exp.flux.err', 'cmodel.exp.flux.flags', 'cmodel.exp.flux.flags.apcorr', 'cmodel.exp.flux.inner', 'cmodel.flux', 'cmodel.flux.apcorr', 'cmodel.flux.apcorr.err', 'cmodel.flux.err', 'cmodel.flux.flags', 'cmodel.flux.flags.apcorr', 'cmodel.flux.inner', 'cmodel.initial.flux', 'cmodel.initial.flux.err', 'cmodel.initial.flux.flags', 'cmodel.initial.flux.inner', 'deblend.has.stray.flux', 'deblend.psf-flux', 'flux.aperture', 'flux.aperture.err', 'flux.aperture.flags', 'flux.aperture.nInterpolatedPixel', 'flux.aperture.nProfile', 'flux.gaussian', 'flux.gaussian.apcorr', 'flux.gaussian.apcorr.err', 'flux.gaussian.err', 'flux.gaussian.flags', 'flux.gaussian.flags.apcorr', 'flux.kron', 'flux.kron.apcorr', 'flux.kron.apcorr.err', 'flux.kron.err', 'flux.kron.flags', 'flux.kron.flags.apcorr', 'flux.kron.flags.badShape', 'flux.kron.flags.edge', 'flux.kron.flags.radius', 'flux.kron.flags.smallRadius', 'flux.kron.flags.usedMinimumRadius', 'flux.kron.flags.usedPsfRadius', 'flux.kron.psfRadius', 'flux.kron.radius', 'flux.kron.radiusForRadius', 'flux.naive', 'flux.naive.err', 'flux.naive.flags', 'flux.psf', 'flux.psf.apcorr', 'flux.psf.apcorr.err', 'flux.psf.err', 'flux.psf.flags', 'flux.psf.flags.apcorr', 'flux.sinc', 'flux.sinc.err', 'flux.sinc.flags', 'force.cmodel.dev.flux', 'force.cmodel.dev.flux.apcorr', 'force.cmodel.dev.flux.apcorr.err', 'force.cmodel.dev.flux.err', 'force.cmodel.exp.flux', 'force.cmodel.exp.flux.apcorr', 'force.cmodel.exp.flux.apcorr.err', 'force.cmodel.exp.flux.err', 'force.cmodel.flux', 'force.cmodel.flux.apcorr', 'force.cmodel.flux.apcorr.err', 'force.cmodel.flux.err', 'force.flux.kron', 'force.flux.kron.apcorr', 'force.flux.kron.apcorr.err', 'force.flux.kron.err', 'force.flux.psf', 'force.flux.psf.apcorr', 'force.flux.psf.apcorr.err', 'force.flux.psf.err']

# Columns related to mag
['flags.pixel.offimage', 'cmodel.dev.mag', 'cmodel.dev.mag.err', 'cmodel.dev.mag.apcorr', 'cmodel.dev.mag.apcorr.err', 'cmodel.exp.mag', 'cmodel.exp.mag.err', 'cmodel.exp.mag.apcorr', 'cmodel.exp.mag.apcorr.err', 'cmodel.mag', 'cmodel.mag.err', 'cmodel.mag.apcorr', 'cmodel.mag.apcorr.err', 'cmodel.initial.mag', 'cmodel.initial.mag.err', 'mag.aperture', 'mag.aperture.err', 'mag.gaussian', 'mag.gaussian.err', 'mag.gaussian.apcorr', 'mag.gaussian.apcorr.err', 'mag.kron', 'mag.kron.err', 'mag.kron.apcorr', 'mag.kron.apcorr.err', 'mag.naive', 'mag.naive.err', 'mag.psf', 'mag.psf.err', 'mag.psf.apcorr', 'mag.psf.apcorr.err', 'mag.sinc', 'mag.sinc.err', 'force.cmodel.dev.mag', 'force.cmodel.dev.mag.err', 'force.cmodel.dev.mag.apcorr', 'force.cmodel.dev.mag.apcorr.err', 'force.cmodel.exp.mag', 'force.cmodel.exp.mag.err', 'force.cmodel.exp.mag.apcorr', 'force.cmodel.exp.mag.apcorr.err', 'force.cmodel.mag', 'force.cmodel.mag.err', 'force.cmodel.mag.apcorr', 'force.cmodel.mag.apcorr.err', 'force.mag.kron', 'force.mag.kron.err', 'force.mag.kron.apcorr', 'force.mag.kron.apcorr.err', 'force.mag.psf', 'force.mag.psf.err', 'force.mag.psf.apcorr', 'force.mag.psf.apcorr.err', 'mag_I', 'mag_G', 'mag_R', 'mag_Z', 'mag_Y', 'mag']

# Columns related to coord
['coord_ra', 'coord_dec']

# Columns related to fake
['fakeClosest', 'fakeId', 'fakeOffR', 'fakeOffX', 'fakeOffY']

# Columns related to id
['blendedness.flags.noCentroid', 'centroid.distorted', 'centroid.naive', 'centroid.naive.err', 'centroid.naive.flags', 'centroid.sdss', 'centroid.sdss.err', 'centroid.sdss.flags', 'cmodel.flags.badCentroid', 'fakeId', 'flags.badcentroid', 'id', 'shape.hsm.moments.centroid', 'shape.hsm.moments.centroid.err', 'shape.hsm.moments.centroid.flags', 'shape.hsm.psfMoments.centroid', 'shape.hsm.psfMoments.centroid.err', 'shape.hsm.psfMoments.centroid.flags', 'shape.sdss.centroid', 'shape.sdss.centroid.err', 'shape.sdss.centroid.flags', 'modelID']

# Columns related to blend
['blendedness.abs.flux', 'blendedness.abs.flux.child', 'blendedness.abs.flux.parent', 'blendedness.abs.shape.child_a', 'blendedness.abs.shape.child_q', 'blendedness.abs.shape.child_theta', 'blendedness.abs.shape.parent_a', 'blendedness.abs.shape.parent_q', 'blendedness.abs.shape.parent_theta', 'blendedness.flags', 'blendedness.flags.noCentroid', 'blendedness.flags.noShape', 'blendedness.old', 'blendedness.raw.flux', 'blendedness.raw.flux.child', 'blendedness.raw.flux.parent', 'blendedness.raw.shape.child_a', 'blendedness.raw.shape.child_q', 'blendedness.raw.shape.child_theta', 'blendedness.raw.shape.parent_a', 'blendedness.raw.shape.parent_q', 'blendedness.raw.shape.parent_theta', 'deblend.deblended-as-psf', 'deblend.has.stray.flux', 'deblend.masked', 'deblend.nchild', 'deblend.parent-too-big', 'deblend.patched.template', 'deblend.psf-center', 'deblend.psf-flux', 'deblend.ramped.template', 'deblend.skipped', 'deblend.too-many-peaks', 'force.deblend.nchild']

# Columns related to flag
['blendedness.flags', 'blendedness.flags.noCentroid', 'blendedness.flags.noShape', 'centroid.naive.flags', 'centroid.sdss.flags', 'cmodel.dev.flags.maxIter', 'cmodel.dev.flags.numericError', 'cmodel.dev.flags.trSmall', 'cmodel.dev.flux.flags', 'cmodel.dev.flux.flags.apcorr', 'cmodel.exp.flags.maxIter', 'cmodel.exp.flags.numericError', 'cmodel.exp.flags.trSmall', 'cmodel.exp.flux.flags', 'cmodel.exp.flux.flags.apcorr', 'cmodel.flags.badCentroid', 'cmodel.flags.noCalib', 'cmodel.flags.noPsf', 'cmodel.flags.noShape', 'cmodel.flags.noWcs', 'cmodel.flags.region.maxArea', 'cmodel.flags.region.maxBadPixelFraction', 'cmodel.flags.region.usedFootprintArea', 'cmodel.flags.region.usedInitialEllipseMax', 'cmodel.flags.region.usedInitialEllipseMin', 'cmodel.flags.region.usedPsfArea', 'cmodel.flags.smallShape', 'cmodel.flux.flags', 'cmodel.flux.flags.apcorr', 'cmodel.initial.flags.maxIter', 'cmodel.initial.flags.numericError', 'cmodel.initial.flags.trSmall', 'cmodel.initial.flux.flags', 'flags.badcentroid', 'flags.negative', 'flags.pixel.bad', 'flags.pixel.bright.object.any', 'flags.pixel.bright.object.center', 'flags.pixel.clipped.any', 'flags.pixel.cr.any', 'flags.pixel.cr.center', 'flags.pixel.edge', 'flags.pixel.interpolated.any', 'flags.pixel.interpolated.center', 'flags.pixel.offimage', 'flags.pixel.saturated.any', 'flags.pixel.saturated.center', 'flags.pixel.suspect.any', 'flags.pixel.suspect.center', 'flux.aperture.flags', 'flux.gaussian.flags', 'flux.gaussian.flags.apcorr', 'flux.kron.flags', 'flux.kron.flags.apcorr', 'flux.kron.flags.badShape', 'flux.kron.flags.edge', 'flux.kron.flags.radius', 'flux.kron.flags.smallRadius', 'flux.kron.flags.usedMinimumRadius', 'flux.kron.flags.usedPsfRadius', 'flux.naive.flags', 'flux.psf.flags', 'flux.psf.flags.apcorr', 'flux.sinc.flags', 'multishapelet.psf.flags', 'multishapelet.psf.flags.constraint.q', 'multishapelet.psf.flags.constraint.r', 'multishapelet.psf.flags.maxiter', 'multishapelet.psf.flags.tinystep', 'shape.hsm.moments.centroid.flags', 'shape.hsm.moments.flags', 'shape.hsm.psfMoments.centroid.flags', 'shape.hsm.psfMoments.flags', 'shape.hsm.regauss.flags', 'shape.sdss.centroid.flags', 'shape.sdss.flags', 'shape.sdss.flags.maxiter', 'shape.sdss.flags.psf', 'shape.sdss.flags.shift', 'shape.sdss.flags.unweighted', 'shape.sdss.flags.unweightedbad']

Selection of Useful Galaxies

Remove the ambiguously blended galaxies

Generate new catalogs





In [8]:



    


# HSC-G
print("# HSC-G band")
ambG1 = oriG1[oriG1['id'] > 0]
ambG2 = oriG2[oriG2['id'] > 0]
print("## galaxy1: %d fake galaxies got ambiguously blended with real objects" % len(ambG1))
print("## galaxy2: %d fake galaxies got ambiguously blended with real objects" % len(ambG2))

# HSC-R
print("\n# HSC-R band")
ambR1 = oriR1[oriR1['id'] > 0]
ambR2 = oriR2[oriR2['id'] > 0]
print("## galaxy1: %d fake galaxies got ambiguously blended with real objects" % len(ambR1))
print("## galaxy2: %d fake galaxies got ambiguously blended with real objects" % len(ambR2))

# HSC-I
print("\n# HSC-I band")
ambI1 = oriI1[oriI1['id'] > 0]
ambI2 = oriI2[oriI2['id'] > 0]
print("## galaxy1: %d fake galaxies got ambiguously blended with real objects" % len(ambI1))
print("## galaxy2: %d fake galaxies got ambiguously blended with real objects" % len(ambI2))

# HSC-Z
print("\n# HSC-Z band")
ambZ1 = oriZ1[oriZ1['id'] > 0]
ambZ2 = oriZ2[oriZ2['id'] > 0]
print("## galaxy1: %d fake galaxies got ambiguously blended with real objects" % len(ambZ1))
print("## galaxy2: %d fake galaxies got ambiguously blended with real objects" % len(ambZ2))

# HSC-Y
print("\n# HSC-Y band")
ambY1 = oriY1[oriY1['id'] > 0]
ambY2 = oriY2[oriY2['id'] > 0]
print("## galaxy1: %d fake galaxies got ambiguously blended with real objects" % len(ambY1))
print("## galaxy2: %d fake galaxies got ambiguously blended with real objects" % len(ambY2))



    


















    






# HSC-G band
## galaxy1: 1549 fake galaxies got ambiguously blended with real objects
## galaxy2: 1279 fake galaxies got ambiguously blended with real objects

# HSC-R band
## galaxy1: 1516 fake galaxies got ambiguously blended with real objects
## galaxy2: 1262 fake galaxies got ambiguously blended with real objects

# HSC-I band
## galaxy1: 1517 fake galaxies got ambiguously blended with real objects
## galaxy2: 1266 fake galaxies got ambiguously blended with real objects

# HSC-Z band
## galaxy1: 1499 fake galaxies got ambiguously blended with real objects
## galaxy2: 1254 fake galaxies got ambiguously blended with real objects

# HSC-Y band
## galaxy1: 1544 fake galaxies got ambiguously blended with real objects
## galaxy2: 1217 fake galaxies got ambiguously blended with real objects

















In [28]:



    


useG1 = copy.deepcopy(catG1)
useG2 = copy.deepcopy(catG2)
useG1.remove_rows(np.nonzero(np.in1d(useG1['fakeId'], ambG1['fakeId']))[0])
useG2.remove_rows(np.nonzero(np.in1d(useG2['fakeId'], ambG2['fakeId']))[0])

useR1 = copy.deepcopy(catR1)
useR2 = copy.deepcopy(catR2)
useR1.remove_rows(np.nonzero(np.in1d(useR1['fakeId'], ambR1['fakeId']))[0])
useR2.remove_rows(np.nonzero(np.in1d(useR2['fakeId'], ambR2['fakeId']))[0])

useI1 = copy.deepcopy(catI1)
useI2 = copy.deepcopy(catI2)
useI1.remove_rows(np.nonzero(np.in1d(useI1['fakeId'], ambI1['fakeId']))[0])
useI2.remove_rows(np.nonzero(np.in1d(useI2['fakeId'], ambI2['fakeId']))[0])

useZ1 = copy.deepcopy(catZ1)
useZ2 = copy.deepcopy(catZ2)
useZ1.remove_rows(np.nonzero(np.in1d(useZ1['fakeId'], ambZ1['fakeId']))[0])
useZ2.remove_rows(np.nonzero(np.in1d(useZ2['fakeId'], ambZ2['fakeId']))[0])

useY1 = copy.deepcopy(catY1)
useY2 = copy.deepcopy(catY2)
useY1.remove_rows(np.nonzero(np.in1d(useY1['fakeId'], ambY1['fakeId']))[0])
useY2.remove_rows(np.nonzero(np.in1d(useY2['fakeId'], ambY2['fakeId']))[0])



    











In [13]:



    


useG1m = copy.deepcopy(catG1m)
useG1m.remove_rows(np.nonzero(np.in1d(useG1m['fakeId'], ambG1['fakeId']))[0])

useR1m = copy.deepcopy(catR1m)
useR1m.remove_rows(np.nonzero(np.in1d(useR1m['fakeId'], ambR1['fakeId']))[0])

useI1m = copy.deepcopy(catI1m)
useI1m.remove_rows(np.nonzero(np.in1d(useI1m['fakeId'], ambI1['fakeId']))[0])

useZ1m = copy.deepcopy(catZ1m)
useZ1m.remove_rows(np.nonzero(np.in1d(useZ1m['fakeId'], ambZ1['fakeId']))[0])

useY1m = copy.deepcopy(catY1m)
useY1m.remove_rows(np.nonzero(np.in1d(useY1m['fakeId'], ambY1['fakeId']))[0])

Write the tables





In [30]:



    


useG1.write(os.path.join(galaxyDir, 'gal_all1_HSC-G_use.fits'), format='fits', overwrite=True)
useG2.write(os.path.join(galaxyDir, 'gal_all2_HSC-G_use.fits'), format='fits', overwrite=True)

useR1.write(os.path.join(galaxyDir, 'gal_all1_HSC-R_use.fits'), format='fits', overwrite=True)
useR2.write(os.path.join(galaxyDir, 'gal_all2_HSC-R_use.fits'), format='fits', overwrite=True)

useI1.write(os.path.join(galaxyDir, 'gal_all1_HSC-I_use.fits'), format='fits', overwrite=True)
useI2.write(os.path.join(galaxyDir, 'gal_all2_HSC-I_use.fits'), format='fits', overwrite=True)

useZ1.write(os.path.join(galaxyDir, 'gal_all1_HSC-Z_use.fits'), format='fits', overwrite=True)
useZ2.write(os.path.join(galaxyDir, 'gal_all2_HSC-Z_use.fits'), format='fits', overwrite=True)

useY1.write(os.path.join(galaxyDir, 'gal_all1_HSC-Y_use.fits'), format='fits', overwrite=True)
useY2.write(os.path.join(galaxyDir, 'gal_all2_HSC-Y_use.fits'), format='fits', overwrite=True)



    











In [31]:



    


useG1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-G_use_multipix.fits'), format='fits', overwrite=True)
useR1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-R_use_multipix.fits'), format='fits', overwrite=True)
useI1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-I_use_multipix.fits'), format='fits', overwrite=True)
useZ1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-Z_use_multipix.fits'), format='fits', overwrite=True)
useY1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-Y_use_multipix.fits'), format='fits', overwrite=True)

Read the available tables





In [5]:



    


useG1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-G_use.fits'), format='fits')
useG2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-G_use.fits'), format='fits')

useR1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-R_use.fits'), format='fits')
useR2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-R_use.fits'), format='fits')

useI1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-I_use.fits'), format='fits')
useI2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-I_use.fits'), format='fits')

useZ1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-Z_use.fits'), format='fits')
useZ2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-Z_use.fits'), format='fits')

useY1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-Y_use.fits'), format='fits')
useY2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-Y_use.fits'), format='fits')



    











In [6]:



    


useG1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-G_use_multipix.fits'), format='fits')
useR1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-R_use_multipix.fits'), format='fits')
useI1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-I_use_multipix.fits'), format='fits')
useZ1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-Z_use_multipix.fits'), format='fits')
useY1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-Y_use_multipix.fits'), format='fits')

Select matched ones, isolate the unmatched ones

Generate new catalogs





In [14]:



    


# Galaxy1 
print("# HSC-G band")
print("# Galaxy1 catalog has %d objects" % len(useG1))

# Matched ones 
matchG1 = useG1[useG1['id'] > 0]
unMatchG1 = useG1[useG1['id'] == 0]
print("\n## %d synthetic galaxies are matched" % len(matchG1))
print("## %d synthetic galaxies are not matched" % len(unMatchG1))
# Median I-band magntiudes of the unmatched and matched ones
print("## Median magnitude for unmatched and matched galaxies: %5.3f, %5.3f" % (np.median(unMatchG1['mag_I']), 
                                                                                np.median(matchG1['mag_I'])))

# Matched ones 
print("\n# Multipix ones!")
matchG1m = useG1m[useG1m['id'] > 0]
unMatchG1m = useG1m[useG1m['id'] == 0]
print("## %d synthetic galaxies are matched" % len(matchG1m))
print("## %d synthetic galaxies are not matched" % len(unMatchG1m))
# Median I-band magntiudes of the unmatched and matched ones
print("## Median magnitude for unmatched and matched galaxies: %5.3f, %5.3f" % (np.median(unMatchG1m['mag_I']), 
                                                                                np.median(matchG1m['mag_I'])))

# Galaxy2
print("\n #Galaxy2 catalog has %d objects" % len(useG2))

# Matched ones 
matchG2 = useG2[useG2['id'] > 0]
unMatchG2 = useG2[useG2['id'] == 0]
print("## %d synthetic galaxies are matched" % len(matchG2))
print("## %d synthetic galaxies are not matched" % len(unMatchG2))
# Median I-band magntiudes of the unmatched and matched ones
print("## Median magnitude for unmatched and matched galaxies: %5.3f, %5.3f" % (np.median(unMatchG2['mag_I']), 
                                                                                np.median(matchG2['mag_I'])))



    


















    






# HSC-G band
# Galaxy1 catalog has 54559 objects

## 52463 synthetic galaxies are matched
## 2096 synthetic galaxies are not matched
## Median magnitude for unmatched and matched galaxies: 24.579, 24.021

# Multipix ones!
## 52463 synthetic galaxies are matched
## 2095 synthetic galaxies are not matched
## Median magnitude for unmatched and matched galaxies: 24.580, 24.021

 #Galaxy2 catalog has 54587 objects
## 51630 synthetic galaxies are matched
## 2957 synthetic galaxies are not matched
## Median magnitude for unmatched and matched galaxies: 24.606, 23.999

















In [15]:



    


# Galaxy1 
print("# HSC-R band")
print("# Galaxy1 catalog has %d objects" % len(useR1))

# Matched ones 
matchR1 = useR1[useR1['id'] > 0]
unMatchR1 = useR1[useR1['id'] == 0]
print("\n## %d synthetic galaxies are matched" % len(matchR1))
print("## %d synthetic galaxies are not matched" % len(unMatchR1))
# Median I-band magntiudes of the unmatched and matched ones
print("## Median magnitude for unmatched and matched galaxies: %5.3f, %5.3f" % (np.median(unMatchR1['mag_I']), 
                                                                                np.median(matchR1['mag_I'])))

# Matched ones 
print("\n# Multipix ones!")
matchR1m = useR1m[useR1m['id'] > 0]
unMatchR1m = useR1m[useR1m['id'] == 0]
print("\n## %d synthetic galaxies are matched" % len(matchR1m))
print("## %d synthetic galaxies are not matched" % len(unMatchR1m))
# Median I-band magntiudes of the unmatched and matched ones
print("## Median magnitude for unmatched and matched galaxies: %5.3f, %5.3f" % (np.median(unMatchR1m['mag_I']), 
                                                                                np.median(matchR1m['mag_I'])))

# Galaxy2
print("\n# Galaxy2 catalog has %d objects" % len(useR2))

# Matched ones 
matchR2 = useR2[useR2['id'] > 0]
unMatchR2 = useR2[useR2['id'] == 0]
print("## %d synthetic galaxies are matched" % len(matchR2))
print("## %d synthetic galaxies are not matched" % len(unMatchR2))
# Median I-band magntiudes of the unmatched and matched ones
print("## Median magnitude for unmatched and matched galaxies: %5.3f, %5.3f" % (np.median(unMatchR2['mag_I']), 
                                                                                np.median(matchR2['mag_I'])))



    


















    






# HSC-R band
# Galaxy1 catalog has 54988 objects

## 53256 synthetic galaxies are matched
## 1732 synthetic galaxies are not matched
## Median magnitude for unmatched and matched galaxies: 24.560, 24.027

# Multipix ones!

## 53256 synthetic galaxies are matched
## 1732 synthetic galaxies are not matched
## Median magnitude for unmatched and matched galaxies: 24.560, 24.027

# Galaxy2 catalog has 54866 objects
## 52262 synthetic galaxies are matched
## 2604 synthetic galaxies are not matched
## Median magnitude for unmatched and matched galaxies: 24.616, 24.000

















In [16]:



    


# Galaxy1 
print("# HSC-I band")
print("# Galaxy1 catalog has %d objects" % len(useI1))

# Matched ones 
matchI1 = useI1[useI1['id'] > 0]
unMatchI1 = useI1[useI1['id'] == 0]
print("## %d synthetic galaxies are matched" % len(matchI1))
print("## %d synthetic galaxies are not matched" % len(unMatchI1))
# Median I-band magntiudes of the unmatched and matched ones
print("## Median magnitude for unmatched and matched galaxies: %5.3f, %5.3f" % (np.median(unMatchI1['mag_I']), 
                                                                                np.median(matchI1['mag_I'])))

# Matched ones 
print("\n# Multipix ones!")
matchI1m = useI1m[useI1m['id'] > 0]
unMatchI1m = useI1m[useI1m['id'] == 0]
print("## %d synthetic galaxies are matched" % len(matchI1m))
print("## %d synthetic galaxies are not matched" % len(unMatchI1m))
# Median I-band magntiudes of the unmatched and matched ones
print("## Median magnitude for unmatched and matched galaxies: %5.3f, %5.3f" % (np.median(unMatchI1m['mag_I']), 
                                                                                np.median(matchI1m['mag_I'])))

# Galaxy2
print("\n# Galaxy2 catalog has %d objects" % len(useI2))

# Matched ones 
matchI2 = useI2[useI2['id'] > 0]
unMatchI2 = useI2[useI2['id'] == 0]
print("## %d synthetic galaxies are matched" % len(matchI2))
print("## %d synthetic galaxies are not matched" % len(unMatchI2))
# Median I-band magntiudes of the unmatched and matched ones
print("## Median magnitude for unmatched and matched galaxies: %5.3f, %5.3f" % (np.median(unMatchI2['mag_I']), 
                                                                                np.median(matchI2['mag_I'])))



    


















    






# HSC-I band
# Galaxy1 catalog has 55031 objects
## 53312 synthetic galaxies are matched
## 1719 synthetic galaxies are not matched
## Median magnitude for unmatched and matched galaxies: 24.554, 24.027

# Multipix ones!
## 53313 synthetic galaxies are matched
## 1719 synthetic galaxies are not matched
## Median magnitude for unmatched and matched galaxies: 24.554, 24.027

# Galaxy2 catalog has 54979 objects
## 52481 synthetic galaxies are matched
## 2498 synthetic galaxies are not matched
## Median magnitude for unmatched and matched galaxies: 24.619, 24.001

















In [17]:



    


# Galaxy1 
print("# HSC-Z band")
print("# Galaxy1 catalog has %d objects" % len(useZ1))

# Matched ones 
matchZ1 = useZ1[useZ1['id'] > 0]
unMatchZ1 = useZ1[useZ1['id'] == 0]
print("\n## %d synthetic galaxies are matched" % len(matchZ1))
print("## %d synthetic galaxies are not matched" % len(unMatchZ1))
# Median I-band magntiudes of the unmatched and matched ones
print("## Median magnitude for unmatched and matched galaxies: %5.3f, %5.3f" % (np.median(unMatchZ1['mag_I']), 
                                                                                np.median(matchZ1['mag_I'])))

# Matched ones 
print("\n# Multipix ones!")
matchZ1m = useZ1m[useZ1m['id'] > 0]
unMatchZ1m = useZ1m[useZ1m['id'] == 0]
print("## %d synthetic galaxies are matched" % len(matchZ1m))
print("## %d synthetic galaxies are not matched" % len(unMatchZ1m))
# Median I-band magntiudes of the unmatched and matched ones
print("## Median magnitude for unmatched and matched galaxies: %5.3f, %5.3f" % (np.median(unMatchZ1m['mag_I']), 
                                                                                np.median(matchZ1m['mag_I'])))

# Galaxy2
print("\n# Galaxy2 catalog has %d objects" % len(useZ2))

# Matched ones 
matchZ2 = useZ2[useZ2['id'] > 0]
unMatchZ2 = useZ2[useZ2['id'] == 0]
print("## %d synthetic galaxies are matched" % len(matchZ2))
print("## %d synthetic galaxies are not matched" % len(unMatchZ2))
# Median I-band magntiudes of the unmatched and matched ones
print("## Median magnitude for unmatched and matched galaxies: %5.3f, %5.3f" % (np.median(unMatchZ2['mag_I']), 
                                                                                np.median(matchZ2['mag_I'])))



    


















    






# HSC-Z band
# Galaxy1 catalog has 55039 objects

## 53350 synthetic galaxies are matched
## 1689 synthetic galaxies are not matched
## Median magnitude for unmatched and matched galaxies: 24.568, 24.027

# Multipix ones!
## 53349 synthetic galaxies are matched
## 1689 synthetic galaxies are not matched
## Median magnitude for unmatched and matched galaxies: 24.568, 24.027

# Galaxy2 catalog has 55015 objects
## 52412 synthetic galaxies are matched
## 2603 synthetic galaxies are not matched
## Median magnitude for unmatched and matched galaxies: 24.650, 23.997

















In [18]:



    


# Galaxy1 
print("# HSC-Y band")
print("# Galaxy1 catalog has %d objects" % len(useY1))

# Matched ones 
matchY1 = useY1[useY1['id'] > 0]
unMatchY1 = useY1[useY1['id'] == 0]
print("\n## %d synthetic galaxies are matched" % len(matchY1))
print("## %d synthetic galaxies are not matched" % len(unMatchY1))
# Median I-band magntiudes of the unmatched and matched ones
print("## Median magnitude for unmatched and matched galaxies: %5.3f, %5.3f" % (np.median(unMatchY1['mag_I']), 
                                                                                np.median(matchY1['mag_I'])))

# Matched ones 
print("\n# Multipix ones!")
matchY1m = useY1m[useY1m['id'] > 0]
unMatchY1m = useY1m[useY1m['id'] == 0]
print("## %d synthetic galaxies are matched" % len(matchY1m))
print("## %d synthetic galaxies are not matched" % len(unMatchY1m))
# Median I-band magntiudes of the unmatched and matched ones
print("## Median magnitude for unmatched and matched galaxies: %5.3f, %5.3f" % (np.median(unMatchY1m['mag_I']), 
                                                                                np.median(matchY1m['mag_I'])))


# Galaxy2
print("\n# Galaxy2 catalog has %d objects" % len(useY2))

# Matched ones 
matchY2 = useY2[useY2['id'] > 0]
unMatchY2 = useY2[useY2['id'] == 0]
print("## %d synthetic galaxies are matched" % len(matchY2))
print("## %d synthetic galaxies are not matched" % len(unMatchY2))
# Median I-band magntiudes of the unmatched and matched ones
print("## Median magnitude for unmatched and matched galaxies: %5.3f, %5.3f" % (np.median(unMatchY2['mag_I']), 
                                                                                np.median(matchY2['mag_I'])))



    


















    






# HSC-Y band
# Galaxy1 catalog has 54793 objects

## 52459 synthetic galaxies are matched
## 2334 synthetic galaxies are not matched
## Median magnitude for unmatched and matched galaxies: 24.663, 24.007

# Multipix ones!
## 52459 synthetic galaxies are matched
## 2333 synthetic galaxies are not matched
## Median magnitude for unmatched and matched galaxies: 24.663, 24.007

# Galaxy2 catalog has 54724 objects
## 51086 synthetic galaxies are matched
## 3638 synthetic galaxies are not matched
## Median magnitude for unmatched and matched galaxies: 24.703, 23.967

Write the catalogs





In [38]:



    


matchG1.write(os.path.join(galaxyDir, 'gal_all1_HSC-G_match.fits'), format='fits', overwrite=True)
matchG2.write(os.path.join(galaxyDir, 'gal_all2_HSC-G_match.fits'), format='fits', overwrite=True)
unMatchG1.write(os.path.join(galaxyDir, 'gal_all1_HSC-G_unmatch.fits'), format='fits', overwrite=True)
unMatchG2.write(os.path.join(galaxyDir, 'gal_all2_HSC-G_unmatch.fits'), format='fits', overwrite=True)

matchR1.write(os.path.join(galaxyDir, 'gal_all1_HSC-R_match.fits'), format='fits', overwrite=True)
matchR2.write(os.path.join(galaxyDir, 'gal_all2_HSC-R_match.fits'), format='fits', overwrite=True)
unMatchR1.write(os.path.join(galaxyDir, 'gal_all1_HSC-R_unmatch.fits'), format='fits', overwrite=True)
unMatchR2.write(os.path.join(galaxyDir, 'gal_all2_HSC-R_unmatch.fits'), format='fits', overwrite=True)

matchI1.write(os.path.join(galaxyDir, 'gal_all1_HSC-I_match.fits'), format='fits', overwrite=True)
matchI2.write(os.path.join(galaxyDir, 'gal_all2_HSC-I_match.fits'), format='fits', overwrite=True)
unMatchI1.write(os.path.join(galaxyDir, 'gal_all1_HSC-I_unmatch.fits'), format='fits', overwrite=True)
unMatchI2.write(os.path.join(galaxyDir, 'gal_all2_HSC-I_unmatch.fits'), format='fits', overwrite=True)

matchZ1.write(os.path.join(galaxyDir, 'gal_all1_HSC-Z_match.fits'), format='fits', overwrite=True)
matchZ2.write(os.path.join(galaxyDir, 'gal_all2_HSC-Z_match.fits'), format='fits', overwrite=True)
unMatchZ1.write(os.path.join(galaxyDir, 'gal_all1_HSC-Z_unmatch.fits'), format='fits', overwrite=True)
unMatchZ2.write(os.path.join(galaxyDir, 'gal_all2_HSC-Z_unmatch.fits'), format='fits', overwrite=True)

matchY1.write(os.path.join(galaxyDir, 'gal_all1_HSC-Y_match.fits'), format='fits', overwrite=True)
matchY2.write(os.path.join(galaxyDir, 'gal_all2_HSC-Y_match.fits'), format='fits', overwrite=True)
unMatchY1.write(os.path.join(galaxyDir, 'gal_all1_HSC-Y_unmatch.fits'), format='fits', overwrite=True)
unMatchY2.write(os.path.join(galaxyDir, 'gal_all2_HSC-Y_unmatch.fits'), format='fits', overwrite=True)



    











In [39]:



    


matchG1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-G_match_multipix.fits'), format='fits', overwrite=True)
unMatchG1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-G_unmatch_multipix.fits'), format='fits', overwrite=True)

matchR1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-R_match_multipix.fits'), format='fits', overwrite=True)
unMatchR1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-R_unmatch_multipix.fits'), format='fits', overwrite=True)

matchI1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-I_match_multipix.fits'), format='fits', overwrite=True)
unMatchI1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-I_unmatch_multipix.fits'), format='fits', overwrite=True)

matchZ1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-Z_match_multipix.fits'), format='fits', overwrite=True)
unMatchZ1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-Z_unmatch_multipix.fits'), format='fits', overwrite=True)

matchY1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-Y_match_multipix.fits'), format='fits', overwrite=True)
unMatchY1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-Y_unmatch_multipix.fits'), format='fits', overwrite=True)

Read the catalogs





In [7]:



    


matchG1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-G_match.fits'), format='fits')
matchG2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-G_match.fits'), format='fits')
unMatchG1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-G_unmatch.fits'), format='fits')
unMatchG2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-G_unmatch.fits'), format='fits')

matchR1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-R_match.fits'), format='fits')
matchR2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-R_match.fits'), format='fits')
unMatchR1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-R_unmatch.fits'), format='fits')
unMatchR2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-R_unmatch.fits'), format='fits')

matchI1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-I_match.fits'), format='fits')
matchI2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-I_match.fits'), format='fits')
unMatchI1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-I_unmatch.fits'), format='fits')
unMatchI2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-I_unmatch.fits'), format='fits')

matchZ1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-Z_match.fits'), format='fits')
matchZ2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-Z_match.fits'), format='fits')
unMatchZ1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-Z_unmatch.fits'), format='fits')
unMatchZ2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-Z_unmatch.fits'), format='fits')

matchY1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-Y_match.fits'), format='fits')
matchY2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-Y_match.fits'), format='fits')
unMatchY1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-Y_unmatch.fits'), format='fits')
unMatchY2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-Y_unmatch.fits'), format='fits')



    











In [8]:



    


matchG1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-G_match_multipix.fits'), format='fits')
unMatchG1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-G_unmatch_multipix.fits'), format='fits')

matchR1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-R_match_multipix.fits'), format='fits')
unMatchR1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-R_unmatch_multipix.fits'), format='fits')

matchI1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-I_match_multipix.fits'), format='fits')
unMatchI1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-I_unmatch_multipix.fits'), format='fits')

matchZ1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-Z_match_multipix.fits'), format='fits')
unMatchZ1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-Z_unmatch_multipix.fits'), format='fits')

matchY1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-Y_match_multipix.fits'), format='fits')
unMatchY1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-Y_unmatch_multipix.fits'), format='fits')

Select the primary and non-parent detections

Generate new catalogs





In [9]:



    


print([col for col in catG1.colnames if 'primary' in col.lower()])
print([col for col in catG1.colnames if 'match' in col.lower()])
print([col for col in catG1.colnames if 'child' in col.lower()])
print([col for col in catG1.colnames if 'close' in col.lower()])
print([col for col in catG1.colnames if 'extend' in col.lower()])



    


















    






['detect.is-primary', 'nPrimary']
['nMatched', 'rMatched']
['blendedness.abs.flux.child', 'blendedness.abs.shape.child_a', 'blendedness.abs.shape.child_q', 'blendedness.abs.shape.child_theta', 'blendedness.raw.flux.child', 'blendedness.raw.shape.child_a', 'blendedness.raw.shape.child_q', 'blendedness.raw.shape.child_theta', 'deblend.nchild', 'force.deblend.nchild', 'nNoChild']
['fakeClosest']
['classification.extendedness', 'force.classification.extendedness']

















In [8]:



    


print("# HSC-G band")
primaryG1 = matchG1[(matchG1['detect.is-primary']) &
                    (matchG1['deblend.nchild'] == 0) & 
                    (matchG1['fakeClosest'])] 
print("#  1: %d galaxies are primary, non-parent, closest matches" % len(primaryG1))

primaryG1m = matchG1m[(matchG1m['detect.is-primary']) &
                      (matchG1m['deblend.nchild'] == 0) & 
                      (matchG1m['fakeClosest'])] 
print("# 1m: %d galaxies are primary, non-parent, closest matches" % len(primaryG1m))

primaryG2 = matchG2[(matchG2['detect.is-primary']) &
                    (matchG2['deblend.nchild'] == 0) & 
                    (matchG2['fakeClosest'])] 
print("#  2: %d galaxies are primary, non-parent, closest matches" % len(primaryG2))



    


















    






# HSC-G band
#  1: 31460 galaxies are primary, non-parent, closest matches
# 1m: 31461 galaxies are primary, non-parent, closest matches
#  2: 31347 galaxies are primary, non-parent, closest matches

















In [9]:



    


print("# HSC-R band")
primaryR1 = matchR1[(matchR1['detect.is-primary']) &
                    (matchR1['deblend.nchild'] == 0) & 
                    (matchR1['fakeClosest'])] 
print("#  1: %d galaxies are primary, non-parent, closest matches" % len(primaryR1))

primaryR1m = matchR1m[(matchR1m['detect.is-primary']) &
                      (matchR1m['deblend.nchild'] == 0) & 
                      (matchR1m['fakeClosest'])] 
print("# 1m: %d galaxies are primary, non-parent, closest matches" % len(primaryR1m))

primaryR2 = matchR2[(matchR2['detect.is-primary']) &
                    (matchR2['deblend.nchild'] == 0) & 
                    (matchR2['fakeClosest'])] 
print("#  2: %d galaxies are primary, non-parent, closest matches" % len(primaryR2))



    


















    






# HSC-R band
#  1: 31543 galaxies are primary, non-parent, closest matches
# 1m: 31543 galaxies are primary, non-parent, closest matches
#  2: 31646 galaxies are primary, non-parent, closest matches

















In [10]:



    


print("# HSC-I band")
primaryI1 = matchI1[(matchI1['detect.is-primary']) &
                    (matchI1['deblend.nchild'] == 0) & 
                    (matchI1['fakeClosest'])] 
print("#  1: %d galaxies are primary, non-parent, closest matches" % len(primaryI1))

primaryI1m = matchI1m[(matchI1m['detect.is-primary']) &
                      (matchI1m['deblend.nchild'] == 0) & 
                      (matchI1m['fakeClosest'])] 
print("# 1m: %d galaxies are primary, non-parent, closest matches" % len(primaryI1m))

primaryI2 = matchI2[(matchI2['detect.is-primary']) &
                    (matchI2['deblend.nchild'] == 0) & 
                    (matchI2['fakeClosest'])] 
print("#  2: %d galaxies are primary, non-parent, closest matches" % len(primaryI2))



    


















    






# HSC-I band
#  1: 31130 galaxies are primary, non-parent, closest matches
# 1m: 31130 galaxies are primary, non-parent, closest matches
#  2: 31562 galaxies are primary, non-parent, closest matches

















In [11]:



    


print("# HSC-Z band")
primaryZ1 = matchZ1[(matchZ1['detect.is-primary']) &
                    (matchZ1['deblend.nchild'] == 0) & 
                    (matchZ1['fakeClosest'])] 
print("#  1: %d galaxies are primary, non-parent, closest matches" % len(primaryZ1))

primaryZ1m = matchZ1m[(matchZ1m['detect.is-primary']) &
                      (matchZ1m['deblend.nchild'] == 0) & 
                      (matchZ1m['fakeClosest'])] 
print("# 1m: %d galaxies are primary, non-parent, closest matches" % len(primaryZ1m))

primaryZ2 = matchZ2[(matchZ2['detect.is-primary']) &
                    (matchZ2['deblend.nchild'] == 0) & 
                    (matchZ2['fakeClosest'])] 
print("#  2: %d galaxies are primary, non-parent, closest matches" % len(primaryZ2))



    


















    






# HSC-Z band
#  1: 31609 galaxies are primary, non-parent, closest matches
# 1m: 31609 galaxies are primary, non-parent, closest matches
#  2: 31421 galaxies are primary, non-parent, closest matches

















In [12]:



    


print("# HSC-Y band")
primaryY1 = matchY1[(matchY1['detect.is-primary']) &
                    (matchY1['deblend.nchild'] == 0) & 
                    (matchY1['fakeClosest'])] 
print("#  1: %d galaxies are primary, non-parent, closest matches" % len(primaryY1))

primaryY1m = matchY1m[(matchY1m['detect.is-primary']) &
                      (matchY1m['deblend.nchild'] == 0) & 
                      (matchY1m['fakeClosest'])] 
print("# 1m: %d galaxies are primary, non-parent, closest matches" % len(primaryY1m))

primaryY2 = matchY2[(matchY2['detect.is-primary']) &
                    (matchY2['deblend.nchild'] == 0) & 
                    (matchY2['fakeClosest'])] 
print("#  2: %d galaxies are primary, non-parent, closest matches" % len(primaryY2))



    


















    






# HSC-Y band
#  1: 30631 galaxies are primary, non-parent, closest matches
# 1m: 30631 galaxies are primary, non-parent, closest matches
#  2: 30228 galaxies are primary, non-parent, closest matches

Write the catalogs





In [47]:



    


primaryG1.write(os.path.join(galaxyDir, 'gal_all1_HSC-G_primary.fits'), format='fits', overwrite=True)
primaryG2.write(os.path.join(galaxyDir, 'gal_all2_HSC-G_primary.fits'), format='fits', overwrite=True)

primaryR1.write(os.path.join(galaxyDir, 'gal_all1_HSC-R_primary.fits'), format='fits', overwrite=True)
primaryR2.write(os.path.join(galaxyDir, 'gal_all2_HSC-R_primary.fits'), format='fits', overwrite=True)

primaryI1.write(os.path.join(galaxyDir, 'gal_all1_HSC-I_primary.fits'), format='fits', overwrite=True)
primaryI2.write(os.path.join(galaxyDir, 'gal_all2_HSC-I_primary.fits'), format='fits', overwrite=True)

primaryZ1.write(os.path.join(galaxyDir, 'gal_all1_HSC-Z_primary.fits'), format='fits', overwrite=True)
primaryZ2.write(os.path.join(galaxyDir, 'gal_all2_HSC-Z_primary.fits'), format='fits', overwrite=True)

primaryY1.write(os.path.join(galaxyDir, 'gal_all1_HSC-Y_primary.fits'), format='fits', overwrite=True)
primaryY2.write(os.path.join(galaxyDir, 'gal_all2_HSC-Y_primary.fits'), format='fits', overwrite=True)

primaryG1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-G_primary_multipix.fits'), format='fits', overwrite=True)
primaryR1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-R_primary_multipix.fits'), format='fits', overwrite=True)
primaryI1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-I_primary_multipix.fits'), format='fits', overwrite=True)
primaryZ1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-Z_primary_multipix.fits'), format='fits', overwrite=True)
primaryY1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-Y_primary_multipix.fits'), format='fits', overwrite=True)

Read the new catalogs





In [10]:



    


primaryG1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-G_primary.fits'), format='fits')
primaryG2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-G_primary.fits'), format='fits')

primaryR1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-R_primary.fits'), format='fits')
primaryR2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-R_primary.fits'), format='fits')

primaryI1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-I_primary.fits'), format='fits')
primaryI2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-I_primary.fits'), format='fits')

primaryZ1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-Z_primary.fits'), format='fits')
primaryZ2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-Z_primary.fits'), format='fits')

primaryY1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-Y_primary.fits'), format='fits')
primaryY2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-Y_primary.fits'), format='fits')

Quality control

Generate new catalogs





In [11]:



    


print([col for col in catG1.colnames if 'flag' in col.lower()])
print("\n")
print([col for col in catG1.colnames if 'cmodel' in col.lower()])
print("\n")
print([col for col in catG1.colnames if 'kron' in col.lower()])



    


















    






['blendedness.flags', 'blendedness.flags.noCentroid', 'blendedness.flags.noShape', 'centroid.naive.flags', 'centroid.sdss.flags', 'cmodel.dev.flags.maxIter', 'cmodel.dev.flags.numericError', 'cmodel.dev.flags.trSmall', 'cmodel.dev.flux.flags', 'cmodel.dev.flux.flags.apcorr', 'cmodel.exp.flags.maxIter', 'cmodel.exp.flags.numericError', 'cmodel.exp.flags.trSmall', 'cmodel.exp.flux.flags', 'cmodel.exp.flux.flags.apcorr', 'cmodel.flags.badCentroid', 'cmodel.flags.noCalib', 'cmodel.flags.noPsf', 'cmodel.flags.noShape', 'cmodel.flags.noWcs', 'cmodel.flags.region.maxArea', 'cmodel.flags.region.maxBadPixelFraction', 'cmodel.flags.region.usedFootprintArea', 'cmodel.flags.region.usedInitialEllipseMax', 'cmodel.flags.region.usedInitialEllipseMin', 'cmodel.flags.region.usedPsfArea', 'cmodel.flags.smallShape', 'cmodel.flux.flags', 'cmodel.flux.flags.apcorr', 'cmodel.initial.flags.maxIter', 'cmodel.initial.flags.numericError', 'cmodel.initial.flags.trSmall', 'cmodel.initial.flux.flags', 'flags.badcentroid', 'flags.negative', 'flags.pixel.bad', 'flags.pixel.bright.object.any', 'flags.pixel.bright.object.center', 'flags.pixel.clipped.any', 'flags.pixel.cr.any', 'flags.pixel.cr.center', 'flags.pixel.edge', 'flags.pixel.interpolated.any', 'flags.pixel.interpolated.center', 'flags.pixel.offimage', 'flags.pixel.saturated.any', 'flags.pixel.saturated.center', 'flags.pixel.suspect.any', 'flags.pixel.suspect.center', 'flux.aperture.flags', 'flux.gaussian.flags', 'flux.gaussian.flags.apcorr', 'flux.kron.flags', 'flux.kron.flags.apcorr', 'flux.kron.flags.badShape', 'flux.kron.flags.edge', 'flux.kron.flags.radius', 'flux.kron.flags.smallRadius', 'flux.kron.flags.usedMinimumRadius', 'flux.kron.flags.usedPsfRadius', 'flux.naive.flags', 'flux.psf.flags', 'flux.psf.flags.apcorr', 'flux.sinc.flags', 'multishapelet.psf.flags', 'multishapelet.psf.flags.constraint.q', 'multishapelet.psf.flags.constraint.r', 'multishapelet.psf.flags.maxiter', 'multishapelet.psf.flags.tinystep', 'shape.hsm.moments.centroid.flags', 'shape.hsm.moments.flags', 'shape.hsm.psfMoments.centroid.flags', 'shape.hsm.psfMoments.flags', 'shape.hsm.regauss.flags', 'shape.sdss.centroid.flags', 'shape.sdss.flags', 'shape.sdss.flags.maxiter', 'shape.sdss.flags.psf', 'shape.sdss.flags.shift', 'shape.sdss.flags.unweighted', 'shape.sdss.flags.unweightedbad']


['cmodel.center', 'cmodel.dev.ellipse_a', 'cmodel.dev.ellipse_q', 'cmodel.dev.ellipse_theta', 'cmodel.dev.fixed', 'cmodel.dev.flags.maxIter', 'cmodel.dev.flags.numericError', 'cmodel.dev.flags.trSmall', 'cmodel.dev.flux', 'cmodel.dev.flux.apcorr', 'cmodel.dev.flux.apcorr.err', 'cmodel.dev.flux.err', 'cmodel.dev.flux.flags', 'cmodel.dev.flux.flags.apcorr', 'cmodel.dev.flux.inner', 'cmodel.dev.nIter', 'cmodel.dev.nonlinear', 'cmodel.dev.objective', 'cmodel.dev.time', 'cmodel.ellipse_a', 'cmodel.ellipse_q', 'cmodel.ellipse_theta', 'cmodel.exp.ellipse_a', 'cmodel.exp.ellipse_q', 'cmodel.exp.ellipse_theta', 'cmodel.exp.fixed', 'cmodel.exp.flags.maxIter', 'cmodel.exp.flags.numericError', 'cmodel.exp.flags.trSmall', 'cmodel.exp.flux', 'cmodel.exp.flux.apcorr', 'cmodel.exp.flux.apcorr.err', 'cmodel.exp.flux.err', 'cmodel.exp.flux.flags', 'cmodel.exp.flux.flags.apcorr', 'cmodel.exp.flux.inner', 'cmodel.exp.nIter', 'cmodel.exp.nonlinear', 'cmodel.exp.objective', 'cmodel.exp.time', 'cmodel.flags.badCentroid', 'cmodel.flags.noCalib', 'cmodel.flags.noPsf', 'cmodel.flags.noShape', 'cmodel.flags.noWcs', 'cmodel.flags.region.maxArea', 'cmodel.flags.region.maxBadPixelFraction', 'cmodel.flags.region.usedFootprintArea', 'cmodel.flags.region.usedInitialEllipseMax', 'cmodel.flags.region.usedInitialEllipseMin', 'cmodel.flags.region.usedPsfArea', 'cmodel.flags.smallShape', 'cmodel.flux', 'cmodel.flux.apcorr', 'cmodel.flux.apcorr.err', 'cmodel.flux.err', 'cmodel.flux.flags', 'cmodel.flux.flags.apcorr', 'cmodel.flux.inner', 'cmodel.fracDev', 'cmodel.initial.ellipse_a', 'cmodel.initial.ellipse_q', 'cmodel.initial.ellipse_theta', 'cmodel.initial.fixed', 'cmodel.initial.flags.maxIter', 'cmodel.initial.flags.numericError', 'cmodel.initial.flags.trSmall', 'cmodel.initial.flux', 'cmodel.initial.flux.err', 'cmodel.initial.flux.flags', 'cmodel.initial.flux.inner', 'cmodel.initial.nIter', 'cmodel.initial.nonlinear', 'cmodel.initial.objective', 'cmodel.initial.time', 'cmodel.objective', 'cmodel.region.final.ellipse_a', 'cmodel.region.final.ellipse_q', 'cmodel.region.final.ellipse_theta', 'cmodel.region.initial.ellipse_a', 'cmodel.region.initial.ellipse_q', 'cmodel.region.initial.ellipse_theta', 'force.cmodel.dev.flux', 'force.cmodel.dev.flux.apcorr', 'force.cmodel.dev.flux.apcorr.err', 'force.cmodel.dev.flux.err', 'force.cmodel.exp.flux', 'force.cmodel.exp.flux.apcorr', 'force.cmodel.exp.flux.apcorr.err', 'force.cmodel.exp.flux.err', 'force.cmodel.flux', 'force.cmodel.flux.apcorr', 'force.cmodel.flux.apcorr.err', 'force.cmodel.flux.err', 'force.cmodel.fracDev', 'cmodel.dev.mag', 'cmodel.dev.mag.err', 'cmodel.dev.mag.apcorr', 'cmodel.dev.mag.apcorr.err', 'cmodel.exp.mag', 'cmodel.exp.mag.err', 'cmodel.exp.mag.apcorr', 'cmodel.exp.mag.apcorr.err', 'cmodel.mag', 'cmodel.mag.err', 'cmodel.mag.apcorr', 'cmodel.mag.apcorr.err', 'cmodel.initial.mag', 'cmodel.initial.mag.err', 'force.cmodel.dev.mag', 'force.cmodel.dev.mag.err', 'force.cmodel.dev.mag.apcorr', 'force.cmodel.dev.mag.apcorr.err', 'force.cmodel.exp.mag', 'force.cmodel.exp.mag.err', 'force.cmodel.exp.mag.apcorr', 'force.cmodel.exp.mag.apcorr.err', 'force.cmodel.mag', 'force.cmodel.mag.err', 'force.cmodel.mag.apcorr', 'force.cmodel.mag.apcorr.err']


['flux.kron', 'flux.kron.apcorr', 'flux.kron.apcorr.err', 'flux.kron.err', 'flux.kron.flags', 'flux.kron.flags.apcorr', 'flux.kron.flags.badShape', 'flux.kron.flags.edge', 'flux.kron.flags.radius', 'flux.kron.flags.smallRadius', 'flux.kron.flags.usedMinimumRadius', 'flux.kron.flags.usedPsfRadius', 'flux.kron.psfRadius', 'flux.kron.radius', 'flux.kron.radiusForRadius', 'force.flux.kron', 'force.flux.kron.apcorr', 'force.flux.kron.apcorr.err', 'force.flux.kron.err', 'mag.kron', 'mag.kron.err', 'mag.kron.apcorr', 'mag.kron.apcorr.err', 'force.mag.kron', 'force.mag.kron.err', 'force.mag.kron.apcorr', 'force.mag.kron.apcorr.err']

















In [8]:



    


def galQC(cat, kron=False, blendedness=False):
    """Quanlity control cuts."""
    use = cat[(~cat['centroid.sdss.flags']) & 
              (~cat['flags.badcentroid']) & 
              (~cat['flags.pixel.bad']) &
              (~cat['flags.pixel.cr.center']) &
              (~cat['flags.pixel.edge']) & 
              (~cat['flags.pixel.interpolated.center']) &
              (~cat['flags.pixel.offimage']) &
              (~cat['flags.pixel.saturated.center']) & 
              (~cat['flags.pixel.suspect.center']) &
              (~cat['cmodel.flux.flags']) &
              (np.isfinite(cat['force.classification.extendedness'])) &
              (np.isfinite(cat['force.cmodel.mag'])) &
              (np.isfinite(cat['force.cmodel.mag.err']))]
    
    if kron:
        use = use[(np.isfinite(cat['force.mag.kron'])) &
                  (np.isfinite(cat['force.mag.kron.err']))]
    
    if blendedness:
        use = use[~cat['blendedness.flags']]
    
    return use



    











In [16]:



    


print("# HSC-G band")
galUseG1 = galQC(primaryG1)
print("##  1: %d / %d galaxies are useful" % (len(galUseG1), len(primaryG1)))
galUseG1m = galQC(primaryG1m)
print("## 1m: %d / %d galaxies are useful" % (len(galUseG1m), len(primaryG1m)))
galUseG2 = galQC(primaryG2)
print("##  2: %d / %d galaxies are useful" % (len(galUseG2), len(primaryG2)))

print("\n# HSC-R band")
galUseR1 = galQC(primaryR1)
print("##  1: %d / %d galaxies are useful" % (len(galUseR1), len(primaryR1)))
galUseR1m = galQC(primaryR1m)
print("## 1m: %d / %d galaxies are useful" % (len(galUseR1m), len(primaryR1m)))
galUseR2 = galQC(primaryR2)
print("##  2: %d / %d galaxies are useful" % (len(galUseR2), len(primaryR2)))

print("\n# HSC-I band")
galUseI1 = galQC(primaryI1)
print("##  1: %d / %d galaxies are useful" % (len(galUseI1), len(primaryI1)))
galUseI1m = galQC(primaryI1m)
print("## 1m: %d / %d galaxies are useful" % (len(galUseI1m), len(primaryI1m)))
galUseI2 = galQC(primaryI2)
print("##  2: %d / %d galaxies are useful" % (len(galUseI2), len(primaryI2)))

print("\n# HSC-Z band")
galUseZ1 = galQC(primaryZ1)
print("##  1: %d / %d galaxies are useful" % (len(galUseZ1), len(primaryZ1)))
galUseZ1m = galQC(primaryZ1m)
print("## 1m: %d / %d galaxies are useful" % (len(galUseZ1m), len(primaryZ1m)))
galUseZ2 = galQC(primaryZ2)
print("##  2: %d / %d galaxies are useful" % (len(galUseZ2), len(primaryZ2)))

print("\n# HSC-Y band")
galUseY1 = galQC(primaryY1)
print("##  1: %d / %d galaxies are useful" % (len(galUseY1), len(primaryY1)))
galUseY1m = galQC(primaryY1m)
print("## 1m: %d / %d galaxies are useful" % (len(galUseY1m), len(primaryY1m)))
galUseY2 = galQC(primaryY2)
print("##  2: %d / %d galaxies are useful" % (len(galUseY2), len(primaryY2)))



    


















    






# HSC-G band
##  1: 30128 / 31460 galaxies are useful
## 1m: 30089 / 31461 galaxies are useful
##  2: 30186 / 31347 galaxies are useful

# HSC-R band
##  1: 30578 / 31543 galaxies are useful
## 1m: 30569 / 31543 galaxies are useful
##  2: 30873 / 31646 galaxies are useful

# HSC-I band
##  1: 30331 / 31130 galaxies are useful
## 1m: 30329 / 31130 galaxies are useful
##  2: 31119 / 31562 galaxies are useful

# HSC-Z band
##  1: 30597 / 31609 galaxies are useful
## 1m: 30587 / 31609 galaxies are useful
##  2: 30336 / 31421 galaxies are useful

# HSC-Y band
##  1: 26772 / 30631 galaxies are useful
## 1m: 26747 / 30631 galaxies are useful
##  2: 26977 / 30228 galaxies are useful

Write the catalogs





In [50]:



    


galUseG1.write(os.path.join(galaxyDir, 'gal_all1_HSC-G_qc.fits'), format='fits', overwrite=True)
galUseG2.write(os.path.join(galaxyDir, 'gal_all2_HSC-G_qc.fits'), format='fits', overwrite=True)

galUseR1.write(os.path.join(galaxyDir, 'gal_all1_HSC-R_qc.fits'), format='fits', overwrite=True)
galUseR2.write(os.path.join(galaxyDir, 'gal_all2_HSC-R_qc.fits'), format='fits', overwrite=True)

galUseI1.write(os.path.join(galaxyDir, 'gal_all1_HSC-I_qc.fits'), format='fits', overwrite=True)
galUseI2.write(os.path.join(galaxyDir, 'gal_all2_HSC-I_qc.fits'), format='fits', overwrite=True)

galUseZ1.write(os.path.join(galaxyDir, 'gal_all1_HSC-Z_qc.fits'), format='fits', overwrite=True)
galUseZ2.write(os.path.join(galaxyDir, 'gal_all2_HSC-Z_qc.fits'), format='fits', overwrite=True)

galUseY1.write(os.path.join(galaxyDir, 'gal_all1_HSC-Y_qc.fits'), format='fits', overwrite=True)
galUseY2.write(os.path.join(galaxyDir, 'gal_all2_HSC-Y_qc.fits'), format='fits', overwrite=True)

galUseG1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-G_qc_multipix.fits'), format='fits', overwrite=True)
galUseR1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-R_qc_multipix.fits'), format='fits', overwrite=True)
galUseI1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-I_qc_multipix.fits'), format='fits', overwrite=True)
galUseZ1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-Z_qc_multipix.fits'), format='fits', overwrite=True)
galUseY1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-Y_qc_multipix.fits'), format='fits', overwrite=True)

Read the catalogs





In [12]:



    


galUseG1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-G_qc.fits'), format='fits')
galUseG2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-G_qc.fits'), format='fits')

galUseR1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-R_qc.fits'), format='fits')
galUseR2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-R_qc.fits'), format='fits')

galUseI1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-I_qc.fits'), format='fits')
galUseI2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-I_qc.fits'), format='fits')

galUseZ1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-Z_qc.fits'), format='fits')
galUseZ2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-Z_qc.fits'), format='fits')

galUseY1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-Y_qc.fits'), format='fits')
galUseY2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-Y_qc.fits'), format='fits')

galUseG1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-G_qc_multipix.fits'), format='fits')
galUseR1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-R_qc_multipix.fits'), format='fits')
galUseI1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-I_qc_multipix.fits'), format='fits')
galUseZ1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-Z_qc_multipix.fits'), format='fits')
galUseY1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-Y_qc_multipix.fits'), format='fits')

Isolate the galaxies that got misclassified

Generate new catalogs





In [18]:



    


print("# HSC-G band")
misClassG1 = galUseG1[(galUseG1['force.classification.extendedness'] <= 0.1)]
galGoodG1 = galUseG1[(galUseG1['force.classification.extendedness'] > 0.1)]
print("#  1: %d galaxies are mis-classified as extended objects" % len(misClassG1))
print("## median magnitude for misclassified and good galaxies: %5.3f, %5.3f" % (np.nanmedian(misClassG1['mag_G']), 
                                                                                 np.nanmedian(galGoodG1['mag_G'])))
print("## median blendedness for misclassified and good galaxies: %6.4f, %6.4f" % (np.nanmedian(misClassG1['blendedness.abs.flux']), 
                                                                                   np.nanmedian(galGoodG1['blendedness.abs.flux'])))

misClassG1m = galUseG1m[(galUseG1m['force.classification.extendedness'] <= 0.5)]
galGoodG1m = galUseG1m[(galUseG1m['force.classification.extendedness'] > 0.5)]
print("# 1m: %d galaxies are mis-classified as extended objects" % len(misClassG1m))
print("## median magnitude for misclassified and good galaxies: %5.3f, %5.3f" % (np.nanmedian(misClassG1m['mag_G']), 
                                                                                 np.nanmedian(galGoodG1m['mag_G'])))
print("## median blendedness for misclassified and good galaxies: %6.4f, %6.4f" % (np.nanmedian(misClassG1m['blendedness.abs.flux']), 
                                                                                   np.nanmedian(galGoodG1m['blendedness.abs.flux'])))

misClassG2 = galUseG2[(galUseG2['force.classification.extendedness'] <= 0.5)]
galGoodG2 = galUseG2[(galUseG2['force.classification.extendedness'] > 0.5)]
print("#  2: %d galaxies are mis-classified as extended objects" % len(misClassG2))
print("## median magnitude for misclassified and good galaxies: %5.3f, %5.3f" % (np.nanmedian(misClassG2['mag_G']), 
                                                                                 np.nanmedian(galGoodG2['mag_G'])))
print("## median blendedness for misclassified and good galaxies: %6.4f, %6.4f" % (np.nanmedian(misClassG2['blendedness.abs.flux']), 
                                                                                   np.nanmedian(galGoodG2['blendedness.abs.flux'])))



    


















    






# HSC-G band
#  1: 1034 galaxies are mis-classified as extended objects
## median magnitude for misclassified and good galaxies: 25.698, 25.063
## median blendedness for misclassified and good galaxies: 0.0069, 0.0015
# 1m: 978 galaxies are mis-classified as extended objects
## median magnitude for misclassified and good galaxies: 25.710, 25.064
## median blendedness for misclassified and good galaxies: 0.0065, 0.0014
#  2: 2021 galaxies are mis-classified as extended objects
## median magnitude for misclassified and good galaxies: 25.523, 25.045
## median blendedness for misclassified and good galaxies: 0.0002, 0.0005

















In [19]:



    


print("# HSC-R band")
misClassR1 = galUseR1[(galUseR1['force.classification.extendedness'] <= 0.5)]
galGoodR1 = galUseR1[(galUseR1['force.classification.extendedness'] > 0.5)]
print("#  1: %d galaxies are mis-classified as extended objects" % len(misClassR1))
print("## median magnitude for misclassified and good galaxies: %5.3f, %5.3f" % (np.nanmedian(misClassR1['mag_R']), 
                                                                                 np.nanmedian(galGoodR1['mag_R'])))
print("## median blendedness for misclassified and good galaxies: %6.4f, %6.4f" % (np.nanmedian(misClassR1['blendedness.abs.flux']), 
                                                                                   np.nanmedian(galGoodR1['blendedness.abs.flux'])))

misClassR1m = galUseR1m[(galUseR1m['force.classification.extendedness'] <= 0.5)]
galGoodR1m = galUseR1m[(galUseR1m['force.classification.extendedness'] > 0.5)]
print("# 1m: %d galaxies are mis-classified as extended objects" % len(misClassR1m))
print("## median magnitude for misclassified and good galaxies: %5.3f, %5.3f" % (np.nanmedian(misClassR1m['mag_R']), 
                                                                                 np.nanmedian(galGoodR1m['mag_R'])))
print("## median blendedness for misclassified and good galaxies: %6.4f, %6.4f" % (np.nanmedian(misClassR1m['blendedness.abs.flux']), 
                                                                                   np.nanmedian(galGoodR1m['blendedness.abs.flux'])))

misClassR2 = galUseR2[(galUseR2['force.classification.extendedness'] <= 0.5)]
galGoodR2 = galUseR2[(galUseR2['force.classification.extendedness'] > 0.5)]
print("#  2: %d galaxies are mis-classified as extended objects" % len(misClassR2))
print("## median magnitude for misclassified and good galaxies: %5.3f, %5.3f" % (np.nanmedian(misClassR2['mag_R']), 
                                                                                 np.nanmedian(galGoodR2['mag_R'])))
print("## median blendedness for misclassified and good galaxies: %6.4f, %6.4f" % (np.nanmedian(misClassR2['blendedness.abs.flux']), 
                                                                                   np.nanmedian(galGoodR2['blendedness.abs.flux'])))



    


















    






# HSC-R band
#  1: 368 galaxies are mis-classified as extended objects
## median magnitude for misclassified and good galaxies: 25.232, 24.644
## median blendedness for misclassified and good galaxies: 0.0028, 0.0001
# 1m: 298 galaxies are mis-classified as extended objects
## median magnitude for misclassified and good galaxies: 25.238, 24.646
## median blendedness for misclassified and good galaxies: 0.0058, 0.0001
#  2: 1908 galaxies are mis-classified as extended objects
## median magnitude for misclassified and good galaxies: 25.176, 24.589
## median blendedness for misclassified and good galaxies: 0.0000, 0.0002

















In [20]:



    


print("# HSC-I band")
misClassI1 = galUseI1[(galUseI1['force.classification.extendedness'] <= 0.5)]
galGoodI1 = galUseI1[(galUseI1['force.classification.extendedness'] > 0.5)]
print("#  1: %d galaxies are mis-classified as extended objects" % len(misClassI1))
print("## median magnitude for misclassified and good galaxies: %5.3f, %5.3f" % (np.nanmedian(misClassI1['mag_I']), 
                                                                                 np.nanmedian(galGoodI1['mag_I'])))
print("## median blendedness for misclassified and good galaxies: %6.4f, %6.4f" % (np.nanmedian(misClassI1['blendedness.abs.flux']), 
                                                                                   np.nanmedian(galGoodI1['blendedness.abs.flux'])))

misClassI1m = galUseI1m[(galUseI1m['force.classification.extendedness'] <= 0.5)]
galGoodI1m = galUseI1m[(galUseI1m['force.classification.extendedness'] > 0.5)]
print("# 1m: %d galaxies are mis-classified as extended objects" % len(misClassI1m))
print("## median magnitude for misclassified and good galaxies: %5.3f, %5.3f" % (np.nanmedian(misClassI1m['mag_I']), 
                                                                                 np.nanmedian(galGoodI1m['mag_I'])))
print("## median blendedness for misclassified and good galaxies: %6.4f, %6.4f" % (np.nanmedian(misClassI1m['blendedness.abs.flux']), 
                                                                                   np.nanmedian(galGoodI1m['blendedness.abs.flux'])))

misClassI2 = galUseI2[(galUseI2['force.classification.extendedness'] <= 0.5)]
galGoodI2 = galUseI2[(galUseI2['force.classification.extendedness'] > 0.5)]
print("#  2: %d galaxies are mis-classified as extended objects" % len(misClassI2))
print("## median magnitude for misclassified and good galaxies: %5.3f, %5.3f" % (np.nanmedian(misClassI2['mag_I']), 
                                                                                 np.nanmedian(galGoodI2['mag_I'])))
print("## median blendedness for misclassified and good galaxies: %6.4f, %6.4f" % (np.nanmedian(misClassI2['blendedness.abs.flux']), 
                                                                                   np.nanmedian(galGoodI2['blendedness.abs.flux'])))



    


















    






# HSC-I band
#  1: 337 galaxies are mis-classified as extended objects
## median magnitude for misclassified and good galaxies: 24.942, 24.213
## median blendedness for misclassified and good galaxies: 0.0013, 0.0000
# 1m: 270 galaxies are mis-classified as extended objects
## median magnitude for misclassified and good galaxies: 24.962, 24.216
## median blendedness for misclassified and good galaxies: 0.0043, 0.0000
#  2: 2037 galaxies are mis-classified as extended objects
## median magnitude for misclassified and good galaxies: 24.848, 24.159
## median blendedness for misclassified and good galaxies: 0.0000, 0.0002

















In [21]:



    


print("# HSC-Z band")
misClassZ1 = galUseZ1[(galUseZ1['force.classification.extendedness'] <= 0.5)]
galGoodZ1 = galUseZ1[(galUseZ1['force.classification.extendedness'] > 0.5)]
print("#  1: %d galaxies are mis-classified as extended objects" % len(misClassZ1))
print("## median magnitude for misclassified and good galaxies: %5.3f, %5.3f" % (np.nanmedian(misClassZ1['mag_Z']), 
                                                                                 np.nanmedian(galGoodZ1['mag_Z'])))
print("## median blendedness for misclassified and good galaxies: %6.4f, %6.4f" % (np.nanmedian(misClassZ1['blendedness.abs.flux']), 
                                                                                   np.nanmedian(galGoodZ1['blendedness.abs.flux'])))

misClassZ1m = galUseZ1m[(galUseZ1m['force.classification.extendedness'] <= 0.5)]
galGoodZ1m = galUseZ1m[(galUseZ1m['force.classification.extendedness'] > 0.5)]
print("# 1m: %d galaxies are mis-classified as extended objects" % len(misClassZ1m))
print("## median magnitude for misclassified and good galaxies: %5.3f, %5.3f" % (np.nanmedian(misClassZ1m['mag_Z']), 
                                                                                 np.nanmedian(galGoodZ1m['mag_Z'])))
print("## median blendedness for misclassified and good galaxies: %6.4f, %6.4f" % (np.nanmedian(misClassZ1m['blendedness.abs.flux']), 
                                                                                   np.nanmedian(galGoodZ1m['blendedness.abs.flux'])))

misClassZ2 = galUseZ2[(galUseZ2['force.classification.extendedness'] <= 0.5)]
galGoodZ2 = galUseZ2[(galUseZ2['force.classification.extendedness'] > 0.5)]
print("#  2: %d galaxies are mis-classified as extended objects" % len(misClassZ2))
print("## median magnitude for misclassified and good galaxies: %5.3f, %5.3f" % (np.nanmedian(misClassZ2['mag_Z']), 
                                                                                 np.nanmedian(galGoodZ2['mag_Z'])))
print("## median blendedness for misclassified and good galaxies: %6.4f, %6.4f" % (np.nanmedian(misClassZ2['blendedness.abs.flux']), 
                                                                                   np.nanmedian(galGoodZ2['blendedness.abs.flux'])))



    


















    






# HSC-Z band
#  1: 303 galaxies are mis-classified as extended objects
## median magnitude for misclassified and good galaxies: 24.632, 23.870
## median blendedness for misclassified and good galaxies: 0.0016, 0.0000
# 1m: 248 galaxies are mis-classified as extended objects
## median magnitude for misclassified and good galaxies: 24.668, 23.872
## median blendedness for misclassified and good galaxies: 0.0023, 0.0000
#  2: 1729 galaxies are mis-classified as extended objects
## median magnitude for misclassified and good galaxies: 24.536, 23.801
## median blendedness for misclassified and good galaxies: 0.0003, 0.0001

















In [22]:



    


print("# HSC-Y band")
misClassY1 = galUseY1[(galUseY1['force.classification.extendedness'] <= 0.5)]
galGoodY1 = galUseY1[(galUseY1['force.classification.extendedness'] > 0.5)]
print("#  1: %d galaxies are mis-classified as extended objects" % len(misClassY1))
print("## median magnitude for misclassified and good galaxies: %5.3f, %5.3f" % (np.nanmedian(misClassY1['mag_Y']), 
                                                                                 np.nanmedian(galGoodY1['mag_Y'])))
print("## median blendedness for misclassified and good galaxies: %6.4f, %6.4f" % (np.nanmedian(misClassY1['blendedness.abs.flux']), 
                                                                                   np.nanmedian(galGoodY1['blendedness.abs.flux'])))

misClassY1m = galUseY1m[(galUseY1m['force.classification.extendedness'] <= 0.5)]
galGoodY1m = galUseY1m[(galUseY1m['force.classification.extendedness'] > 0.5)]
print("# 1m: %d galaxies are mis-classified as extended objects" % len(misClassY1m))
print("## median magnitude for misclassified and good galaxies: %5.3f, %5.3f" % (np.nanmedian(misClassY1m['mag_Y']), 
                                                                                 np.nanmedian(galGoodY1m['mag_Y'])))
print("## median blendedness for misclassified and good galaxies: %6.4f, %6.4f" % (np.nanmedian(misClassY1m['blendedness.abs.flux']), 
                                                                                   np.nanmedian(galGoodY1m['blendedness.abs.flux'])))

misClassY2 = galUseY2[(galUseY2['force.classification.extendedness'] <= 0.5)]
galGoodY2 = galUseY2[(galUseY2['force.classification.extendedness'] > 0.5)]
print("#  2: %d galaxies are mis-classified as extended objects" % len(misClassY2))
print("## median magnitude for misclassified and good galaxies: %5.3f, %5.3f" % (np.nanmedian(misClassY2['mag_Y']), 
                                                                                 np.nanmedian(galGoodY2['mag_Y'])))
print("## median blendedness for misclassified and good galaxies: %6.4f, %6.4f" % (np.nanmedian(misClassY2['blendedness.abs.flux']), 
                                                                                   np.nanmedian(galGoodY2['blendedness.abs.flux'])))



    


















    






# HSC-Y band
#  1: 551 galaxies are mis-classified as extended objects
## median magnitude for misclassified and good galaxies: 24.542, 23.535
## median blendedness for misclassified and good galaxies: 0.0339, 0.0006
# 1m: 526 galaxies are mis-classified as extended objects
## median magnitude for misclassified and good galaxies: 24.568, 23.535
## median blendedness for misclassified and good galaxies: 0.0381, 0.0006
#  2: 1606 galaxies are mis-classified as extended objects
## median magnitude for misclassified and good galaxies: 24.336, 23.497
## median blendedness for misclassified and good galaxies: 0.0109, 0.0008

Write the catalogs





In [23]:



    


galGoodG1.write(os.path.join(galaxyDir, 'gal_all1_HSC-G_good.fits'), format='fits', overwrite=True)
galGoodG2.write(os.path.join(galaxyDir, 'gal_all2_HSC-G_good.fits'), format='fits', overwrite=True)

misClassG1.write(os.path.join(galaxyDir, 'gal_all1_HSC-G_misclass.fits'), format='fits', overwrite=True)
misClassG2.write(os.path.join(galaxyDir, 'gal_all2_HSC-G_misclass.fits'), format='fits', overwrite=True)

galGoodR1.write(os.path.join(galaxyDir, 'gal_all1_HSC-R_good.fits'), format='fits', overwrite=True)
galGoodR2.write(os.path.join(galaxyDir, 'gal_all2_HSC-R_good.fits'), format='fits', overwrite=True)

misClassR1.write(os.path.join(galaxyDir, 'gal_all1_HSC-R_misclass.fits'), format='fits', overwrite=True)
misClassR2.write(os.path.join(galaxyDir, 'gal_all2_HSC-R_misclass.fits'), format='fits', overwrite=True)

galGoodI1.write(os.path.join(galaxyDir, 'gal_all1_HSC-I_good.fits'), format='fits', overwrite=True)
galGoodI2.write(os.path.join(galaxyDir, 'gal_all2_HSC-I_good.fits'), format='fits', overwrite=True)

misClassI1.write(os.path.join(galaxyDir, 'gal_all1_HSC-I_misclass.fits'), format='fits', overwrite=True)
misClassI2.write(os.path.join(galaxyDir, 'gal_all2_HSC-I_misclass.fits'), format='fits', overwrite=True)

galGoodZ1.write(os.path.join(galaxyDir, 'gal_all1_HSC-Z_good.fits'), format='fits', overwrite=True)
galGoodZ2.write(os.path.join(galaxyDir, 'gal_all2_HSC-Z_good.fits'), format='fits', overwrite=True)

misClassZ1.write(os.path.join(galaxyDir, 'gal_all1_HSC-Z_misclass.fits'), format='fits', overwrite=True)
misClassZ2.write(os.path.join(galaxyDir, 'gal_all2_HSC-Z_misclass.fits'), format='fits', overwrite=True)

galGoodY1.write(os.path.join(galaxyDir, 'gal_all1_HSC-Y_good.fits'), format='fits', overwrite=True)
galGoodY2.write(os.path.join(galaxyDir, 'gal_all2_HSC-Y_good.fits'), format='fits', overwrite=True)

misClassY1.write(os.path.join(galaxyDir, 'gal_all1_HSC-Y_misclass.fits'), format='fits', overwrite=True)
misClassY2.write(os.path.join(galaxyDir, 'gal_all2_HSC-Y_misclass.fits'), format='fits', overwrite=True)

galGoodG1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-G_good_multipix.fits'), format='fits', overwrite=True)
misClassG1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-G_misclass_multipix.fits'), format='fits', overwrite=True)

galGoodR1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-R_good_multipix.fits'), format='fits', overwrite=True)
misClassR1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-R_misclass_multipix.fits'), format='fits', overwrite=True)

galGoodI1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-I_good_multipix.fits'), format='fits', overwrite=True)
misClassI1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-I_misclass_multipix.fits'), format='fits', overwrite=True)

galGoodZ1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-Z_good_multipix.fits'), format='fits', overwrite=True)
misClassZ1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-Z_misclass_multipix.fits'), format='fits', overwrite=True)

galGoodY1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-Y_good_multipix.fits'), format='fits', overwrite=True)
misClassY1m.write(os.path.join(galaxyDir, 'gal_all1_HSC-Y_misclass_multipix.fits'), format='fits', overwrite=True)

Read the catalogs





In [13]:



    


galGoodG1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-G_good.fits'), format='fits')
galGoodG2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-G_good.fits'), format='fits')

misClassG1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-G_misclass.fits'), format='fits')
misClassG2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-G_misclass.fits'), format='fits')

galGoodR1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-R_good.fits'), format='fits')
galGoodR2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-R_good.fits'), format='fits')

misClassR1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-R_misclass.fits'), format='fits')
misClassR2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-R_misclass.fits'), format='fits')

galGoodI1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-I_good.fits'), format='fits')
galGoodI2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-I_good.fits'), format='fits')

misClassI1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-I_misclass.fits'), format='fits')
misClassI2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-I_misclass.fits'), format='fits')

galGoodZ1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-Z_good.fits'), format='fits')
galGoodZ2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-Z_good.fits'), format='fits')

misClassZ1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-Z_misclass.fits'), format='fits')
misClassZ2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-Z_misclass.fits'), format='fits')

galGoodY1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-Y_good.fits'), format='fits')
galGoodY2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-Y_good.fits'), format='fits')

misClassY1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-Y_misclass.fits'), format='fits')
misClassY2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-Y_misclass.fits'), format='fits')

galGoodG1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-G_good_multipix.fits'), format='fits')
galGoodR1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-R_good_multipix.fits'), format='fits')
galGoodI1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-I_good_multipix.fits'), format='fits')
galGoodZ1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-Z_good_multipix.fits'), format='fits')
galGoodY1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-Y_good_multipix.fits'), format='fits')

misClassG1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-G_misclass_multipix.fits'), format='fits')
misClassR1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-R_misclass_multipix.fits'), format='fits')
misClassI1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-I_misclass_multipix.fits'), format='fits')
misClassZ1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-Z_misclass_multipix.fits'), format='fits')
misClassY1m = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-Y_misclass_multipix.fits'), format='fits')



    











In [14]:



    


galGoodG1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-G_good.fits'), format='fits')
galGoodG2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-G_good.fits'), format='fits')

galGoodI1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-I_good.fits'), format='fits')
galGoodI2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-I_good.fits'), format='fits')

misClassI1 = Table.read(os.path.join(galaxyDir, 'gal_all1_HSC-I_misclass.fits'), format='fits')
misClassI2 = Table.read(os.path.join(galaxyDir, 'gal_all2_HSC-I_misclass.fits'), format='fits')

Basic distributions of magnitudes and blendedness

Distributions of input galaxies





In [59]:



    


xtickFormat, ytickFormat = r'$\mathbf{%4.1f}$', r'$\mathbf{%4.1f}\ $'
ytickFormat2 = r'$\mathbf{%3.1f}$'

# --------------------------------------------------------------------------------------- #
## Setup up figure
fig = plt.figure(figsize=(30,6))
fig.subplots_adjust(left=0.02, right=0.99, bottom=0.15, top=0.99, 
                    hspace=0.0, wspace=0.0)

# --------------------------------------------------------------------------------------- #
## Ax1
ax1 = fig.add_subplot(151)
ax1.grid()
ax1 = songPlotSetup(ax1, ylabel=25, xlabel=25, border=4.0,
                    xtickFormat=xtickFormat, ytickFormat=ytickFormat)

## X, Y limits
xmin, xmax = 17.1, 27.99
ymin, ymax = -0.001, 0.589

# Good seeing
xx = galGoodI1['mag_G']
kde = KernelDensity(kernel='gaussian', bandwidth=0.05).fit(xx[:, None])
xPlot = np.linspace(xmin, 28.88, 500)
log_dens = kde.score_samples(xPlot[:, None])
ax1.fill(xPlot[:, None], np.exp(log_dens), facecolor=GCMAP(0.5), edgecolor=GCMAP(0.8), 
         alpha=0.6, linewidth=4.0, label=r'$\mathrm{Good\ Seeing}$')

# Bad seeing
xx = galGoodI2['mag_G']
kde = KernelDensity(kernel='gaussian', bandwidth=0.05).fit(xx[:, None])
xPlot = np.linspace(xmin, xmax, 500)
log_dens = kde.score_samples(xPlot[:, None])
ax1.plot(xPlot, np.exp(log_dens), color=GCMAP(0.9), alpha=0.8, linewidth=6.0, 
         linestyle='--', dashes=(40,8), label=r'$\mathrm{Bad\ Seeing}$')

ax1.yaxis.set_major_formatter(NullFormatter())

xlabel = r'$g\ \mathrm{Band\ Input\ Magnitude\ (mag)}$'
ylabel = r'$\#$'
ax1.set_xlabel(xlabel, size=26)
ax1.set_ylabel(ylabel, size=26)

ax1.set_xlim(xmin, xmax)
ax1.set_ylim(ymin, ymax)

# --------------------------------------------------------------------------------------- #
## Legend
l_handles, l_labels = ax1.get_legend_handles_labels()
ax1.legend(l_handles, l_labels, loc=(0.04, 0.72), shadow=True, fancybox=True, 
           numpoints=1, fontsize=24, scatterpoints=1, 
           markerscale=1.2, borderpad=0.2, handletextpad=0.2)
legend = ax1.get_legend()
legend.legendHandles[0].set_color(GCMAP(0.6))
legend.legendHandles[1].set_color(GCMAP(0.9))
# --------------------------------------------------------------------------------------- #

# --------------------------------------------------------------------------------------- #
## Ax2
ax2 = fig.add_subplot(152)
ax2.grid()
ax2 = songPlotSetup(ax2, ylabel=25, xlabel=25, border=4.0,
                    xtickFormat=xtickFormat, ytickFormat=ytickFormat)

# Good seeing
xx = galGoodI1['mag_R']
kde = KernelDensity(kernel='gaussian', bandwidth=0.05).fit(xx[:, None])
xPlot = np.linspace(xmin, xmax, 500)
log_dens = kde.score_samples(xPlot[:, None])
ax2.fill(xPlot[:, None], np.exp(log_dens), facecolor=RCMAP(0.5), edgecolor=RCMAP(0.8), 
         alpha=0.6, linewidth=4.0)

# Bad seeing
xx = galGoodI2['mag_R']
kde = KernelDensity(kernel='gaussian', bandwidth=0.05).fit(xx[:, None])
xPlot = np.linspace(xmin, xmax, 500)
log_dens = kde.score_samples(xPlot[:, None])
ax2.plot(xPlot, np.exp(log_dens), color=RCMAP(0.9), alpha=0.8, linewidth=6.0, 
         linestyle='--', dashes=(40,8))

ax2.set_ylim(ax1.get_ylim())
ax2.yaxis.set_major_formatter(NullFormatter())

xlabel = r'$r\ \mathrm{Band\ Input\ Magnitude\ (mag)}$'
ax2.set_xlabel(xlabel, size=26)

## X, Y limits
xmin, xmax = 17.2, 26.9
ax2.set_xlim(xmin, xmax)

# --------------------------------------------------------------------------------------- #
## Ax3
ax3 = fig.add_subplot(153)
ax3.grid()
ax3 = songPlotSetup(ax3, ylabel=25, xlabel=25, border=4.0,
                    xtickFormat=xtickFormat, ytickFormat=ytickFormat)

# Good seeing
xx = galGoodI1['mag_I']
kde = KernelDensity(kernel='gaussian', bandwidth=0.05).fit(xx[:, None])
xPlot = np.linspace(xmin, xmax, 500)
log_dens = kde.score_samples(xPlot[:, None])
ax3.fill(xPlot[:, None], np.exp(log_dens), facecolor=ICMAP(0.5), edgecolor=ICMAP(0.8), 
         alpha=0.6, linewidth=4.0)

# Bad seeing
xx = galGoodI2['mag_I']
kde = KernelDensity(kernel='gaussian', bandwidth=0.05).fit(xx[:, None])
xPlot = np.linspace(xmin, xmax, 500)
log_dens = kde.score_samples(xPlot[:, None])
ax3.plot(xPlot, np.exp(log_dens), color=ICMAP(0.9), alpha=0.8, linewidth=6.0, 
         linestyle='--', dashes=(40,8))

ax3.set_ylim(ax1.get_ylim())
ax3.yaxis.set_major_formatter(NullFormatter())

xlabel = r'$i\ \mathrm{Band\ Input\ Magnitude\ (mag)}$'
ax3.set_xlabel(xlabel, size=26)

## X, Y limits
xmin, xmax = 17.2, 26.9
ax3.set_xlim(xmin, xmax)

# --------------------------------------------------------------------------------------- #
## Ax4
ax4 = fig.add_subplot(154)
ax4.grid()
ax4 = songPlotSetup(ax4, ylabel=25, xlabel=25, border=4.0,
                    xtickFormat=xtickFormat, ytickFormat=ytickFormat)

# Good seeing
xx = galGoodI1['mag_Z']
kde = KernelDensity(kernel='gaussian', bandwidth=0.05).fit(xx[:, None])
xPlot = np.linspace(xmin, xmax, 500)
log_dens = kde.score_samples(xPlot[:, None])
ax4.fill(xPlot[:, None], np.exp(log_dens), facecolor=ZCMAP(0.5), edgecolor=ZCMAP(0.8), 
         alpha=0.6, linewidth=4.0)

# Bad seeing
xx = galGoodI2['mag_Z']
kde = KernelDensity(kernel='gaussian', bandwidth=0.05).fit(xx[:, None])
xPlot = np.linspace(xmin, xmax, 500)
log_dens = kde.score_samples(xPlot[:, None])
ax4.plot(xPlot, np.exp(log_dens), color=ZCMAP(0.9), alpha=0.8, linewidth=6.0, 
         linestyle='--', dashes=(40,8))

ax4.set_ylim(ax1.get_ylim())
ax4.yaxis.set_major_formatter(NullFormatter())

xlabel = r'$z\ \mathrm{Band\ Input\ Magnitude\ (mag)}$'
ax4.set_xlabel(xlabel, size=26)

## X, Y limits
xmin, xmax = 17.2, 26.9
ax4.set_xlim(xmin, xmax)

# --------------------------------------------------------------------------------------- #
## Ax5
ax5 = fig.add_subplot(155)
ax5.grid()
ax5 = songPlotSetup(ax5, ylabel=25, xlabel=25, border=4.0,
                    xtickFormat=xtickFormat, ytickFormat=ytickFormat)

# Good seeing
xx = galGoodI1['mag_Y']
kde = KernelDensity(kernel='gaussian', bandwidth=0.05).fit(xx[:, None])
xPlot = np.linspace(xmin, xmax, 500)
log_dens = kde.score_samples(xPlot[:, None])
ax5.fill(xPlot[:, None], np.exp(log_dens), facecolor=YCMAP(0.5), edgecolor=YCMAP(0.8), 
         alpha=0.6, linewidth=4.0)

# Bad seeing
xx = galGoodI2['mag_Y']
kde = KernelDensity(kernel='gaussian', bandwidth=0.05).fit(xx[:, None])
xPlot = np.linspace(xmin, xmax, 500)
log_dens = kde.score_samples(xPlot[:, None])
ax5.plot(xPlot, np.exp(log_dens), color=YCMAP(0.9), alpha=0.8, linewidth=6.0, 
         linestyle='--', dashes=(40,8))

ax5.set_ylim(ax1.get_ylim())
ax5.yaxis.set_major_formatter(NullFormatter())

xlabel = r'$Y\ \mathrm{Band\ Input\ Magnitude\ (mag)}$'
ax5.set_xlabel(xlabel, size=26)

## X, Y limits
xmin, xmax = 17.2, 26.9
ax5.set_xlim(xmin, xmax)

# --------------------------------------------------------------------------------------- #
fig.savefig(os.path.join(figDir, 'gal1_hist_input_mag.pdf'), 
            format='pdf', dpi=90)



    


















    
























In [18]:



    


# --------------------------------------------------------------------------------------- #
## Prepare the data

xmin, xmax = -0.79, 2.79
ymin, ymax = -0.59, 1.79

xx = inputI1['mag_G'] - inputI1['mag_R']
yy = inputI1['mag_R'] - inputI1['mag_I']
xTemp, yTemp = copy.deepcopy(xx), copy.deepcopy(yy)
xx = xx[(xTemp >= xmin) & (xTemp <= xmax) &
        (yTemp >= ymin) & (yTemp <= ymax)]
yy = yy[(xTemp >= xmin) & (xTemp <= xmax) &
        (yTemp >= ymin) & (yTemp <= ymax)]

xx2 = inputI2['mag_G'] - inputI2['mag_R']
yy2 = inputI2['mag_R'] - inputI2['mag_I']
xTemp, yTemp = copy.deepcopy(xx2), copy.deepcopy(yy2)
xx2 = xx2[(xTemp >= xmin) & (xTemp <= xmax) &
          (yTemp >= ymin) & (yTemp <= ymax)]
yy2 = yy2[(xTemp >= xmin) & (xTemp <= xmax) &
          (yTemp >= ymin) & (yTemp <= ymax)]

xlabel = r'$(g-r)_{\mathrm{Input}}} \mathrm{(mag)}$'
ylabel = r'$(r-i)_{\mathrm{Input}}} \mathrm{(mag)}$'

xtickFormat, ytickFormat = r'$\mathbf{%5.1f}$', r'$\mathbf{%5.1f}\ $'
ytickFormat2 = r'$\mathbf{%3.1f}$'

# --------------------------------------------------------------------------------------- #
## Setup up figure
fig = plt.figure(figsize=(15, 15))
ax1 = plt.axes(EC1)
ax2 = plt.axes(EC2)
ax3 = plt.axes(EC3)

ax1 = songPlotSetup(ax1, ylabel=50, xlabel=50, border=7.0,
                    xtickFormat=xtickFormat, ytickFormat=ytickFormat)
ax2 = songPlotSetup(ax2, ylabel=50, xlabel=50, border=7.0, 
                    xtickFormat=xtickFormat, ytickFormat=ytickFormat2)
ax3 = songPlotSetup(ax3, ylabel=50, xlabel=50, border=7.0, 
                    xtickFormat=xtickFormat, ytickFormat=ytickFormat2)
# --------------------------------------------------------------------------------------- #

# --------------------------------------------------------------------------------------- #
## Ax1
ax1.grid()

# HexBin
HB = ax1.hexbin(xx, yy, cmap=RCMAP, mincnt=2,
                alpha=0.7, gridsize=[50, 30], label=r'$\mathrm{Good\ Seeing}$',
                marginals=False, vmin=10, vmax=300, 
                edgecolor=RCMAP(0.3))

## Contour
H, xbins, ybins = np.histogram2d(xx2, yy2, bins=(200, 200))
H = gaussian_filter(H, 2)
levels = np.linspace(1, H.max(), 6)
Nx, Ny = len(xbins), len(ybins)
CS = ax1.contour(H.T, extent=[xbins[0], xbins[-1], ybins[0], ybins[-1]], 
                 colors=(RCMAP(0.6), RCMAP(0.6)),
                 linewidths=6, levels=levels, alpha=0.7)

# --------------------------------------------------------------------------------------- #
ax1.set_xlabel(xlabel, size=54)
ax1.set_ylabel(ylabel, size=54)

## X, Y limits
ax1.set_xlim(xmin, xmax)
ax1.set_ylim(ymin, ymax)

# --------------------------------------------------------------------------------------- #
## Legend
l_handles, l_labels = ax1.get_legend_handles_labels()
l_handles[0] = mpatches.RegularPolygon((0,0), 6)
l_handles.append(CS.collections[-1])
l_labels[0] = '$\mathrm{Good\ Seeing}$'
l_labels.append('$\mathrm{Bad\ Seeing}$')
ax1.legend(l_handles, l_labels, loc=(0.52, 0.05), shadow=True, fancybox=True, 
           numpoints=1, fontsize=40, scatterpoints=1, 
           markerscale=1.2, borderpad=0.2, handletextpad=0.3)
legend = ax1.get_legend()
legend.legendHandles[0].set_color(RCMAP(0.3))
legend.legendHandles[1].set_color(RCMAP(0.8))
# --------------------------------------------------------------------------------------- #

# --------------------------------------------------------------------------------------- #
## Ax2

ax2.grid()
# Gaussian KDE
kde = KernelDensity(kernel='gaussian', bandwidth=0.015).fit(xx[:, None])
xPlot = np.linspace(xmin, xmax, 500)
log_dens = kde.score_samples(xPlot[:, None])
ax2.fill(xPlot[:, None], np.exp(log_dens), facecolor=RCMAP(0.3), edgecolor=RCMAP(0.6), 
         alpha=0.7, linewidth=4.0)

kde = KernelDensity(kernel='gaussian', bandwidth=0.015).fit(xx2[:, None])
log_dens = kde.score_samples(xPlot[:, None])
ax2.fill(xPlot[:, None], np.exp(log_dens), facecolor='None', edgecolor=RCMAP(0.9), 
         alpha=0.6, linewidth=7.0, linestyle='-')

ax2.set_xlim(ax1.get_xlim())
ax2.yaxis.set_major_formatter(NullFormatter())
ax2.xaxis.set_major_formatter(NullFormatter())
# --------------------------------------------------------------------------------------- #

# --------------------------------------------------------------------------------------- #
## Ax3

ax3.grid()
# Gaussian KDE
kde = KernelDensity(kernel='gaussian', bandwidth=0.015).fit(yy[:, None])
yPlot = np.linspace(ymin, ymax, 500)
log_dens = kde.score_samples(yPlot[:, None])
ax3.fill_betweenx(yPlot, 0.0, np.exp(log_dens), facecolor=RCMAP(0.3), edgecolor=RCMAP(0.6),
                  alpha=0.7, linewidth=4.0)

kde = KernelDensity(kernel='gaussian', bandwidth=0.015).fit(yy2[:, None])
log_dens = kde.score_samples(yPlot[:, None])
ax3.fill_betweenx(yPlot, 0.0, np.exp(log_dens), facecolor='None', edgecolor=RCMAP(0.9),
                  alpha=0.6, linewidth=7.0, linestyle='-')

ax3.set_ylim(ax1.get_ylim())
ax3.yaxis.set_major_formatter(NullFormatter())
ax3.xaxis.set_major_formatter(NullFormatter())
# --------------------------------------------------------------------------------------- #

fig.savefig(os.path.join(figDir, 'galaxy1_HSC-G-R-input.pdf'), 
            format='pdf', dpi=90)