LRG Templates

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:

  • 0.2 < z < 1.2
  • nobs_[g,r,z] >= 3
  • wisemask == False
  • brightstarinblob == False
  • measured flux in all five grzW1W2 bands

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:

  1. 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.

  1. Next, the notebook below documents how we select a "representative" subset of these galaxies to serve as our LRG template set.
  1. Finally, we generate the templates themselves using the IDL script build_lrg_basis_templates.

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

Establish the random and pathnames.


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

Read the iSEDfit and absolute magnitude catalogs.


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)


Reading /Users/ioannis/work/desi/spectro/templates/lrg_templates/isedfit/v2.0/lrg_parent.fits.gz
Reading /Users/ioannis/work/desi/spectro/templates/lrg_templates/isedfit/v2.0/legacysurvey_lrg_ckc14z_kroupa01_charlot_sfhgrid01.fits.gz
Reading /Users/ioannis/work/desi/spectro/templates/lrg_templates/isedfit/v2.0/legacysurvey_lrg_ckc14z_kroupa01_charlot_sfhgrid01_kcorr.z0.0.fits.gz
Read 48817 galaxies with chi2 < 10 and S/N > 3.0 in all 5 photometric bands.

In [10]:
print(sorted(cat.keys()))


['Mg', 'Mr', 'Mstar', 'Mz', 'gr', 'grobs', 'redshift', 'rz', 'rzobs', 'weight', 'zW1obs']

In [11]:
ised[:2]


Out[11]:
Table length=2
ISEDFIT_IDRADECZMAGGIES [5]IVARMAGGIES [5]BESTMAGGIES [5]CHUNKINDXMODELINDXDELAYEDBURSTTYPECHI2_1TOTALMASSTOTALMASS_ERRMSTARAGESFRAGETAUZMETALAVMUOIIIHBNLYCSFRSFR100B100B1000EWOIIEWOIIIHBEWNIIHANBURSTTRUNCTAUTBURSTDTBURSTFBURSTMSTAR_50AGE_50SFRAGE_50TAU_50ZMETAL_50AV_50MU_50OIIIHB_50SFR_50SFR100_50B100_50B1000_50EWOII_50EWOIIIHB_50EWNIIHA_50MSTAR_AVGAGE_AVGSFRAGE_AVGTAU_AVGZMETAL_AVGAV_AVGMU_AVGOIIIHB_AVGSFR_AVGSFR100_AVGB100_AVGB1000_AVGEWOII_AVGEWOIIIHB_AVGEWNIIHA_AVGMSTAR_ERRAGE_ERRSFRAGE_ERRTAU_ERRZMETAL_ERRAV_ERRMU_ERROIIIHB_ERRSFR_ERRSFR100_ERRB100_ERRB1000_ERREWOII_ERREWOIIIHB_ERREWNIIHA_ERRCHI2_2FLAM_1500CFLUX_3727KCORRECT [3]ABSMAG [3]IVARABSMAG [3]SYNTH_ABSMAG [3]ABSMAG_FILTERLIST [3]
int32float64float64float32float32float32float32int16int16int32int32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32int16float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32bytes11
0145.9106933214981-0.66993890371635560.369256685.069725e-09 .. 7.102788e-082.0459266e+20 .. 3.567438e+175.123763e-09 .. 7.3073075e-08089119.932455337586130000.0987901400.011.3249675.37183674.27835460.548214140.026790530.1007801150.26123187-1.052.577934-0.47493076-1.4258919-3.224153-2.8096416-1.0-1.0-1.00-1.0-1.0-1.0-1.011.3253425.37183674.27835460.42949420.026790530.1007801150.26123187-1.0-0.47455597-1.4255171-3.224153-2.8096416-1.0-1.0-1.011.3239915.2335744.37639670.429618030.0235510580.25241980.3406924-1.0-1.4260209-2.2226713-5.2463737-4.7529306-1.0-1.0-1.00.0439102650.62998180.6226340.149637650.00282235680.22359980.171691690.00.88276340.660506253.08060313.16537140.00.00.09.9324558.01557e-183.4443556e-17-0.25023428 .. 0.10797201-22.282818 .. -23.52427935804.156 .. 41486.633-22.2573 .. -23.544544decam_g.par .. decam_z.par
1145.26691496916993-0.9573229077668040.405937673.6528445e-09 .. 5.560954e-083.4718932e+20 .. 3.883426e+173.854969e-09 .. 5.622779e-080622117.9159555456607100000.01268900500.011.44453057.9281776.68880460.619828940.016427250.325471730.19406769-1.052.102215-1.5808468-2.6522436-4.412614-4.0516644-1.0-1.0-1.00-1.0-1.0-1.0-1.011.4479297.2391016.68880460.244609860.0163039290.275765270.2542269-1.0-3.3289576-3.3289576-11.289713-10.224919-1.0-1.0-1.011.4460547.355166.60618160.37456790.01618810.2930070.2579056-1.0-2.5812752-3.0376954-9.357711-8.683728-1.0-1.0-1.00.0137176510.60130620.309186820.147061840.00075990380.072658960.072231670.00.505159140.242822172.7059142.80077460.00.00.07.91595551.4474227e-182.9058928e-17-0.19238524 .. 0.09691145-22.286253 .. -23.5272441834.586 .. 52586.44-22.2679 .. -23.537163decam_g.par .. decam_z.par

Visualize the sample in the rest-frame and the observed frame.


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


Now let's look at observed color-color space.


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


Choose N templates with probability equal to the statistical weight.

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


Selected 3000 templates from 48817 galaxies.

Visualize the results.


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)

Symbol size is proportional to "responsibility" (nperbin) and color is proportional to rest-frame g-r color.


In [22]:
qa_templates()