<jupyter_start><jupyter_text>LegacyHalos CoaddsThe goal of this notebook is to generate LegacySurvey (*grzW1W2*) coadds of all the central galaxies (and their surrounding groups and clusters) in the legacyhalos parent sample (*legacyhalos-upenn-parent.fits*).Please note:* This notebook must be run from the jupyter-dev.nersc.gov notebook server in order to access the LegacySurvey data.* We assume that the \$LEGACYHALOS_DIR directory exists (in $CSCRATCH) with all the latest parent catalogs.* We do not yet build unWISE coadds.* We assume the DECam pixel scale (since we use just DR5 data).* The code will not handle clusters that span one or more bricks -- a ToDo would be to define custom bricks and generate custom Tractor catalogs.  Imports<jupyter_code>import os, sys, subprocess
from astropy.io import fits
from astropy.cosmology import WMAP9 as cosmo
from legacypipe.survey import LegacySurveyData
from astrometry.util.fits import fits_table
from astrometry.util.multiproc import multiproc
mp = multiproc()
source /project/projectdirs/desi/software/desi_environment.sh {}
nodes = 15
logfilename = os.path.join(reduxdir, 'run', 'logs', 'redrock.log')
print('Running redrock batch job')
print('Should take ~30 minutes; starting at {}'.format(time.asctime()))
print('Logging to {}'.format(logfilename))
with open(logfilename, 'w') as logfile:
    cmd = 'srun -A desi -N {nodes} -n {nodes} -c 64 -t 1:00:00 -C haswell --qos interactive --cpu_bind=cores'
    cmd += ' {minitestdir}/wrap-redrock --mpi --ncpu 32'
    cmd = cmd.format(nodes=nodes, minitestdir=minitestdir)
    err = subprocess.call(cmd.split(), stdout=logfile, stderr=logfile)
    if err != 0:
        print('redrock returned error {}; see {}'.format(err, logfilename))
    else:
        print('done')<jupyter_output><empty_output><jupyter_text>Paths, filenames, and other preliminaries.<jupyter_code>legacyhalos_dir = os.path.join( os.getenv('SCRATCH'), 'legacyhalos' )
dr5_dir = os.path.join(os.sep, 'project', 'projectdirs', 'cosmo', 'data', 'legacysurvey', 'dr5')
parentfile = os.path.join(legacyhalos_dir, 'legacyhalos-upenn-parent.fits')
coaddsdir = os.path.join(legacyhalos_dir, 'coadds')
if not os.path.exists(coaddsdir):
    os.makedirs(coaddsdir)
bands = ['g', 'r', 'z']
survey = LegacySurveyData(cache_dir=dr5_dir)<jupyter_output><empty_output><jupyter_text>Read the Tractor photometry from the parent catalog.<jupyter_code>bcgphot = fits_table(parentfile, ext='LSPHOT')
redcat = fits_table(parentfile, ext='REDMAPPER')
print('Read {} galaxies from {}'.format(len(bcgphot), parentfile))<jupyter_output>Converted brickname from |S8 to <U8
Converted type from |S4 to <U4
Converted wise_coadd_id from |S8 to <U8
Read 1562 galaxies from /global/cscratch1/sd/ioannis/legacyhalos/legacyhalos-upenn-parent.fits<jupyter_text>Get the desired radius of the cutout for each central.To get a sensible radius we use the quantity R_LAMBDA from the redmapper cluster catalog, which quantifies the "richness radius" in h^-1 Mpc and convert it to pixels using a standard cosmology, times a fudge factor (currently 1.5).  Note that we assume the DECam pixel scale of 0.262 arcsec/pix!Finally, we bound the radius to the interval [50, 500] pixels.<jupyter_code>def coadds_radius(redcat, pixscale=0.262, factor=1.5, rmin=50, rmax=500):
    """Get the desired radius of each cluster in pixels 
    using the R_LAMBDA, which is the richness radius in 
    h^-1 Mpc.
    
    """
    print('NB: Assuming DECam data with pixel scale = {:.3f} arcsec/pix'.format(pixscale))
    radius = redcat.r_lambda * 1e3 * cosmo.h # cluster radius in kpc
    rad_arcsec = [factor * pixscale * rad * cosmo.arcsec_per_kpc_proper(cat.z).value for 
                  rad, cat in zip(radius, redcat)]
    rad = np.array(rad_arcsec).astype('int16')
    
    rad[rad < rmin] = rmin
    rad[rad > rmax] = rmax

    return rad
radius = coadds_radius(redcat)
fig, ax = plt.subplots()
ax.scatter(redcat.z, radius, s=15, alpha=0.7)
ax.set_xlabel('Redshift')
ax.set_ylabel('Radius (pixels)')
ax.axhline(y=50, ls='--', color='k')<jupyter_output><empty_output><jupyter_text>Generate image coadds centered on each central galaxy.<jupyter_code>def coadds_stage_tims(bcg, radius=100):
    """Initialize the first step of the pipeline, returning
    a dictionary with the following keys:
    
    ['brickid', 'target_extent', 'version_header', 
     'targetrd', 'brickname', 'pixscale', 'bands', 
     'survey', 'brick', 'ps', 'H', 'ccds', 'W', 
     'targetwcs', 'tims']

    """
    from legacypipe.runbrick import stage_tims
    
    bbox = [bcg.bx-radius, bcg.bx+radius, bcg.by-radius, bcg.by+radius]
    P = stage_tims(brickname=bcg.brickname, survey=survey, target_extent=bbox,
                   pixPsf=True, hybridPsf=True, depth_cut=False, mp=mp)
    return P
def read_tractor(bcg, targetwcs, verbose=True):
    """Read the full Tractor catalog for a given brick 
    and remove the BCG.
    
    """
    H, W = targetwcs.shape

    # Read the full Tractor catalog.
    fn = survey.find_file('tractor', brick=bcg.brickname)
    cat = fits_table(fn)
    if verbose:
        print('Read {} sources from {}'.format(len(cat), fn))
    
    # Restrict to just the sources in our little box. 
    ok, xx, yy = targetwcs.radec2pixelxy(cat.ra, cat.dec)
    cat.cut(np.flatnonzero((xx > 0) * (xx < W) * (yy > 0) * (yy < H)))
    if verbose:
        print('Cut to {} sources within our box.'.format(len(cat)))
    
    # Remove the central galaxy.
    cat.cut(np.flatnonzero(cat.objid != bcg.objid))
    if verbose:
        print('Removed central galaxy with objid = {}'.format(bcg.objid))
        
    return cat
def build_model_image(cat, tims, verbose=True):
    """Generate a model image by rendering each source.
    
    """
    from legacypipe.catalog import read_fits_catalog
    from legacypipe.runbrick import _get_mod
    
    if verbose:
        print('Creating tractor sources...')
    srcs = read_fits_catalog(cat, fluxPrefix='')
    
    if False:
        print('Sources:')
        [print(' ', src) for src in srcs]

    if verbose:
        print('Rendering model images...')
    mods = [_get_mod((tim, srcs)) for tim in tims]

    return mods
def build_coadds(bcg, targetwcs, tims, mods, version_header, verbose=True):
    """Generate individual-band FITS and color coadds for each central.
    
    """
    from legacypipe.coadds import make_coadds, write_coadd_images
    from legacypipe.runbrick import rgbkwargs, rgbkwargs_resid
    from legacypipe.survey import get_rgb, imsave_jpeg
    
    if verbose:
        print('Producing coadds...')
    C = make_coadds(tims, bands, targetwcs, mods=mods, mp=mp,
                    callback=write_coadd_images,
                    callback_args=(survey, bcg.brickname, version_header, 
                                   tims, targetwcs)
                    )
    
    # Move (rename) the coadds into the desired output directory.
    for suffix in ('chi2', 'image', 'invvar', 'model'):
        for band in bands:
            oldfile = os.path.join(survey.output_dir, 'coadd', bcg.brickname[:3], 
                                   bcg.brickname, 'legacysurvey-{}-{}-{}.fits.fz'.format(
                    bcg.brickname, suffix, band))
            newfile = os.path.join(survey.output_dir, '{:05d}-{}-{}.fits.fz'.format(bcg.objid, suffix, band))
            shutil.move(oldfile, newfile)
    shutil.rmtree(os.path.join(survey.output_dir, 'coadd'))
    
    # Build png postage stamps of the coadds.
    coadd_list = [('image', C.coimgs,   rgbkwargs),
                  ('model', C.comods,   rgbkwargs),
                  ('resid', C.coresids, rgbkwargs_resid)]

    for name, ims, rgbkw in coadd_list:
        rgb = get_rgb(ims, bands, **rgbkw)
        kwa = {}
        outfn = os.path.join(survey.output_dir, '{:05d}-{}.jpg'.format(bcg.objid, name))
        if verbose:
            print('Writing {}'.format(outfn))
        imsave_jpeg(outfn, rgb, origin='lower', **kwa)
        del rgb
#for ii in range(bcgphot):
for ii in range(2):
    survey.output_dir = os.path.join(coaddsdir, '{:05d}'.format(bcgphot[ii].objid))

    # Step 1 - Set up the first stage of the pipeline.
    P = coadds_stage_tims(bcgphot[ii], radius=radius[ii])
    
    # Step 2 - Read the Tractor catalog for this brick and remove the central.
    cat = read_tractor(bcgphot[ii], P['targetwcs'])
           
    # Step 3 - Render the model images without the central.
    mods = build_model_image(cat, tims=P['tims'])

    # Step 4 - Generate and write out the coadds.
    build_coadds(bcgphot[ii], P['targetwcs'], P['tims'], mods, P['version_header'])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12267098715426455, 'fraction_numerical': 0.009952551787987502, 'mean_word_length': 3.4569365652398143, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 23, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '617167', 'license': 'MIT', 'n_tokens_mistral': 3059, 'n_tokens_neox': 2907, 'n_words': 800, 'path': 'doc/nb/legacyhalos-coadds.ipynb', 'repo_name': 'ghadeeralsheshakly/legacyhalos', 'size': '8641'}