John Moustakas
Siena College
2018 October 5
The goal of this notebook is to use observations of a large number of luminous red galaxies (LRGs) with DECaLS/DR7 grzW1W2 photometry and ancillary spectroscopic redshifts to generate a set of spectroscopic templates suitable for spectroscopic simulations, targeting investigations, and other DESI-related projects.
The parent sample consists of 111,114 LRGs selected by Rongpu Zhou (Pittsburgh) with spectroscopic redshifts from SDSS/Main, (e)BOSS, AGES, DEEP2, GAMA, OzDES, VIPERS, VVDS, and COSMOS (precise medium-band photometric redshifts).
We cull this sample to 75,315 LRGs by applying the following additional cuts:
To ensure this sample is representative of the parent sample of >3.3M LRGs with precise photometric redshifts, we count the number of photometric LRGs near each LRG with a spectroscopic redshift (in redshift-apparent magnitude-color space) and use this below as a statistical weight.
With this sample in hand, we generate spectral templates using the following procedure:
We use the Bayesian spectral energy distribution (SED) modeling code iSEDfit and the CKC14z simple stellar population models of Conroy et al. to model the observed SEDs of the parent sample of galaxies. The outputs of this code are the best-fitting (maximum likelihood) spectroscopic template and the marginalized posterior distributions on a large number of quantities of interest, stored in the file legacysurvey_lrg_ckc14z_kroupa01_charlot_sfhgrid01.fits.
In addition, we compute K-corrected absolute magnitudes in the DECaLS and BASS+MzLS g-, r-, and z-band bandpasses and store them in the file legacysurvey_lrg_ckc14z_kroupa01_charlot_sfhgrid01_kcorr.z0.0.fits. The code that generates these files is desi_legacysurvey_lrg_isedfit.
In [1]:
import os
import numpy as np
from astropy.table import Table, hstack
from scipy.spatial import cKDTree as KDTree
In [2]:
import matplotlib.pyplot as plt
from corner import corner
In [3]:
plt.style.use('seaborn-talk')
%matplotlib inline
In [4]:
seed = 123
rand = np.random.RandomState(seed)
In [5]:
isedfit_version = 'v2.0'
templates_version = 'v2.1'
In [6]:
isedfit_dir = os.path.join(os.getenv('DESI_ROOT'), 'spectro', 'templates', 'lrg_templates', 'isedfit', isedfit_version)
In [7]:
metadata_outfile = os.path.join(isedfit_dir, 'lrg-templates-isedfit-{}.fits'.format(isedfit_version))
In [8]:
def read_isedfit():
"""Read the iSEDfit fitting results."""
parentfile = os.path.join(isedfit_dir, 'lrg_parent.fits.gz')
isedfile = os.path.join(isedfit_dir, 'legacysurvey_lrg_ckc14z_kroupa01_charlot_sfhgrid01.fits.gz')
kcorrfile = os.path.join(isedfit_dir, 'legacysurvey_lrg_ckc14z_kroupa01_charlot_sfhgrid01_kcorr.z0.0.fits.gz')
print('Reading {}'.format(parentfile))
parent = Table.read(parentfile)
print('Reading {}'.format(isedfile))
ised = Table.read(isedfile)
print('Reading {}'.format(kcorrfile))
kcorr = Table.read(kcorrfile)
snrmin = 3.0
chi2min = 10
keep = np.where(
(ised['CHI2'] < chi2min) *
(np.sum(ised['MAGGIES'] * np.sqrt(ised['IVARMAGGIES']) > snrmin, axis=1) == 5)
)[0]
print('Read {} galaxies with chi2 < {} and S/N > {} in all 5 photometric bands.'.format(
len(keep), chi2min, snrmin))
cat = dict()
cat['weight'] = len(keep) * parent['COUNT'][keep].data / np.sum(parent['COUNT'][keep].data)
cat['redshift'] = kcorr['Z'][keep].data
cat['Mstar'] = ised['MSTAR_50'][keep].data
cat['Mg'] = kcorr['ABSMAG'][keep, 0].data
cat['Mr'] = kcorr['ABSMAG'][keep, 1].data
cat['Mz'] = kcorr['ABSMAG'][keep, 2].data
cat['gr'] = cat['Mg'] - cat['Mr']
cat['rz'] = cat['Mr'] - cat['Mz']
with np.errstate(invalid='ignore'):
cat['grobs'] = -2.5 * np.log10( ised['MAGGIES'][keep, 0].data / ised['MAGGIES'][keep, 1].data )
cat['rzobs'] = -2.5 * np.log10( ised['MAGGIES'][keep, 1].data / ised['MAGGIES'][keep, 2].data )
cat['zW1obs'] = -2.5 * np.log10( ised['MAGGIES'][keep, 2].data / ised['MAGGIES'][keep, 3].data )
#mm = - 2.5 * np.log10(ised['MAGGIES'][keep, 0].data)
#_ = plt.hist(mm, bins=100)
kcorr.remove_columns(['Z', 'ISEDFIT_ID', 'MAGGIES', 'IVARMAGGIES'])
out = hstack([ised[keep], kcorr[keep]])
return cat, out
In [9]:
cat, ised = read_isedfit()
ngal = len(ised)
In [10]:
print(sorted(cat.keys()))
In [11]:
ised[:2]
Out[11]:
In [12]:
zlim, Mrlim, grlim, rzlim = (0.2, 1.1), (-21, -25), (0.2, 1.2), (0.2, 0.8)
In [13]:
def qa_rest():
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 10))
ax1.hexbin(cat['redshift'], cat['Mr'], mincnt=1, bins='log',
C=cat['weight'], reduce_C_function=np.sum,
cmap=plt.cm.get_cmap('RdYlBu'))
ax1.set_ylim(Mrlim)
ax1.set_xlim(zlim)
ax1.set_xlabel('Redshift')
ax1.set_ylabel(r'$M_{0.0r}$')
ax2.hexbin(cat['redshift'], cat['gr'], mincnt=1, bins='log',
C=cat['weight'], reduce_C_function=np.sum,
cmap=plt.cm.get_cmap('RdYlBu'))
ax2.set_xlim(zlim)
ax2.set_ylim(grlim)
ax2.set_xlabel('Redshift')
ax2.set_ylabel(r'$^{0.0}(g - r)$')
ax3.hexbin(cat['rz'], cat['gr'], mincnt=1, bins='log',
C=cat['weight'], reduce_C_function=np.sum,
cmap=plt.cm.get_cmap('RdYlBu'))
ax3.set_xlabel(r'$^{0.0}(r - z)$')
ax3.set_ylabel(r'$^{0.0}(g - r)$')
ax3.set_ylim(grlim)
ax3.set_xlim(rzlim)
hb = ax4.hexbin(cat['Mr'], cat['gr'], mincnt=1, bins='log',
C=cat['weight'], reduce_C_function=np.sum,
cmap=plt.cm.get_cmap('RdYlBu'))
ax4.set_xlabel(r'$M_{0.0r}$')
ax4.set_ylabel(r'$^{0.0}(g - r)$')
ax4.set_xlim(Mrlim)
ax4.set_ylim(grlim)
cax = fig.add_axes([0.88, 0.15, 0.02, 0.7])
fig.colorbar(hb, cax=cax, label=r'$\log_{10}$ (Weighted Number of Galaxies)')
for aa in (ax1, ax2, ax3, ax4):
aa.grid(True)
plt.subplots_adjust(wspace=0.2, right=0.85)
In [14]:
qa_rest()
In [15]:
grobslim, rzobslim = (0.5, 3), (0.8, 2.2)
In [16]:
def qa_obs():
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
ax1.hexbin(cat['rzobs'], cat['grobs'], mincnt=1, bins='log',
C=cat['weight'], reduce_C_function=np.sum,
cmap=plt.cm.get_cmap('RdYlBu'))
ax1.set_xlabel(r'$(r - z)_{\rm obs}$')
ax1.set_ylabel(r'$(g - r)_{\rm obs}$')
hb = ax2.hexbin(cat['rzobs'], cat['zW1obs'], mincnt=1, bins='log',
C=cat['weight'], reduce_C_function=np.sum,
cmap=plt.cm.get_cmap('RdYlBu'))
cb = plt.colorbar(hb)
cb.set_label(r'$\log_{10}$ (Number of Galaxies)')
ax2.set_ylabel(r'$(r - z)_{\rm obs}$')
ax2.set_xlabel(r'$(z - W1)_{\rm obs}$')
ax1.set_xlim(0.7, 2.3)
ax1.set_ylim(0.3, 3.0)
ax2.set_xlim(0.7, 2.3)
ax2.set_ylim(0.2, 2.5)
for aa in (ax1, ax2):
aa.grid(True)
In [17]:
qa_obs()
The question is: what minimum set templates can be used to describe this population of LRGs. This problem is NP hard, although there are methods for solving it (see, e.g., the SetCoverPy algorithm).
Here we simply choose the desired number of templates and choose them with probability equal to the (known) statistical weight.
In [18]:
ntemplates = 3000
In [19]:
iarch = rand.choice(ngal, ntemplates, p=cat['weight']/np.sum(cat['weight']))
resp = ngal * cat['weight'] / np.sum(cat['weight'])
print('Selected {} templates from {} galaxies.'.format(ntemplates, ngal))
In [20]:
def _markers():
size = 50 * (1+(resp - resp.min()) / resp.ptp())
shade = (cat['gr'][iarch] - cat['gr'][iarch].min()) / cat['gr'][iarch].ptp()
col = plt.cm.coolwarm(shade)
return size, col
In [21]:
def qa_templates():
"""Generate a color-color plot with the symbol size scaled by the responsibility.
"""
size, col = _markers()
fig, ((ax1, ax2), (ax3, ax4), (ax5, ax6)) = plt.subplots(3, 2, figsize=(15, 15))
ax1.scatter(cat['rzobs'], cat['grobs'], s=10, c='lightgray', alpha=0.5, edgecolor='k')
ax1.scatter(cat['rzobs'][iarch], cat['grobs'][iarch], c=col, marker='o', s=size,
alpha=0.9, edgecolor='k')
ax1.set_xlabel(r'$(r - z)_{\rm obs}$')
ax1.set_ylabel(r'$(g - r)_{\rm obs}$')
ax2.scatter(cat['zW1obs'], cat['rzobs'], s=10, c='lightgray', alpha=0.5, edgecolor='k')
ax2.scatter(cat['zW1obs'][iarch], cat['rzobs'][iarch], c=col, marker='o', s=size,
alpha=0.9, edgecolor='k')
ax2.set_ylabel(r'$(r - z)_{\rm obs}$')
ax2.set_xlabel(r'$(z - W1)_{\rm obs}$')
ax3.scatter(cat['rz'], cat['gr'], s=10, c='lightgray', alpha=0.5, edgecolor='k')
ax3.scatter(cat['rz'][iarch], cat['gr'][iarch], c=col, marker='o', s=size,
alpha=0.9, edgecolor='k')
ax3.set_xlabel(r'$^{0.0}(r - z)$')
ax3.set_ylabel(r'$^{0.0}(g - r)$')
ax3.set_ylim(grlim)
ax3.set_xlim(rzlim)
ax4.scatter(cat['Mr'], cat['gr'], s=10, c='lightgray', alpha=0.5, edgecolor='k')
ax4.scatter(cat['Mr'][iarch], cat['gr'][iarch], c=col, marker='o', s=size,
alpha=0.9, edgecolor='k')
ax4.set_xlabel(r'$M_{0.0r}$')
ax4.set_ylabel(r'$^{0.0}(g - r)$')
ax4.set_xlim(Mrlim)
ax4.set_ylim(grlim)
ax5.scatter(cat['redshift'], cat['Mr'], s=10, c='lightgray', alpha=0.5, edgecolor='k')
ax5.scatter(cat['redshift'][iarch], cat['Mr'][iarch], c=col, marker='o', s=size,
alpha=0.9, edgecolor='k')
ax5.set_xlabel('Redshift')
ax5.set_ylabel(r'$M_{0.0r}$')
ax5.set_xlim(zlim)
ax5.set_ylim(Mrlim)
ax6.scatter(cat['redshift'], cat['gr'], s=10, c='lightgray', alpha=0.5, edgecolor='k')
ax6.scatter(cat['redshift'][iarch], cat['gr'][iarch], c=col, marker='o', s=size,
alpha=0.9, edgecolor='k')
ax6.set_xlabel('Redshift')
ax6.set_ylabel(r'$^{0.0}(g - r)$')
ax6.set_xlim(zlim)
ax6.set_ylim(grlim)
for aa in (ax1, ax2, ax3, ax4, ax5, ax6):
aa.grid(True)
fig.subplots_adjust(wspace=0.2, hspace=0.3)
In [22]:
qa_templates()
SFR100 is the star formation rate (SFR) averaged over the previous 100 Myr. Note that templates with higher SFRs are somewhat over-represented in the template set because we choose templates with equal probability, not weighted by the number of galaxies in each cell.
In [23]:
def qa_template_properties():
size, col = _markers()
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(12, 10))
ax1.hist(ised['MSTAR_50'], weights=cat['weight'], density=True, bins=30, label='Full Sample')
ax1.hist(ised['MSTAR_50'][iarch], density=True, bins=30, alpha=0.5,
label='Templates')
ax1.set_xlabel(r'$\log_{10}\ (M\ /\ M_{\odot})$')
ax1.set_ylabel('Relative Fraction')
ax1.set_xlim((10, 12.5))
ax1.yaxis.set_ticks([])
ax2.hist(ised['SFR100_50'], weights=cat['weight'], density=True, bins=30, label='Full Sample')
ax2.hist(ised['SFR100_50'][iarch], density=True, bins=30, alpha=0.5,
label='Templates')
ax2.set_xlabel(r'$\log_{10}\ (SFR_{100}\ /\ M_{\odot}\ {\rm yr}^{-1})$')
ax2.set_xlim((-4.5, 2))
ax2.yaxis.set_ticks([])
ax2.legend(loc='upper right')
ax3.scatter(cat['redshift'], ised['MSTAR_50'], s=10, c='lightgray', alpha=0.5, edgecolor='k')
ax3.scatter(cat['redshift'][iarch], ised['MSTAR_50'][iarch], c=col, marker='o', s=size,
alpha=0.9, edgecolor='k')
ax3.set_xlabel('Redshift')
ax3.set_ylabel(r'$\log_{10}\ (M\ /\ M_{\odot})$')
ax3.set_xlim(zlim)
ax3.set_ylim((10, 12.5))
ax4.scatter(cat['redshift'], ised['SFR100_50'], s=10, c='lightgray', alpha=0.5, edgecolor='k')
ax4.scatter(cat['redshift'][iarch], ised['SFR100_50'][iarch], c=col, marker='o', s=size,
alpha=0.9, edgecolor='k')
ax4.set_xlabel('Redshift')
ax4.set_ylabel(r'$\log_{10}\ (SFR_{100}\ /\ M_{\odot}\ {\rm yr}^{-1})$')
ax4.set_xlim(zlim)
ax4.set_ylim((-4.5, 2))
plt.subplots_adjust(wspace=0.25, hspace=0.2)
In [24]:
qa_template_properties()
In [25]:
print('Writing {}'.format(metadata_outfile))
ised[iarch].write(metadata_outfile, overwrite=True)
The spectra have to built in IDL using the script build_lrg_basis_templates, which we assume has happend.
In [26]:
from desisim.io import read_basis_templates
In [27]:
flux, wave, meta = read_basis_templates('LRG')
nt = len(meta)
print('Number of templates = {}'.format(nt))
In [28]:
meta
Out[28]:
In [29]:
def plot_subset(nplot=25, ncol=5):
"""Plot a random sampling of the basis templates."""
nspec, npix = flux.shape
nrow = np.ceil(nplot / ncol).astype('int')
these = rand.choice(nspec, nplot, replace=False)
fig, ax = plt.subplots(nrow, ncol, figsize=(2.2*ncol, 2.2*nrow), sharey=True, sharex=True)
for thisax, indx in zip(ax.flat, these):
thisax.plot(wave, flux[indx, :])
thisax.text(0.95, 0.93, '{:0d}'.format(indx), ha='right',
va='top', transform=thisax.transAxes, fontsize=11)
thisax.xaxis.set_major_locator(plt.MaxNLocator(3))
thisax.set_xscale('log')
thisax.set_yscale('log')
fig.subplots_adjust(wspace=0.05, hspace=0.05)
In [30]:
plot_subset()
In [31]:
def zgrid_colors():
"""Compute the colors of the templates on a fixed redshift grid."""
from speclite import filters
filt = filters.load_filters('decam2014-g', 'decam2014-r', 'decam2014-z', 'wise2010-W1')
zmin, zmax, dz = 0.2, 1.2, 0.1
#zmin, zmax, dz = 0.0, 1.0, 0.1
nz = np.round( (zmax - zmin) / dz ).astype('i2')
print('Number of redshift points = {}'.format(nz))
cc = dict(
redshift = np.linspace(zmin, zmax, nz),
gr = np.zeros( (nt, nz) ),
rz = np.zeros( (nt, nz) ),
zW1 = np.zeros( (nt, nz), )
)
for iz, red in enumerate(cc['redshift']):
zwave = wave.astype('float') * (1 + red)
phot = filt.get_ab_maggies(flux, zwave, mask_invalid=False)
cc['gr'][:, iz] = -2.5 * np.log10( phot['decam2014-g'] / phot['decam2014-r'] )
cc['rz'][:, iz] = -2.5 * np.log10( phot['decam2014-r'] / phot['decam2014-z'] )
cc['zW1'][:, iz] = -2.5 * np.log10( phot['decam2014-z'] / phot['wise2010-W1'] )
return cc
In [32]:
%time cc = zgrid_colors()
In [33]:
grrange = (0.0, 3.0)
rzrange = (0.0, 2.5)
zW1range = (-1, 4)
ntspace = 20 # spacing between model curves
In [34]:
def grz(pngfile=None):
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
hb = ax1.hexbin(cat['rzobs'], cat['grobs'], bins='log',
mincnt=1, extent=rzrange+grrange,
C=cat['weight'], reduce_C_function=np.sum,
cmap=plt.cm.get_cmap('RdYlBu'))
ax1.set_xlabel(r'$(r - z)_{\rm obs}$')
ax1.set_ylabel(r'$(g - r)_{\rm obs}$')
ax1.set_xlim(rzrange)
ax1.set_ylim(grrange)
ax1.text(0.05, 0.9, 'Data', ha='left', va='bottom',
transform=ax1.transAxes, fontsize=14)
ax1.grid(True)
#cb = fig.colorbar(hb, ax=ax1)
#cb.set_label(r'log$_{10}$ (Number of Galaxies)')
for tt in np.arange(0, nt, ntspace):
ax2.plot(cc['rz'][tt, :], cc['gr'][tt, :], marker='s',
markersize=5, ls='-', alpha=0.5)
for tt in np.arange(0, nt, ntspace):
ax2.scatter(cc['rz'][tt, 0], cc['gr'][tt, 0], marker='o',
facecolors='none', s=80, edgecolors='k',
linewidth=2)
ax2.text(0.1, 0.1, 'z=0.2', ha='left', va='bottom',
transform=ax2.transAxes, fontsize=14)
ax2.text(0.05, 0.9, 'Models (z=0.2-1.2, dz=0.1)', ha='left', va='bottom',
transform=ax2.transAxes, fontsize=14)
ax2.set_xlim(rzrange)
ax2.set_ylim(grrange)
ax2.set_xlabel(r'$(r - z)_{\rm obs}$')
ax2.set_ylabel(r'$(g - r)_{\rm obs}$')
ax2.yaxis.set_label_position('right')
ax2.yaxis.tick_right()
ax2.grid(True)
plt.subplots_adjust(wspace=0.05)
if pngfile:
fig.savefig(pngfile)
In [35]:
grz()
In [36]:
def rzW1(pngfile=None):
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
hb = ax1.hexbin(cat['zW1obs'], cat['rzobs'], bins='log',
mincnt=1, extent=rzrange+grrange,
C=cat['weight'], reduce_C_function=np.sum,
cmap=plt.cm.get_cmap('RdYlBu'))
ax1.set_ylabel(r'$(r - z)_{\rm obs}$')
ax1.set_xlabel(r'$(z - W1)_{\rm obs}$')
ax1.set_xlim(zW1range)
ax1.set_ylim(rzrange)
ax1.text(0.05, 0.9, 'Data', ha='left', va='bottom',
transform=ax1.transAxes, fontsize=14)
ax1.grid(True)
#cb = fig.colorbar(hb, ax=ax1)
#cb.set_label(r'log$_{10}$ (Number of Galaxies)')
for tt in np.arange(0, nt, ntspace):
ax2.plot(cc['zW1'][tt, :], cc['rz'][tt, :], marker='s',
markersize=5, ls='-', alpha=0.5)
for tt in np.arange(0, nt, ntspace):
ax2.scatter(cc['zW1'][tt, 0], cc['rz'][tt, 0], marker='o',
facecolors='none', s=80, edgecolors='k',
linewidth=2)
ax2.text(0.03, 0.3, 'z=0.2', ha='left', va='bottom',
transform=ax2.transAxes, fontsize=14)
ax2.text(0.05, 0.9, 'Models (z=0.2-1.2, dz=0.1)', ha='left', va='bottom',
transform=ax2.transAxes, fontsize=14)
ax2.set_xlim(zW1range)
ax2.set_ylim(rzrange)
ax2.set_xlabel(r'$(z - W1)_{\rm obs}$')
ax2.set_ylabel(r'$(r - z)_{\rm obs}$')
ax2.yaxis.set_label_position('right')
ax2.yaxis.tick_right()
ax2.grid(True)
plt.subplots_adjust(wspace=0.05)
if pngfile:
fig.savefig(pngfile)
In [37]:
rzW1()
In [ ]: