Updated LRG Templates

ToDo: Update to DR5

The goal of this notebook is to use observations of a large number of luminous red galaxies (LRGs) with DECaLS/DR3 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 2,511,248 LRGs selected from the official DECaLS/DR3.1 target catalog targets-dr3.1-0.14.0.fits using the selection criterion

DESI_TARGET & desi_mask.LRG != 0

Subsequently, Rongpu Zhou (Pittsburgh) assembled spectroscopic redshifts for the sample from SDSS/Main, SDSS/BOSS, AGES, DEEP2, and VIPERS, supplemented with COSMOS medium-band photometric redshifts, yielding precise redshifts for 62,177 objects (2.4% of the parent sample). Using these data, Zhou used a random forest technique (described elsewhere) to compute photometric redshifts for the complete sample, producing the catalog dr3.1-0.14.0-lrg-rf-photoz-0.2.fits, which we use as input into the generation of the spectral templates.

Specifically, we generate templates using the following procedure:

  1. We focus exclusively on the subset of 62,177 LRGs with spectroscopic redshifts since we find that the training sample is fairly representative of the full sample in the four-dimensional space of g-r, r-z, z-band magnitude, and redshift. (Future versions of these templates may utilize the random forest photometric redshifts.)
  2. 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 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.
  3. Next, the notebook below documents how we select a "representative" subset of these galaxies to serve as our LRG template set.
  4. 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 = 'v1.0'
templates_version = 'v2.0'

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."""
    
    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(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['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()


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

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


dict_keys(['Mz', 'redshift', 'rzobs', 'Mg', 'Mr', 'zW1obs', 'gr', 'Mstar', 'rz', 'grobs'])

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]
int32float64float64float32float32float32float32int16int16int32int32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32int16float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32str11
07.10186257679-3.841642608270.4933051.7658e-09 .. 2.19342e-082.10202e+20 .. 1.63934e+171.95907e-09 .. 2.72933e-0803005113.40343.79652e+112.49584e+0911.35728.433798.250550.09161570.007751650.2084180.506265-1.052.1077-3.42061-3.42061-15.0-15.0-1.0-1.0-1.00-1.0-1.0-1.0-1.011.33047.997436.936510.1240430.008730570.3142160.357754-1.0-3.42019-3.42019-15.0-13.4455-1.0-1.0-1.011.3197.492196.911420.2906270.00894760.366890.349315-1.0-2.85242-3.17479-10.7666-10.4481-1.0-1.0-1.00.0592420.8562110.864430.1995980.002190140.3238970.09187060.00.7350660.4933062.934093.001250.00.00.03.40348.62409e-191.93438e-170.0124708 .. 0.0714089-22.1503 .. -23.29988778.52 .. 9104.59-22.1392 .. -23.3014decam_g.par .. decam_z.par
17.23270630501-3.739003233830.5961241.21947e-09 .. 3.73271e-083.41204e+20 .. 1.61952e+171.38542e-09 .. 4.38815e-0801746114.812937.38799e+113.82096e+0911.65276.852196.734490.05884570.01442480.1744840.521158-1.052.1897-3.13147-3.13147-15.0-15.0-1.0-1.0-1.00-1.0-1.0-1.0-1.011.64976.919776.728350.1009160.01442480.4067310.254227-1.0-3.12706-3.13482-15.0-15.0-1.0-1.0-1.011.65237.065516.693760.1858810.01410040.3726660.320061-1.0-2.96507-3.12227-12.9086-12.6098-1.0-1.0-1.00.02305530.3532240.3654860.1137960.001430680.1888520.07775560.00.2490840.0244882.430042.541530.00.00.04.812936.44898e-191.93495e-170.415429 .. 0.191121-22.7973 .. -24.05729283.86 .. 18576.3-22.7938 .. -24.0487decam_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_corner():
    fig = corner(np.vstack((cat['redshift'], cat['Mr'], cat['gr'], cat['rz'])).T, 
                 labels=('Redshift', r'$M_{0.0r}$', r'$^{0.0}(g - r)$', r'$^{0.0}(r - z)$'),
                 range=(zlim, Mrlim, grlim, rzlim))

Corner plot showing the inter-relationship of various rest-frame quantities.


In [14]:
qa_rest_corner()



In [15]:
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', 
                    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', 
               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', 
               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', 
                    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}$ (Number of Galaxies)')

    for aa in (ax1, ax2, ax3, ax4):
        aa.grid(True)
    
    plt.subplots_adjust(wspace=0.2, right=0.85)

Another version of the previous plot but focusing on particular relationships.


In [16]:
qa_rest()


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


In [17]:
grobslim, rzobslim = (0.5, 3), (0.8, 2.2)

In [18]:
def qa_obs():
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
    ax1.hexbin(cat['rzobs'], cat['grobs'], mincnt=1, bins='log', 
               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', 
                    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}$')
    
    for aa in (ax1, ax2):
        aa.grid(True)

In [19]:
qa_obs()


Group similar templates.

Here's the black magic. 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 opt for a much simpler procedure.

Based on the preceding plots, we have opted to use the rest-frame g-r and r-z color together with the r-band absolute magnitude to identify "like" galaxies. Specifically, we simply choose a random subset of NPERBIN galaxies in each 3D bin (containing at least MINCOUNT objects) of rest-frame color-color-magnitude space. This procedure has the benefit of being relatively simple and yielding diversity of templates.

In addition, the NPERBIN and MINCOUNT parameters (as well as the size of the bins) can be tuned to yield the desired number of total templates. Here, we opt for approximately ~1100 templates.


In [20]:
def group_templates(X_in, lim_in, delta_in):
    """Use scipy's binned_statistic_dd method to assigned galaxies to 
    fixed input bins.    
    """
    from scipy.stats import binned_statistic_dd
    
    nbins = [np.ceil(np.ptp(lim) / delta).astype('int') for lim, delta in zip(lim_in, delta_in)]

    X = np.vstack( X_in ).T

    stat, edges, num = binned_statistic_dd(X, _, statistic='count', bins=nbins, range=lim_in)
    
    index = [np.where(num == ii)[0] for ii in np.unique(num)]
    count = np.array([len(aa) for aa in index])
    
    return index, count

The unused code in the next cell demonstrates how, e.g., redshift can be used in the template selection.


In [21]:
if False:
    X_in = (cat['redshift'], cat['gr'], cat['rz'])
    lim_in = (zlim, grlim, rzlim)
    delta_in = (0.05, 0.05, 0.05)
else:
    X_in = (cat['Mr'], cat['gr'], cat['rz'])
    lim_in = (np.flip(Mrlim, 0), grlim, rzlim)
    delta_in = (0.1, 0.05, 0.05)

In [22]:
%time index, count = group_templates(X_in, lim_in, delta_in)


CPU times: user 201 ms, sys: 10.2 ms, total: 212 ms
Wall time: 232 ms

In [23]:
print('Identified {} color-color bins with at least one galaxy:'.format(len(index)))
for mnmx in ((1, 5), (5, 10), (10, 30), (30, 100), (100, count.max()+1)):
    cnt = np.sum( (count >= mnmx[0]) * (count < mnmx[1]) ).astype('int')
    print('  {:03} bins with {} - {}'.format(cnt, mnmx[0], mnmx[1]))


Identified 924 color-color bins with at least one galaxy:
  457 bins with 1 - 5
  111 bins with 5 - 10
  141 bins with 10 - 30
  110 bins with 30 - 100
  105 bins with 100 - 1194

In [24]:
def choose_templates(mincount=3, nperbin=2):
    """Return the indices of the chosen templates and the responsibility, 
    or the number of galaxies in that template's bin.
    
    """
    if nperbin > mincount:
        print('NPERBIN should be less than or equal to MINCOUNT.')
        
    iarch = list()
    norm = list()
    good = np.where(count >= mincount)[0]
    for gg in good:
        iarch.append(rand.choice(index[gg], nperbin, replace=False))
        
    iarch = np.hstack(iarch)
    resp = count[good]
    
    print('Chose {} templates to represent {} galaxies.'.format(len(iarch), np.sum(resp)))
    
    return iarch, resp

In [25]:
iarch, resp = choose_templates()


Chose 1142 templates to represent 48294 galaxies.

In [26]:
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 [27]:
def qa_templates():
    """Generate a color-color plot with the symbol size scaled by the responsibility.
    
    """
    
    size, col = _markers()

    fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 10))

    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)
    
    for aa in (ax1, ax2, ax3, ax4):
        aa.grid(True)    
    
    fig.subplots_adjust(wspace=0.3)

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


In [28]:
qa_templates()


Some additional diagnostics of the 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 [29]:
def qa_template_properties():
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5), sharey=True)
    
    ax1.hist(ised['MSTAR_50'], normed=True, bins=30, label='Full Sample')
    ax1.hist(ised['MSTAR_50'][iarch], normed=True, bins=30, alpha=0.5,
            label='Templates')
    ax1.set_xlabel(r'$\log_{10}\ (M\ /\ M_{\odot})$')
    ax1.set_ylabel('Relative Number')
    
    ax2.hist(ised['SFR100_50'], normed=True, bins=30, label='Full Sample')
    ax2.hist(ised['SFR100_50'][iarch], normed=True, bins=30, alpha=0.5,
            label='Templates')    
    ax2.set_xlabel(r'$\log_{10}\ (SFR_{100}\ /\ M_{\odot}\ {\rm yr}^{-1})$')
    ax2.legend(loc='upper right')
    
    plt.subplots_adjust(wspace=0.02)

In [30]:
qa_template_properties()


Write out the metadata.


In [31]:
print('Writing {}'.format(metadata_outfile))
ised[iarch].write(metadata_outfile, overwrite=True)


Writing /Users/ioannis/work/desi/spectro/templates/lrg_templates/isedfit/v1.0/lrg-templates-isedfit-v1.0.fits

Visualize the resulting templates.

The spectra have to built in IDL using the script build_lrg_basis_templates, which we assume has happend.


In [32]:
from desisim.io import read_basis_templates

In [33]:
flux, wave, meta = read_basis_templates('LRG')
nt = len(meta)
print('Number of templates = {}'.format(nt))


INFO:io.py:967:read_basis_templates: Reading /Users/ioannis/work/desi/spectro/templates/basis_templates/v2.4/lrg_templates_v2.0.fits
Number of templates = 1142

In [34]:
meta


Out[34]:
<Table length=1142>
TEMPLATEIDISEDFIT_IDRADECZMAGGIES [5]IVARMAGGIES [5]LOGMSTARLOGSFRLOGSFR100AV_ISMAGESFRAGEHBETA_CONTINUUMD4000
int32int32float64float64float32float32float32float32float32float32float32float32float32float32float32
038770161.6607979221.973973676190.717822.08524e-09 .. 8.69486e-081.98737e+20 .. 1.28632e+1712.0261-1.66163-2.249980.1280756.010995.148971.99862e-111.66315
12372939.20931880941.593000227740.7176512.0286e-09 .. 8.40202e-082.14052e+20 .. 2.16765e+1711.9925-2.72429-2.767250.002455025.816575.732882.07025e-111.68762
251254320.0643084941.059139902420.7186371.92813e-09 .. 9.23302e-087.5455e+19 .. 1.16828e+1711.9723-0.873839-1.542650.0548856.24944.951562.04717e-111.70254
340678189.85453137516.86523256330.6880993.12021e-09 .. 9.45608e-081.001e+19 .. 2.975e+1612.0011-0.562495-1.245660.2495447.138265.446551.8156e-111.60199
42060939.5509326937-1.047029364380.7504911.19407e-09 .. 7.81926e-081.78903e+20 .. 1.65252e+1711.9544-1.67601-2.353850.06098445.864375.029412.06431e-111.6996
532908.535977446772.365875901130.7113371.74082e-09 .. 7.51581e-083.75509e+20 .. 9.77913e+1612.0359-2.44352-2.674180.131936.991696.954641.63681e-111.70158
645556233.4678909192.631263716210.654193.23746e-09 .. 8.56976e-081.6822e+20 .. 1.1953e+1711.8795-0.11013-0.870150.2282516.536154.904062.10307e-111.53727
756606.836837250889.887478923450.6592862.74467e-09 .. 1.16089e-073.76665e+19 .. 1.01018e+1711.8445-0.551618-1.338460.02803582.734572.263763.34019e-111.68822
81821534.0672896261-0.574316869950.6696211.9005e-09 .. 6.34756e-086.16894e+19 .. 5.56466e+1611.9387-2.45091-2.683140.131936.991696.954641.63682e-111.70158
91128015.20874245971.366696320670.8315716.4923e-10 .. 6.3998e-083.3673e+20 .. 1.29444e+1711.8964-1.51695-2.293930.06098445.864375.029412.06431e-111.6996
.............................................
11321386426.1353388741-0.28184161480.3857661.46689e-09 .. 4.78138e-083.04167e+20 .. 2.03254e+1711.129-1.49694-1.924151.357674.743363.316658.95616e-121.49408
113336046150.5353590511.473235789660.39091.23679e-09 .. 3.12848e-084.06701e+20 .. 1.68091e+1711.2664-2.20425-2.622221.603226.429356.390284.65853e-121.61182
11342396133.64979495172.660692058620.3477881.8939e-09 .. 3.97278e-081.46867e+20 .. 1.83394e+1711.3004-3.20817-3.403580.9636128.874958.304645.50758e-121.89368
113541983.341939802972.228939501250.3448411.92395e-09 .. 3.78447e-083.41046e+20 .. 1.95181e+1711.3949-3.30314-3.329071.603226.429356.390284.65853e-121.61183
113641872219.840478089-0.96662024890.3223392.57094e-09 .. 2.63113e-086.14136e+19 .. 2.64669e+1711.0722-3.38744-3.588150.09489129.041168.856911.11887e-111.97237
113739511160.91202105824.08682970950.3682831.50836e-09 .. 2.50217e-084.11457e+20 .. 2.57469e+1711.0444-2.73832-3.190280.2495836.350786.240021.1681e-111.95586
113835079149.23279370417.69180121760.3680611.46792e-09 .. 2.61672e-082.44238e+20 .. 1.96902e+1711.2281-3.00041-3.307760.7809866.794276.476137.29186e-121.897
113941678.155856320151.223477397750.3578291.72951e-09 .. 3.82299e-085.1636e+20 .. 1.7327e+1711.1867-0.953373-1.634280.6966868.427566.196758.44085e-121.75777
114036777150.1757478952.686414401570.44986.68522e-10 .. 1.98241e-084.39116e+20 .. 2.15053e+1711.2452-2.7689-3.246440.7496838.022016.884996.63798e-122.0129
11411738827.48703500560.2597459434570.3570061.63789e-09 .. 3.61699e-085.67175e+20 .. 2.41533e+1711.2612-3.2071-3.446110.9636128.874958.304645.50758e-121.89369

In [35]:
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 [36]:
plot_subset()


Finally, compare the colors of the templates to the data.

This test brings this analysis full circle.


In [37]:
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.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(0.0, 2.0, 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 [38]:
%time cc = zgrid_colors()


Number of redshift points = 10
CPU times: user 3.93 s, sys: 2.45 s, total: 6.38 s
Wall time: 8.07 s

In [39]:
grrange = (0.0, 3.0)
rzrange = (0.0, 2.5)
zW1range = (-1, 4)
ntspace = 20 # spacing between model curves

In [40]:
def grz(pngfile=None):
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4), sharey=True)

    hb = ax1.hexbin(cat['rzobs'], cat['grobs'], bins='log', 
                   mincnt=1, extent=rzrange+grrange,
                   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)

    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.05, 'z=0', ha='left', va='bottom',
             transform=ax2.transAxes, fontsize=14)
    ax2.text(0.05, 0.9, 'Models (z=0-1, 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()
    
    plt.subplots_adjust(wspace=0.15)
    
    if pngfile:
        fig.savefig(pngfile)

In [41]:
grz()



In [42]:
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,
                   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)

    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.01, 0.05, 'z=0', ha='left', va='bottom',
             transform=ax2.transAxes, fontsize=14)
    ax2.text(0.05, 0.9, 'Models (z=0-1, 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()
    
    plt.subplots_adjust(wspace=0.15)
    
    if pngfile:
        fig.savefig(pngfile)

In [43]:
rzW1()