In [2]:
%load_ext autoreload
%reload_ext autoreload
%autoreload 2

% matplotlib inline
from __future__ import (division, 
                        print_function)

import os
import sys
import copy
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
import cPickle as pickle
from scipy import ndimage

# Astropy
from astropy.io import fits
from astropy    import units as u
from astropy.stats import sigma_clip
from astropy.table import Table, Column
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.pyplot as plt
import matplotlib.patches as mpatches
import matplotlib.gridspec as gridspec
from matplotlib.patches import Ellipse
from matplotlib.ticker import NullFormatter, MaxNLocator, FormatStrFormatter
from matplotlib.collections import PatchCollection
tickFormat = FormatStrFormatter('$\mathbf{%g}$') 

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

# SEP 
import sep

# 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
import galSBP as sbp
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 = [left,   bottom, width, height]
recHist = [right,  bottom,  0.20, height]
SBP1 = [0.124, 0.085, 0.865, 0.33]
SBP2 = [0.124, 0.41, 0.865, 0.55]
EC1 = [0.135, 0.072, 0.862, 0.295]
EC2 = [0.135, 0.366, 0.862, 0.295]
EC3 = [0.135, 0.666, 0.862, 0.295]
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'

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


# Cubehelix
from palettable.cubehelix import Cubehelix
CUBE = Cubehelix.make(start=0.3, rotation=-0.5, n=16).mpl_colormap

import itertools


The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

Catalog and File Names


In [3]:
cat = Table.read('catalog_highres_z0.37.fits', format='fits')

mstar = 'stellar_mass_maps'
mtot = 'total_mass_maps'

print(cat.colnames)


['id', 'id_subhalo', 'central', 'logM500', 'r500', 'logM200', 'r200', 'logM30_3D', 'logM100_3D']

Estimate the average shape and ellipticity profiles


In [4]:
#isophote = '/iraf/iraf/extern/stsdas/bin.macosx/x_isophote.e'
#xttools = '/iraf/iraf/extern/tables/bin.macosx/x_ttools.e'
isophote = '/Users/song/anaconda/pkgs/iraf-2.16.1-0/variants/common/iraf/stsci_iraf/stsdas/bin.macosx/x_isophote.e'
xttools = '/Users/song/anaconda/pkgs/iraf-2.16.1-0/variants/common/iraf/stsci_iraf/tables/bin.macosx/x_ttools.e'

# The Gaussian kernal used for convolution
kernel1 = np.asarray([[0.560000, 0.980000, 0.560000],
                      [0.980000, 1.000000, 0.980000],
                      [0.560000, 0.980000, 0.560000]])

kernel2 = np.asarray([[0.000000, 0.220000, 0.480000, 0.220000, 0.000000],
                      [0.220000, 0.990000, 1.000000, 0.990000, 0.220000],
                      [0.480000, 1.000000, 1.000000, 1.000000, 0.480000],
                      [0.220000, 0.990000, 1.000000, 0.990000, 0.220000],
                      [0.000000, 0.220000, 0.480000, 0.220000, 0.000000]])

kernel3 = np.asarray([[0.092163, 0.221178, 0.296069, 0.221178, 0.092163],
                      [0.221178, 0.530797, 0.710525, 0.530797, 0.221178],
                      [0.296069, 0.710525, 0.951108, 0.710525, 0.296069],
                      [0.221178, 0.530797, 0.710525, 0.530797, 0.221178],
                      [0.092163, 0.221178, 0.296069, 0.221178, 0.092163]])

In [17]:
len(cat) / 7


Out[17]:
19.0

Center and Average Shape Using SEP


In [5]:
# Prepare the plot
fig = plt.figure(figsize=(14, 38))
fig.subplots_adjust(left=0.0, right=1.0, 
                    bottom=0.0, top=1.0,
                    wspace=0.00, hspace=0.00)
gs = gridspec.GridSpec(19, 7)
gs.update(wspace=0.0, hspace=0.00)

xArr, yArr, baArr, paArr, m10Arr, m30Arr, m100Arr, mtotArr = [], [], [], [], [], [], [], []
aperProf, aperMass = [], []
r50Arr = []

with ProgressBar(len(cat)) as bar:
    for item in enumerate(cat):
        ii, col = item
        galID = str(col['id'])
        bar.update()
        #print("## Deal with galaxy : %s" % galID)

        # Convert the mass density map to just mass maps
        imgTest = os.path.join(mstar, 'galaxy_' + galID + '_stars_no_sats.fits')
        imgSave = os.path.join(mstar, 'galaxy_' + galID + '_stars_no_sats_ms.fits')

        imgN = fits.open(imgTest)[0].data / 1.0E6 * 25.0 
        hdu = fits.PrimaryHDU(imgN)
        hdulist = fits.HDUList([hdu])
        hdulist.writeto(imgSave, clobber=True)

        # Extract objects
        bkg = sep.Background(imgN, bw=20, bh=20, fw=4, fh=4)
        imgSub = imgN - bkg
        objects = sep.extract(imgSub, 15.0, filter_kernel=kernel3)

        # Find the object at the center
        if len(objects) == 1: 
            index = 0
        else: 
            index = np.argmin(np.sqrt((objects['x'] - 50.0) ** 2.0 + (objects['y'] - 50.0) ** 2.0))

        #print(len(objects), index)

        # Get the naive ba, theta, xcen, ycen
        ba = objects[index]['b'] / objects[index]['a']
        theta = objects[index]['theta']
        xcen, ycen = objects[index]['x'], objects[index]['y']
        baArr.append(ba)
        paArr.append(theta * 180.0 / np.pi)
        xArr.append(xcen)
        yArr.append(ycen)

        # Aperture photometry and naive "profiles"
        rad = np.linspace(1.2, 4.0, 50)
        kpc = (rad ** 4.0)
        pix = (kpc / 5.0)
        area = (np.pi * (kpc ** 2.0) * ba)
        maper = sep.sum_ellipse(imgN, xcen, ycen, pix, pix * ba, theta, 1.0, 
                                bkgann=None, subpix=11)[0]
        mring = maper[1:] - maper[0:-1]
        aring = area[1:] - area[0:-1]
        rho = np.log10(mring / aring)
        aperMass.append(maper)
        aperProf.append(rho)
        
        # Mass within different apertures
        rad = np.asarray([10.0, 30.0, 100.0, 250.0]) / 5.0
        maper = sep.sum_ellipse(imgN, xcen, ycen, rad, rad * ba, theta, 1.0, 
                                bkgann=None, subpix=11)[0]
        m10Arr.append(np.log10(maper[0]))
        m30Arr.append(np.log10(maper[1]))
        m100Arr.append(np.log10(maper[2]))
        mtotArr.append(np.log10(np.sum(imgN)))
        
        # Get r50 
        r50, flag = sep.flux_radius(imgN, xcen, ycen, 30.0, 0.5,
                                    subpix=11)
        r50Arr.append(r50 * 5.0)

        # Show the aperture
        ax = plt.subplot(gs[ii])
        ax.yaxis.set_major_formatter(NullFormatter())
        ax.xaxis.set_major_formatter(NullFormatter())

        m, s = np.mean(imgN), np.std(imgN)
        im = ax.imshow(imgN, interpolation='nearest', cmap=CUBE,
                       vmin=m-1.0*s, vmax=m+10*s, origin='lower')

        # plot an ellipse for each object
        for i in range(len(objects)):
            e = Ellipse(xy=(xcen, ycen),
                        width=20 * 2.0,
                        height=(20 * 2.0 * ba),
                        angle=theta * 180. / np.pi)
            e.set_facecolor('none')
            e.set_edgecolor('r')
            e.set_alpha(0.5)
            e.set_linewidth(1.0)
            ax.add_artist(e)
        ax.set_aspect('equal')
        
        # Show the 2D and 3D mass within 100 kpc 
        ax.text(0.5, 0.85, "3D: %5.2f" % col['logM100_3D'], 
                verticalalignment='center', horizontalalignment='center',
                fontsize=20.0, transform=ax.transAxes, weight='bold', 
                color='r', alpha=0.8)
        ax.text(0.5, 0.15, "2D: %5.2f" % np.log10(maper[2]), 
                verticalalignment='center', horizontalalignment='center',
                fontsize=20.0, transform=ax.transAxes, weight='bold', 
                color='r', alpha=0.8)
    
fig.savefig('bhamas_map_aperture_hres_z0.37.pdf', dpi=140)


/Users/song/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:60: RuntimeWarning: divide by zero encountered in log10