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
In [3]:
cat = Table.read('catalog_highres_z0.37.fits', format='fits')
mstar = 'stellar_mass_maps'
mtot = 'total_mass_maps'
print(cat.colnames)
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]:
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)