Gallery for DR8

The purpose of this notebook is to build a nice gallery of object images for DR8. The theme is intermediate-redshift galaxy clusters.

For future reference: The notebook must be run from https://jupyter.nersc.gov with the following (approximate) activation script:

#!/bin/bash                                                                                                           
version=$1                                                                                                            
connection_file=$2                                                                                                    

desiconda_version=20180512-1.2.5-img
module use /global/common/software/desi/$NERSC_HOST/desiconda/$desiconda_version/modulefiles
module load desiconda

export LEGACYPIPE_DIR=$HOME/repos/git/legacypipe                                                                       

export PATH=$LEGACYPIPE_DIR/bin:${PATH}                                                                               
export PYTHONPATH=$LEGACYPIPE_DIR/py:${PYTHONPATH}                                                                    

export DUST_DIR=/global/project/projectdirs/desi/software/${NERSC_HOST}/dust/v0_0

export UNWISE_COADDS_DIR=/global/projecta/projectdirs/cosmo/work/wise/outputs/merge/neo4/fulldepth:/global/project/projectdirs/cosmo/data/unwise/allwise/unwise-coadds/fulldepth
export UNWISE_COADDS_TIMERESOLVED_DIR=/global/projecta/projectdirs/cosmo/work/wise/outputs/merge/neo4

export GAIA_CAT_DIR=/global/project/projectdirs/cosmo/work/gaia/chunks-gaia-dr2-astrom/
export GAIA_CAT_VER=2

export TYCHO2_KD_DIR=/global/project/projectdirs/cosmo/staging/tycho2

export PS1CAT_DIR=/global/project/projectdirs/cosmo/work/ps1/cats/chunks-qz-star-v3 # calibration only

exec python -m ipykernel -f $connection_file

Imports and paths


In [1]:
import os, sys
import shutil, time, warnings
from contextlib import redirect_stdout
import numpy as np
import numpy.ma as ma
import matplotlib.pyplot as plt

In [2]:
import fitsio
from astropy.io import ascii
import astropy.units as u
from astropy.coordinates import SkyCoord, FK5
from astropy.table import Table, vstack, Column
from astropy.cosmology import FlatLambdaCDM
from PIL import Image, ImageDraw, ImageFont

In [3]:
from legacypipe.survey import LegacySurveyData
from legacypipe.runbrick import run_brick
from astrometry.util.util import Tan
from astrometry.util.fits import merge_tables
from astrometry.libkd.spherematch import match_radec

In [4]:
import multiprocessing
nproc = multiprocessing.cpu_count() // 2

In [5]:
cosmo = FlatLambdaCDM(H0=70, Om0=0.3)

In [6]:
%matplotlib inline

Preliminaries

Define the data release and the various output directories.


In [7]:
PIXSCALE = 0.262
#PIXSCALE = 2.62

In [8]:
dr = 'dr8'
drdir = '/global/project/projectdirs/cosmo/work/legacysurvey/dr8'
gallerydir = os.path.join( os.getenv('SCRATCH'), dr, 'gallery')
samplefile = os.path.join(gallerydir, 'sample-{}.fits'.format(dr))

In [9]:
jpgdir = os.path.join(gallerydir, 'jpg')
if not os.path.isdir(jpgdir):
    os.mkdir(jpgdir)

In [10]:
pngdir = os.path.join(gallerydir, 'png')
if not os.path.isdir(pngdir):
    os.mkdir(pngdir)

Build the sample.

Use the sample of clusters from http://risa.stanford.edu/redmapper/redmapper_dr8_public_v6.3_catalog.fits.gz

Cut to lambda>100, corresponding to M(halo) = (0.1-3) x 10^15 Msun over the redshift range 0.1-0.5.


In [11]:
def read_redmapper(radius_kpc=200, seed=1, npilot=100, 
                   richmin=20, richmax=300, zmin=0.1, 
                   zmax=0.5):
    rmfile = os.path.join(gallerydir, 'redmapper_dr8_public_v6.3_catalog.fits.gz')
    cat = Table(fitsio.read(rmfile, lower=True))
    cut = ((cat['lambda'] > richmin) * (cat['z_lambda'] > zmin) * 
           (cat['z_lambda'] < zmax) * (cat['p_cen'][:, 0] > 0.9))
    cat = cat[cut]
    
    # Choose 100 clusters drawn uniformly in richness
    #redshift = cat['z_lambda']
    richness = cat['lambda']
    nbin = 20
    _xbin = np.linspace(richmin, richmax, nbin)
    idx  = np.digitize(richness, _xbin)

    prob = np.zeros_like(richness)
    for kk in range(nbin):
        ww = idx == kk
        if np.sum(ww) > 1:
            prob[ww] = 1 / np.sum(ww)
    prob /= np.sum(prob)

    rand = np.random.RandomState(seed=seed)
    these = rand.choice(len(cat), npilot, p=prob, replace=False)
    srt = np.argsort(cat['lambda'][these])[::-1]
    cat = cat[these[srt]]
    
    # Get the diameter in pixels corresponding to 500 kpc (comoving).
    arcsec_per_kpc = cosmo.arcsec_per_kpc_proper(cat['z_lambda']).value
    radius_arcsec = radius_kpc * arcsec_per_kpc # [float arcsec]
    diam = np.rint(2 * radius_arcsec / PIXSCALE).astype(int) # [integer/rounded pixels]
    
    cat.add_column(Column(name='diam', data=diam, dtype='i4')) # [pixels]
    cat.add_column(Column(name='run', length=len(cat), dtype='S14'))
    
    return cat

In [12]:
cat = read_redmapper()
cat


Out[12]:
<Table length=100>
idnameradecz_lambdaz_lambda_errlambdalambda_errsz_specobjidimagimag_errmodel_mag_umodel_magerr_umodel_mag_gmodel_magerr_gmodel_mag_rmodel_magerr_rmodel_mag_imodel_magerr_imodel_mag_zmodel_magerr_zilump_cen [5]ra_cen [5]dec_cen [5]id_cen [5]pzbins [21]pz [21]diamrun
int32bytes20float64float64float32float32float32float32float32float32int64float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32int64float32float32int32bytes14
6RMJ133520.1+410004.1203.8337226793197541.00114644090520.231746940.0060674027189.181155.6133591.01003770.228308712376623067224474980.00.019.9958040.10798485617.9016720.00915489216.3668250.004500995415.8506550.004327438315.5079190.011259914128.978870.98141843 .. 3.6546284e-05203.83372 .. 203.8387641.00115 .. 40.99122625032419 .. 250324230.1959384 .. 0.267555458.0235895e-07 .. 5.993135e-06413
27RMJ094951.8+170710.6147.465841075984717.1195974975727370.383245080.013010344185.382667.11442661.02126370.3827856212376677339563953410.00.022.2165850.983209419.6990780.05694011217.8786470.01596755317.1836010.01232830116.915440.031415425126.3268360.99873763 .. 1.4405185e-06147.46584 .. 147.4611817.119598 .. 17.09875339424958 .. 400218260.31329086 .. 0.453199277.051299e-05 .. 4.096737e-06292
215RMJ111514.8+531954.6168.811865624597453.331823595559290.47128370.010392168184.8182214.1453151.59291710.466421212376575897823478910.00.023.0007821.400842520.5103110.06005492418.5997330.01843009517.7054230.01315209117.2800270.031151984106.089760.99795884 .. 7.910885e-06168.81186 .. 168.8166753.331825 .. 53.32002610643753 .. 106438700.41006044 .. 0.53250696.312513e-07 .. 5.772295e-06258
250RMJ222842.7+083924.4337.177784353191448.6567820364482270.417336460.01182232181.1265310.9680411.4284549-1.012376790093157503380.00.020.965210.4600936819.8763430.063508617.89530.01756582217.1981140.014344024516.7907140.03765561115.92240.9999326 .. 4.5067432e-07337.1778 .. 337.17978.656782 .. 8.655086552434029 .. 524340300.3519025 .. 0.482770472.1631968e-05 .. 3.0974295e-06277
5RMJ090912.2+105824.9137.3007446351910510.973594935457740.170462060.0036083707174.704184.94786361.0145435-1.012376712601265769200.00.026.3431223.217321618.1581670.02315970117.1842330.01367530516.5320340.01074851716.2881340.028455446112.058440.94703996 .. 0.00021698933137.30075 .. 137.3293910.973595 .. 10.98569245540350 .. 455403200.15001194 .. 0.190912191.0995842e-05 .. 1.323996e-05526
129RMJ015949.3-084958.929.95555555246912-8.8330347640603540.408139440.012372703167.778438.6551521.19309090.405218212376709567906449410.00.019.9394910.136068819.2299390.02570270417.8157560.01234654117.1803970.0104227616.8070580.028954698102.938860.9989383 .. 3.9697662e-0729.955555 .. 29.944187-8.8330345 .. -8.84102644913646 .. 449137140.34273946 .. 0.473539380.00018539147 .. 5.0218814e-06281
424RMJ110608.5+333339.7166.535353417260633.561023713485920.497296630.011500471167.2444915.7361071.79648140.4881110812376650163130210440.00.026.9847130.585276120.9724750.0837710619.2539120.02810034918.333690.01919185217.7582650.04585198799.7359540.99532783 .. 0.00020922536166.53535 .. 166.5051933.561024 .. 33.5605231097894 .. 310978540.42666465 .. 0.56792861.8874722e-08 .. 4.8688376e-06251
320RMJ221145.9-034944.5332.94131100281584-3.82902869947158830.426461220.011818006159.418710.3009391.481744-1.012376800664161486490.00.021.5699841.34560618.9634590.04026359717.1440750.012984971516.4158630.00961958616.1468010.029088078104.401420.9960301 .. 2.8823106e-05332.9413 .. 332.93222-3.8290286 .. -3.823425854896963 .. 548970640.3610203 .. 0.491902082.1219097e-05 .. 4.711655e-06273
304RMJ012542.3-063442.321.426365942207184-6.578423751845430.435859080.011315561159.3424210.55244351.4463356-1.012376802430703375700.00.021.2910210.654548120.1512430.0767922718.367980.02607430717.5436170.02065530817.1355110.04723501691.02840.9108611 .. 7.153804e-0621.426367 .. 21.42733-6.578424 .. -6.61754955986054 .. 559860520.3731506 .. 0.498567582.9843713e-05 .. 4.797557e-06270
56RMJ121218.5+273255.1183.077010105166827.5486477313799720.353901150.013516471158.244066.03478051.00771140.3506435812376673232562423920.00.022.751721.500492119.1608260.0272529517.390930.01057453516.815860.00948453216.4627380.02650024104.596090.9128537 .. 2.2766733e-05183.07701 .. 183.0783127.548647 .. 27.54095638037007 .. 380370080.27750033 .. 0.430301961.5178048e-06 .. 2.9256653e-06307
................................................................................................
18748RMJ124409.3+013923.7191.03862988376261.65658562264940760.37389720.01585863230.8040073.13947370.999533950.4002188712376712665711250760.00.021.8563580.358693919.904570.03249834518.152530.01148736317.5292170.01071552517.129330.02227587822.396190.95947415 .. 3.0555446e-05191.03864 .. 191.048871.6565856 .. 1.623865845852060 .. 458520200.30031574 .. 0.447478650.0030597614 .. 3.0352016e-06296
17773RMJ093400.1+081821.8143.50033095178788.306051735310790.360870870.01646095330.3311023.1498961.0558510.378248712376604131892925400.00.025.1203253.675136319.7285180.04661588717.8934760.01652014417.3146380.01653361116.9289530.03906598722.713930.9999755 .. 2.5268622e-07143.50034 .. 143.507548.306052 .. 8.29020516263775 .. 162637810.27024 .. 0.451501732.0469513e-06 .. 2.5075133e-06303
28428RMJ001256.8+353148.13.236865514470764535.530018402732150.364454660.01797888428.5310153.40497351.0642918-1.012376805093382560880.00.024.7011452.603279820.728820.0664565618.8801730.0173699718.2223990.01476502517.8842720.03661260818.7283330.992062 .. 5.0772986e-053.2368655 .. 3.240557735.530018 .. 35.51899758247281 .. 582472830.27620688 .. 0.45270240.0002079113 .. 2.3311693e-06301
28698RMJ090909.2+365124.2137.2881690763986236.85671956017370.338520740.01736423625.388492.89623260.998094260.3669912812376603444704052670.00.021.2977940.2548442220.0061990.0370368918.2991260.01497882417.6936150.01457523517.3083760.03321632717.2992730.91638786 .. 2.83402e-07137.28816 .. 137.3151236.85672 .. 36.8096116169608 .. 161696820.22554575 .. 0.45149571.561465e-14 .. 2.7524254e-06316
22830RMJ231826.3+070312.7349.60972376208267.0535369125356750.227792350.00861421823.9850332.9634941.0208337-1.012376696806465210710.00.020.4979840.3324347718.5106620.02068412917.0912880.00921956616.587610.0084255516.2610470.0223798615.8941490.9548585 .. 4.026261e-06349.6097 .. 349.573067.053537 .. 7.02746744068566 .. 440684960.16922867 .. 0.286356062.617821e-14 .. 5.558568e-06419
33706RMJ012251.5-061659.520.714637099409146-6.2832038790955980.33682170.01792341122.962952.9125490.9968275-1.012376802436068803960.00.024.4924583.758366319.9082360.0333163618.279820.01353948617.6775570.01230058117.3109990.03095805316.7867850.99892 .. 0.0001128190420.714638 .. 20.701805-6.283204 .. -6.264387656039440 .. 560395550.23194489 .. 0.44169852.0127823e-12 .. 2.6232542e-06317
31226RMJ150357.4+362756.2225.9893747602662336.465601912209230.336759780.01674640222.383262.61476830.991402570.369834212376618504096320320.00.021.8472750.509415420.002910.04182582718.3506930.01565876817.7407040.01372300217.3616580.03293065417.8860130.92852014 .. 1.6650607e-05225.98938 .. 225.9964636.465603 .. 36.46301320347202 .. 203472630.23341492 .. 0.440104662.1594304e-12 .. 2.525342e-06317
24247RMJ004529.0+081326.511.3709570009083348.224041496854080.274721120.01303967322.0476512.49341730.992919150.2570715512376697031963445060.00.020.3739680.2234056118.6833520.02049750517.2672820.00930830116.7514060.00951492716.421190.02314286116.1527060.99796563 .. 9.754984e-0911.370957 .. 11.3773798.224042 .. 8.20125944362247 .. 443623020.18487516 .. 0.364567072.0349958e-14 .. 3.56539e-06365
20509RMJ083819.5+181435.8129.5812525329628818.2432688958595040.212595660.007575539421.4374522.24885461.019020.2171032612376672115806703750.00.021.5093350.3121999519.1742820.01789001617.7937720.00952519617.2764190.0091533916.9151860.0247817913.1966120.9201768 .. 4.4588403e-05129.58125 .. 129.5793318.243269 .. 18.29815736955171 .. 369551750.1604335 .. 0.26475787.3292494e-14 .. 6.63762e-06441
46629RMJ005812.4+095401.414.5518354363254059.900377125603810.327885630.01846754720.9492932.76763320.990400850.2934065212376789067888070790.00.022.792361.341735420.3546030.04865289518.7549780.02076700918.1256140.01660628817.794440.04686328814.0378950.9686195 .. 0.001240854614.551835 .. 14.5506459.900377 .. 9.89393151930253 .. 519302510.20197922 .. 0.453792049.523486e-20 .. 2.738783e-06323

In [13]:
#_ = plt.hist(cat['mag_10'], bins=30)
_ = plt.hist(cat['z_lambda'], bins=30)



In [14]:
_ = plt.hist(cat['lambda'], bins=30)



In [15]:
_ = plt.hist(cat['diam'], bins=30)


Restrict to objects that are in the DR8 footprint.


In [16]:
def read_all_ccds(dr='dr8'):                                                                                                        
    """Read the CCDs files, treating DECaLS and BASS+MzLS separately.                                                               
                                                                                                                                    
    """                                                                                                                             
    from astrometry.libkd.spherematch import tree_open                                                                              
    #survey = LegacySurveyData()                                                                                                    
                                                                                                                                    
    kdccds_north = []                                                                                                               
    for camera in ('90prime', 'mosaic'):                                                                                            
        ccdsfile = os.path.join(drdir, 'survey-ccds-{}-{}.kd.fits'.format(camera, dr))                                              
        ccds = tree_open(ccdsfile, 'ccds')                                                                                          
        print('Read {} CCDs from {}'.format(ccds.n, ccdsfile))                                                                      
        kdccds_north.append((ccdsfile, ccds))                                                                                       
                                                                                                                                    
    ccdsfile = os.path.join(drdir, 'survey-ccds-decam-{}.kd.fits'.format(dr))                                                       
    ccds = tree_open(ccdsfile, 'ccds')                                                                                              
    print('Read {} CCDs from {}'.format(ccds.n, ccdsfile))                                                                          
    kdccds_south = (ccdsfile, ccds)                                                                                                 
                                                                                                                                    
    return kdccds_north, kdccds_south

In [17]:
kdccds_north, kdccds_south = read_all_ccds()


Read 132066 CCDs from /global/project/projectdirs/cosmo/work/legacysurvey/dr8/survey-ccds-90prime-dr8.kd.fits
Read 206834 CCDs from /global/project/projectdirs/cosmo/work/legacysurvey/dr8/survey-ccds-mosaic-dr8.kd.fits
Read 3913907 CCDs from /global/project/projectdirs/cosmo/work/legacysurvey/dr8/survey-ccds-decam-dr8.kd.fits

In [18]:
survey = LegacySurveyData()
survey.output_dir = gallerydir

In [19]:
def get_name(cat, nice=False):
    name = np.atleast_1d(cat['name'])
    if nice:
        outname = name
    else:
        outname = np.array([nn.replace(' ', '-').lower() for nn in name])
    if len(outname) == 1:
        outname = outname[0]
    return outname

In [20]:
def _build_sample_one(args):
    """Wrapper function for the multiprocessing."""
    return build_sample_one(*args)

In [21]:
def build_sample_one(cat, verbose=False):
    """Determine the "run", i.e., determine whether we should use the BASS+MzLS CCDs
    or the DECaLS CCDs file when running the pipeline. 

    """
    from astrometry.util.fits import fits_table, merge_tables   
    from astrometry.util.util import Tan                                    
    from astrometry.libkd.spherematch import tree_search_radec  
    from legacypipe.survey import ccds_touching_wcs                         

    ra, dec, diam = cat['ra'], cat['dec'], cat['diam']
    if dec < 25:                        
        run = 'decam'                       
    elif dec > 40:                                      
        run = '90prime-mosaic'      
    else:                               
        width = np.ceil(diam * 60 / PIXSCALE).astype(int)
        wcs = Tan(ra, dec, width/2+0.5, width/2+0.5,
                  -PIXSCALE/3600.0, 0.0, 0.0, PIXSCALE/3600.0, 
                  float(width), float(width))  
                        
        # BASS+MzLS 
        TT = []         
        for fn, kd in kdccds_north:
            I = tree_search_radec(kd, ra, dec, 1.0)
            if len(I) == 0:
                continue
            TT.append(fits_table(fn, rows=I))  
        if len(TT) == 0:                    
            inorth = []                     
        else:                           
            ccds = merge_tables(TT, columns='fillzero')     
            inorth = ccds_touching_wcs(wcs, ccds)   
                                    
        # DECaLS            
        fn, kd = kdccds_south  
        I = tree_search_radec(kd, ra, dec, 1.0) 
        if len(I) > 0:          
            ccds = fits_table(fn, rows=I)  
            isouth = ccds_touching_wcs(wcs, ccds)  
        else: 
            isouth = []
        if len(inorth) > len(isouth): 
            run = '90prime-mosaic' 
        else:
            run = 'decam' 
            
        if verbose:
            print('Cluster {}: RA, Dec={:.6f}, {:.6f}: run={} ({} north CCDs, {} south CCDs).'.format(
                cat['name'], ra, dec, run, len(inorth), len(isouth)))
            
        if (len(inorth) == 0) * (len(isouth) == 0):
            run = None
            
    if run:
        cat['run'] = run
        return cat
    else:
        return None

In [22]:
def build_sample(cat, use_nproc=nproc, verbose=False):
    """Build the full sample with grz coverage in DR7."""

    sampleargs = list()
    for cc in cat:
        sampleargs.append( (cc, verbose) )

    if use_nproc > 1:
        p = multiprocessing.Pool(nproc)
        result = p.map(_build_sample_one, sampleargs)
        p.close()
    else:
        result = list()
        for args in sampleargs:
            result.append(_build_sample_one(args))

    #result = list(zip(*result))
    # Remove non-matching objects and write out the sample
    outcat = vstack(list(filter(None, result)))
    print('Found {}/{} clusters in the DR8 footprint.'.format(len(outcat), len(cat)))
    
    return outcat

In [23]:
#rr = build_sample(cat, use_nproc=4, verbose=True)

In [24]:
samplelogfile = os.path.join(gallerydir, 'build-sample.log')
print('Building the sample.')
print('Logging to {}'.format(samplelogfile))
t0 = time.time()
with open(samplelogfile, 'w') as log:
    with redirect_stdout(log):
        sample = build_sample(cat)
print('Found {}/{} objects in the DR8 footprint.'.format(len(sample), len(cat)))
print('Total time = {:.3f} seconds.'.format(time.time() - t0))


Building the sample.
Logging to /global/cscratch1/sd/ioannis/dr8/gallery/build-sample.log
Found 100/100 objects in the DR8 footprint.
Total time = 1.771 seconds.

In [25]:
print('Writing {}'.format(samplefile))
sample.write(samplefile, overwrite=True)


Writing /global/cscratch1/sd/ioannis/dr8/gallery/sample-dr8.fits

In [26]:
#sample

In [27]:
[print(ss['name'].replace(' ', '-'), ss['ra'], ss['dec']) for ss in sample[:10]]
#[print(ss['ra'], ss['dec']) for ss in sample]


RMJ133520.1+410004.1 203.83372267931975 41.0011464409052
RMJ094951.8+170710.6 147.4658410759847 17.119597497572737
RMJ111514.8+531954.6 168.8118656245974 53.33182359555929
RMJ222842.7+083924.4 337.17778435319144 8.656782036448227
RMJ090912.2+105824.9 137.30074463519105 10.97359493545774
RMJ015949.3-084958.9 29.95555555246912 -8.833034764060354
RMJ110608.5+333339.7 166.5353534172606 33.56102371348592
RMJ221145.9-034944.5 332.94131100281584 -3.8290286994715883
RMJ012542.3-063442.3 21.426365942207184 -6.57842375184543
RMJ121218.5+273255.1 183.0770101051668 27.548647731379972
Out[27]:
[None, None, None, None, None, None, None, None, None, None]

In [28]:
def qa_sample():
    fig, ax = plt.subplots()
    ax.scatter(cat['ra'], cat['dec'], alpha=0.5, s=10, label='Abell Catalog')
    ax.scatter(sample['ra'], sample['dec'], s=20, marker='s', label='Objects in DR8 Footprint')
    ax.set_xlabel('RA')
    ax.set_ylabel('Dec')
    ax.legend(loc='upper right')

In [29]:
qa_sample()


Generate the color mosaics for each object.


In [30]:
def custom_brickname(obj, prefix='custom-'): 
    brickname = 'custom-{:06d}{}{:05d}'.format(
        int(1000*obj['ra']), 'm' if obj['dec'] < 0 else 'p', 
        int(1000*np.abs(obj['dec'])))
    return brickname

In [31]:
#name = get_name(sample)
#[print(ii, nn, d1, dd) for ii, (nn, d1, dd) in enumerate(zip(name, sample['diam'], sample['diam'] * 60 / 0.3))]
#print(sample['diam'].data * 60 / 1)

In [32]:
def custom_coadds_one(obj, scale=PIXSCALE, clobber=False, log=None):
    import subprocess
    #from legacypipe.runbrick import run_brick
    #from astrometry.util.multiproc import multiproc
    #from legacypipe.runbrick import stage_tims, run_brick
    #from legacypipe.coadds import make_coadds

    name = get_name(obj)
    jpgfile = os.path.join(jpgdir, '{}.jpg'.format(name))
    if os.path.isfile(jpgfile) and not clobber:
        print('File {} exists...skipping.'.format(jpgfile))
    else:
        size = obj['diam']
        #size = np.rint(obj['diam'] * 60 / scale).astype('int') # [pixels]
        #if size > 500:
        #    size = 500
        #print('Generating mosaic for {} with width={} pixels.'.format(name, size))
        
        bands = ('g', 'r', 'z')
        rgb_kwargs = dict({'Q': 20, 'm': 0.03})
            
        brickname = custom_brickname(obj, prefix='custom-')
        
        cmd = 'python {legacypipe_dir}/py/legacypipe/runbrick.py '
        cmd += '--radec {ra} {dec} --width {width} --height {width} --pixscale {pixscale} '
        cmd += '--run {run} '
        cmd += '--threads {threads} --outdir {outdir} --no-write '
        #cmd += '--apodize '
        cmd += '--stage image_coadds --skip-calibs '
    
        cmd = cmd.format(legacypipe_dir=os.getenv('LEGACYPIPE_DIR'), 
                         ra=obj['ra'], dec=obj['dec'], run=obj['run'].strip(),
                         width=size, pixscale=scale, threads=nproc, 
                         outdir=gallerydir)
    
        print(cmd, flush=True, file=log)
        err = subprocess.call(cmd.split(), stdout=log, stderr=log)
        if err != 0:
            print('Something we wrong; please check the logfile for {}.'.format(name))
            return 0
        else:
            #with warnings.catch_warnings():
            #    warnings.simplefilter("ignore")
            #    run_brick(None, survey, radec=(obj['ra'], obj['dec']), pixscale=scale, 
            #              width=size, height=size, rgb_kwargs=rgb_kwargs, threads=nproc, 
            #              stages=['image_coadds'], splinesky=True, early_coadds=True, pixPsf=True, 
            #              hybridPsf=True, normalizePsf=True, write_pickles=False, depth_cut=False, 
            #              apodize=True, do_calibs=False, ceres=False)
            #sys.stdout.flush()    
            
            _jpgfile = os.path.join(survey.output_dir, 'coadd', 'cus', brickname, 
                                   'legacysurvey-{}-image.jpg'.format(brickname))
            if os.path.isfile(_jpgfile):
                shutil.copy(_jpgfile, jpgfile)
                shutil.rmtree(os.path.join(survey.output_dir, 'coadd'))
                shutil.rmtree(os.path.join(survey.output_dir, 'metrics'))
            else:
                print('File {} not found; check log!'.format(_jpgfile))

In [33]:
#custom_coadds_one(sample[320], scale=2, clobber=True)

In [34]:
def custom_coadds(sample, clobber=False, log=None):
    for obj in sample:
        custom_coadds_one(obj, clobber=clobber, log=log)

In [35]:
coaddslogfile = os.path.join(gallerydir, 'make-coadds.log')
print('Generating the coadds.')
print('Logging to {}'.format(coaddslogfile))
t0 = time.time()
with open(coaddslogfile, 'w') as log:
    #with redirect_stdout(log):
    custom_coadds(sample, clobber=True, log=log)
print('Total time = {:.3f} minutes.'.format((time.time() - t0) / 60))


Generating the coadds.
Logging to /global/cscratch1/sd/ioannis/dr8/gallery/make-coadds.log
File /global/cscratch1/sd/ioannis/dr8/gallery/coadd/cus/custom-226728p34551/legacysurvey-custom-226728p34551-image.jpg not found; check log!
File /global/cscratch1/sd/ioannis/dr8/gallery/coadd/cus/custom-003236p35530/legacysurvey-custom-003236p35530-image.jpg not found; check log!
Total time = 51.255 minutes.

Add labels and a scale bar.


In [57]:
barlen = np.round(15.0 / PIXSCALE).astype('int')
fonttype = os.path.join(gallerydir, 'Georgia-Italic.ttf')

In [74]:
def _add_labels_one(args):
    """Wrapper function for the multiprocessing."""
    return add_labels_one(*args)

In [75]:
def add_labels_one(obj, verbose=False):
    name = get_name(obj)
    nicename = get_name(obj, nice=True)
    
    jpgfile = os.path.join(jpgdir, '{}.jpg'.format(name))
    pngfile = os.path.join(pngdir, '{}.png'.format(name))
    thumbfile = os.path.join(pngdir, 'thumb-{}.png'.format(name))
    
    if os.path.isfile(jpgfile):
        im = Image.open(jpgfile)
        sz = im.size
        fntsize = np.round(sz[0]/28).astype('int')
        if fntsize < 10:
            fntsize = 10
        width = np.round(sz[0]/175).astype('int')
        font = ImageFont.truetype(fonttype, size=fntsize)
        draw = ImageDraw.Draw(im)
        # Label the object name--
        draw.text((0+fntsize*2, 0+fntsize*2), nicename, font=font)
        # Add a scale bar--
        x0, x1, yy = sz[1]-fntsize*2-barlen, sz[1]-fntsize*2, sz[0]-fntsize*2
        draw.line((x0, yy, x1, yy), fill='white', width=width)
        im.save(pngfile)    
        
        # Generate a thumbnail
        if False:
            cmd = '/usr/bin/convert -thumbnail 300x300 {} {}'.format(pngfile, thumbfile)
            os.system(cmd)
    else:
        print('File {} not found.'.format(jpgfile))

In [76]:
def add_labels(sample):
    labelargs = list()
    for obj in sample:
        labelargs.append((obj, False))

    #if False:
    if nproc > 1:
        p = multiprocessing.Pool(nproc)
        res = p.map(_add_labels_one, labelargs)
        p.close()
    else:
        for args in labelargs:
            res = _add_labels_one(args)

In [77]:
%time add_labels(sample)


File /global/cscratch1/sd/ioannis/dr8/gallery/jpg/rmj150654.9+343306.1.jpg not found.
File /global/cscratch1/sd/ioannis/dr8/gallery/jpg/rmj001256.8+353148.1.jpg not found.
CPU times: user 205 ms, sys: 307 ms, total: 513 ms
Wall time: 40.3 s

To test the webpage before release, do

rsync -auvPn --delete /global/cscratch1/sd/ioannis/dr8/gallery/png /global/project/projectdirs/cosmo/www/temp/ioannis/dr8/gallery/
rsync -auvPn /global/cscratch1/sd/ioannis/dr8/gallery/*.html /global/project/projectdirs/cosmo/www/temp/ioannis/dr8/gallery/

and then navigate to

https://portal.nersc.gov/project/cosmo/temp/ioannis/dr8/gallery/


In [78]:
reject = ['rmj150654.9+343306.1', 'rmj001256.8+353148.1']
toss = np.zeros(len(sample), dtype=bool)
name = get_name(sample)
for ii, nn in enumerate(name):
    for rej in np.atleast_1d(reject):
        toss[ii] = rej in nn.lower()
        if toss[ii]:
            break
print('Rejecting {} objects.'.format(np.sum(toss)))
pngkeep = sample[~toss]
if np.sum(toss) > 0:
    pngrej = sample[toss]
else:
    pngrej = []


Rejecting 2 objects.

In [79]:
htmlfile = os.path.join(gallerydir, 'index.html')
htmlfile_reject = os.path.join(gallerydir, 'index-reject.html')
baseurl = 'http://legacysurvey.org/viewer'

In [83]:
def html_rows(pngkeep, nperrow=4):
    nrow = np.ceil(len(pngkeep) / nperrow).astype('int')
    pngsplit = list()
    for ii in range(nrow):
        i1 = nperrow*ii
        i2 = nperrow*(ii+1)
        if i2 > len(pngkeep):
            i2 = len(pngkeep)
        pngsplit.append(pngkeep[i1:i2])
    #pngsplit = np.array_split(pngkeep, nrow)
    print('Splitting the sample into {} rows with {} mosaics per row.'.format(nrow, nperrow))

    html.write('<table class="ls-gallery">\n')
    html.write('<tbody>\n')
    for pngrow in pngsplit:
        html.write('<tr>\n')
        for obj in pngrow:
            name = get_name(obj)
            nicename = get_name(obj, nice=True)
            pngfile = os.path.join('png', '{}.png'.format(name))
            thumbfile = os.path.join('png', '{}.png'.format(name))
            #thumbfile = os.path.join('png', 'thumb-{}.png'.format(name))
            img = 'src="{}" alt="{}"'.format(thumbfile, nicename)
            #img = 'class="ls-gallery" src="{}" alt="{}"'.format(thumbfile, nicename)
            html.write('<td><a href="{}"><img {} width=300px></a></td>\n'.format(pngfile, img))
        html.write('</tr>\n')
        html.write('<tr>\n')
        for obj in pngrow:
            nicename = get_name(obj, nice=True)
            href = '{}/?layer={}&ra={:.8f}&dec={:.8f}&zoom=12'.format(baseurl, dr, obj['ra'], obj['dec'])
            html.write('<td><a href="{}" target="_blank">{}</a></td>\n'.format(href, nicename))
        html.write('</tr>\n')
    html.write('</tbody>\n')            
    html.write('</table>\n')

In [84]:
with open(htmlfile, 'w') as html:
    html.write('<html><head>\n')
    html.write('<style type="text/css">\n')
    html.write('table.ls-gallery {width: 90%;}\n')
    #html.write('img.ls-gallery {display: block;}\n')
    #html.write('td.ls-gallery {width: 100%; height: auto}\n')
    #html.write('td.ls-gallery {width: 100%; word-wrap: break-word;}\n')
    html.write('p.ls-gallery {width: 80%;}\n')
    html.write('</style>\n')
    html.write('</head><body>\n')
    html.write('<h1>DR8 Image Gallery</h1>\n')
    
    html.write("""<p class="ls-gallery">This data release's gallery highlights galaxy clusters, 
    distant systems containing tens to thousands of individual galaxies and usually containing
    a single, large, dominant galaxy at the very center.  The sample presented here is taken
    from the <a href="http://risa.stanford.edu/redmapper/">redMaPPer/v6.3</a> cluster catalog.</p>\n
    
    <p class="ls-gallery">The diameter of each image cutout equals 400 kiloparsecs (more than 1.2 million light-years) at the 
    redshift of the cluster.  The clusters are sorted by increasing richness or dark matter halo mass,
    ranging between approximately 3x10<sup>15</sup> solar masses (top, left-hand side of the page)
    to 1x10<sup>14</sup> solar masses (bottom-right).</p>\n
    
    <p class="ls-gallery">The object name below each thumbnail links to the <a href="http://legacysurvey.org/viewer">Sky Viewer</a>, 
    and the horizontal white bar in the lower-right corner of each image represents 15 arcseconds.</p>\n""")
    
    #<p class="ls-gallery">For more eye candy, please visit the gallery of objects highlighted in the 
    #<a href="http://portal.nersc.gov/project/cosmo/data/legacysurvey/dr7/gallery/">DR7 Gallery.</a></p>\n""")
    
    html_rows(pngkeep)
    html.write('<br />\n')

    html.write('</body></html>\n')


Splitting the sample into 25 rows with 4 mosaics per row.

In [82]:
if len(pngrej) > 0 and False:
    with open(htmlfile_reject, 'w') as html:
        html.write('<html><head>\n')
        html.write('<style type="text/css">\n')
        html.write('img.ls-gallery {display: block;}\n')
        html.write('td.ls-gallery {width: 20%; word-wrap: break-word;}\n')
        html.write('</style>\n')
        html.write('</head><body>\n')
        html.write('<h1>DR8 Image Gallery - Rejected</h1>\n')
        html_rows(pngrej)
        html.write('</body></html>\n')

In [ ]: