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
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
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)
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]:
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)
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()
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))
In [25]:
print('Writing {}'.format(samplefile))
sample.write(samplefile, overwrite=True)
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]
Out[27]:
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()
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))
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)
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 = []
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')
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 [ ]: