The purpose of this notebook is to build the gallery for the fifth Legacy Survey data release, DR5. The theme of this gallery is "galaxy groups."
The parent sample is a diameter-limited (D25>5 arcsec) sample defined and documented as part of the Legacy Survey Large Galaxy Atlas. Here, we use the group catalog constructed in this notebook.
In [38]:
import os
import numpy as np
import fitsio
import matplotlib.pyplot as plt
from astropy.table import Table, Column, vstack
from PIL import Image, ImageDraw, ImageFont
from astrometry.util.util import Tan
from astrometry.util.fits import merge_tables
from legacypipe.survey import ccds_touching_wcs, LegacySurveyData
In [2]:
import multiprocessing
nproc = multiprocessing.cpu_count() // 2
nproc = 1 # hack because multiprocessing doesn't seem to be working in jupyter-dev
In [3]:
%matplotlib inline
In [4]:
dr = 'dr5'
if dr == 'dr5':
PIXSCALE = 0.262
else:
raise NotImplementedError
In [5]:
suffix = '0.05'
In [6]:
mindiameter, maxdiameter, d25min = 1.0, 10.0, 0.4 # [arcmin]
minnmembers = 3
In [7]:
GALAXY_DIAMFACTOR = 2.0
GROUP_DIAMFACTOR = 2.0
In [8]:
gallerydir = os.path.join( os.getenv('SCRATCH'), dr, 'gallery' )
galleryfile = os.path.join(gallerydir, 'gallery-{}.fits'.format(dr))
In [9]:
survey = LegacySurveyData()
In [10]:
bricks = survey.get_bricks()
allccds = survey.get_annotated_ccds()
cut = survey.photometric_ccds(allccds)
if cut is not None:
allccds.cut(cut)
cut = survey.ccd_cuts(allccds)
allccds.cut(cut == 0)
print('Read {} CCDs.'.format(len(allccds)))
In [11]:
def read_groupcat(gallerydir='.', mindiameter=0.0, maxdiameter=1000.0, d25min=0.5,
decmin=-90.0, decmax=+90.0, ramin=0.0, ramax=360.0,
minnmembers=1):
"""Read the LEDA group catalog. Look for individual galaxies or groups with
diameters between [mindiameter, maxdiameter] (arcminutes) but also require the
smallest member of each group to have D(25)>d25min (arcminute). These cuts
ensure that we don't build up mosaics of lots of small separated galaxies --
we require at least one "big" galaxy in each mosaic.
In addition, there are optional cuts on (ra,dec) which are useful for testing.
minnmembers is the minimum number of members in each group.
"""
groupfile = os.path.join(gallerydir, 'leda-logd25-{}-groupcat.fits'.format(suffix))
print('Reading {}'.format(groupfile))
cat = Table.read(groupfile)
these = np.where((cat['DIAMETER'] / 60.0 <= maxdiameter) *
(cat['DIAMETER'] / 60.0 >= mindiameter) *
(cat['D25MIN'] / 60.0 >= d25min) *
(cat['NMEMBERS'] >= minnmembers) *
(cat['DEC'] <= decmax) *
(cat['DEC'] >= decmin) *
(cat['RA'] <= ramax) *
(cat['RA'] >= ramin)
)[0]
outcat = cat[these]
outcat.add_column(Column(name='BRICKNAME', dtype='S17', shape=(6,), length=len(outcat)))
print('Found {}/{} groups with N>{} members, each with D(25)>{}, and in the diameter range={}-{} arcmin.'.format(
len(outcat), len(cat), minnmembers, d25min, mindiameter, maxdiameter))
return outcat
In [12]:
parent = read_groupcat(gallerydir=gallerydir, mindiameter=mindiameter,
maxdiameter=maxdiameter, d25min=d25min,
minnmembers=minnmembers)
parent
Out[12]:
In [13]:
if False:
these = (parent['RA'] > 200) * (parent['RA'] < 210) * (parent['DEC'] > 5) * (parent['DEC'] < 10.0)
parent = parent[these]
print(np.sum(these))
In [14]:
def _uniqccds(ccds):
'''Get the unique set of CCD files.'''
ccdfile = []
[ccdfile.append('{}-{}'.format(expnum, ccdname)) for expnum,
ccdname in zip(ccds.expnum, ccds.ccdname)]
_, indx = np.unique(ccdfile, return_index=True)
return ccds[indx]
In [15]:
def _galwcs(gal, factor=2.0):
'''Build a simple WCS object for a single galaxy.
'''
diam = factor*np.ceil(gal['DIAMETER']/PIXSCALE).astype('int16') # [pixels]
galwcs = Tan(gal['RA'], gal['DEC'], diam/2+0.5, diam/2+0.5,
-PIXSCALE/3600.0, 0.0, 0.0, PIXSCALE/3600.0,
float(diam), float(diam))
return galwcs
In [16]:
def _build_sample_onegalaxy(args):
"""Filler function for the multiprocessing."""
return build_sample_onegalaxy(*args)
In [17]:
def build_sample_onegalaxy(gal, verbose=False):
"""Wrapper function to find overlapping CCDs for a given galaxy.
First generously find the nearest set of CCDs that are near the galaxy and
then demand that there's 3-band coverage in a much smaller region centered
on the galaxy.
"""
#print('Working on {}...'.format(gal['GALAXY'].strip()))
if gal['NMEMBERS'] == 1:
galwcs = _galwcs(gal, factor=GALAXY_DIAMFACTOR)
else:
galwcs = _galwcs(gal, factor=GROUP_DIAMFACTOR)
these = ccds_touching_wcs(galwcs, allccds)
if len(these) > 0:
ccds1 = _uniqccds( allccds[these] )
# Is there 3-band coverage?
galwcs_small = _galwcs(gal, factor=1.0)
these_small = ccds_touching_wcs(galwcs_small, ccds1)
ccds1_small = _uniqccds( ccds1[these_small] )
if 'g' in ccds1_small.filter and 'r' in ccds1_small.filter and 'z' in ccds1_small.filter:
if verbose:
print('For {} found {} CCDs, RA = {:.5f}, Dec = {:.5f}, Diameter={:.4f} arcsec'.format(
gal['GALAXY'].strip(), len(ccds1), gal['RA'], gal['DEC'], gal['DIAMETER']))
# Also get the set of bricks touching this footprint.
rad = gal['DIAMETER'] / 3600 # [degree]
brickindx = survey.bricks_touching_radec_box(bricks,
gal['RA']-rad, gal['RA']+rad,
gal['DEC']-rad, gal['DEC']+rad)
if len(brickindx) == 0 or len(brickindx) > len(gal['BRICKNAME']):
print('This should not happen!')
import pdb ; pdb.set_trace()
gal['BRICKNAME'][:len(brickindx)] = bricks.brickname[brickindx]
return [gal, ccds1]
return None
In [18]:
def build_sample():
"""Build the sample for all galaxies with the option of multiprocessing
(although it doesn't work on jupyter-dev)."""
sampleargs = list()
for cc in parent:
sampleargs.append( (cc, False) ) # the False refers to verbose=False
if nproc > 1:
p = multiprocessing.Pool(nproc)
result = p.map(_build_sample_onegalaxy, sampleargs)
p.close()
else:
result = list()
for args in sampleargs:
result.append(_build_sample_onegalaxy(args))
# Remove non-matching galaxies and write out the sample
result = list(filter(None, result))
result = list(zip(*result))
outcat = vstack(result[0])
outccds = merge_tables(result[1])
print('Found {}/{} galaxies and groups in the DR5 footprint.'.format(len(outcat), len(parent)))
# Add a fracmasked column.
outcat.add_column(Column(name='FRACMASKED', shape=(3,), length=len(outcat), dtype='f2'))
return outcat
In [19]:
%time sample = build_sample()
sample
Out[19]:
In [20]:
fitsdir = os.path.join(gallerydir, 'fits')
if not os.path.isdir(fitsdir):
os.mkdir(fitsdir)
In [21]:
def galaxyname(gal, html=False):
"""We want to prevent the output FITS/png filenames from arbitrarily long
so we use the following convention:
1 member - use the galaxy name
2-3 members - join the individual galaxy names with a hyphen
>3 members - use the groupid with format :08d
"""
if gal['NMEMBERS'] > 500:
galaxy = 'group{:08d}'.format(gal['GROUPID'])
else:
if html:
galaxy = gal['GALAXY'].strip().lower().replace(',', ' ')
else:
galaxy = gal['GALAXY'].strip().lower().replace(',', '-')
return galaxy
In [22]:
def fits_cutouts(verbose=False):
"""Grab cutouts of every galaxy and also quickly compute the fraction of
pixels that are masked in each image (in fracmasked).
"""
for jj, gal in enumerate(sample):
galaxy = galaxyname(gal)
if gal['NMEMBERS'] == 1:
size = np.ceil(GALAXY_DIAMFACTOR*gal['DIAMETER']/PIXSCALE).astype('int16') # [pixels]
else:
size = np.ceil(GROUP_DIAMFACTOR*gal['DIAMETER']/PIXSCALE).astype('int16') # [pixels]
# Get a FITS cutout and then split the file into three individual bands.
baseurl = 'http://legacysurvey.org/viewer-dev/fits-cutout/'
fitsurl = '{}?ra={:.6f}&dec={:.6f}&pixscale={:.3f}&size={:g}&layer=decals-{}'.format(
baseurl, gal['RA'], gal['DEC'], PIXSCALE, size, dr)
fitsfile = os.path.join(fitsdir, '{}.fits'.format(galaxy))
cmd = 'wget --continue -O {:s} "{:s}"' .format(fitsfile, fitsurl)
if os.path.isfile(os.path.join(fitsdir, '{}-{}.fits'.format(galaxy, 'g'))):
if verbose:
print('Skipping galaxy {}'.format(galaxy))
else:
if verbose:
print(cmd)
os.system(cmd)
img, hdr = fitsio.read(fitsfile, ext=0, header=True)
for ii, band in enumerate(['g', 'r', 'z']):
bandfile = os.path.join(fitsdir, '{}-{}.fits'.format(galaxy, band))
#print(' Writing {}'.format(bandfile))
sample['FRACMASKED'][jj][ii] = np.sum(img[ii, :, :] == 0) / img[ii, :, :].size
fitsio.write(bandfile, img[ii, :, :], clobber=True, header=hdr)
#print(' Removing {}'.format(fitsfile))
os.remove(fitsfile)
In [23]:
%time fits_cutouts(verbose=False)
In [24]:
if os.path.isfile(galleryfile):
os.remove(galleryfile)
print('Writing {}'.format(galleryfile))
sample.write(galleryfile)
In [25]:
fraccut = 0.02 # [fraction]
In [26]:
pngdir = os.path.join(gallerydir, 'png')
if not os.path.isdir(pngdir):
os.mkdir(pngdir)
In [27]:
barlen30 = np.round(30.0 / PIXSCALE).astype('int')
fonttype = os.path.join(gallerydir, 'Georgia.ttf')
In [28]:
def make_png(pngsample, verbose=False):
"""Build a png file for each galaxy / group using Trilogy.
"""
for gal in pngsample:
galaxy = galaxyname(gal)
pngfile = os.path.join(pngdir, '{}.png'.format(galaxy))
thumbfile = os.path.join(pngdir, 'thumb-{}.png'.format(galaxy))
paramfile = os.path.join(pngdir, '{}.in'.format(galaxy))
with open(paramfile, 'w') as param:
param.write('B\n')
param.write(os.path.join(fitsdir, '{}-g.fits\n'.format(galaxy)))
param.write('\n')
param.write('G\n')
param.write(os.path.join(fitsdir, '{}-r.fits\n'.format(galaxy)))
param.write('\n')
param.write('R\n')
param.write(os.path.join(fitsdir, '{}-z.fits\n'.format(galaxy)))
param.write('\n')
param.write('indir {}\n'.format(fitsdir))
param.write('outname {}\n'.format(pngfile))
param.write('noiselum 0.3\n')
param.write('satpercent 0.001\n')
param.write('samplesize 1000\n')
param.write('stampsize 5000\n')
#param.write('scaling {}\n'.format(os.path.join(gallerydir, 'levels.txt')))
param.write('colorsatfac 0.95\n')
param.write('deletetests 1\n')
param.write('showstamps 0\n')
param.write('show 0\n')
param.write('legend 0\n')
param.write('testfirst 0\n')
param.close()
trilogy = os.path.join(os.getenv('LEGACYPIPE_DIR'), 'py', 'legacyanalysis', 'trilogy.py')
cmd = 'python {} {}'.format(trilogy, paramfile)
#print(cmd)
os.system(cmd)
for txt in ('levels.txt', 'trilogyfilterlog.txt'):
thisfile = os.path.join(os.getenv('LEGACYPIPE_DIR'), 'doc', 'nb', txt)
os.remove(thisfile)
os.remove(os.path.join(pngdir, '{}_filters.txt'.format(galaxy)))
os.remove(os.path.join(pngdir, '{}.in'.format(galaxy)))
im = Image.open(pngfile)
sz = im.size
fntsize = np.round(sz[0]/28).astype('int')
width = np.round(sz[0]/175).astype('int')
font = ImageFont.truetype(fonttype, size=fntsize)
draw = ImageDraw.Draw(im)
# Label the galaxy name--
draw.text((0+fntsize*2, 0+fntsize*2), galaxy.upper(), font=font)
#draw.text((sz[0]-fntsize*6, sz[1]-fntsize*3), galaxy.upper(), font=font)
# Add a 30 arcsec scale bar--
x0, x1, yy = sz[1]-fntsize*2-barlen30, sz[1]-fntsize*2, sz[0]-fntsize*2
draw.line((x0, yy, x1, yy), fill='white', width=width)
im.save(pngfile)
# Generate a thumbnail
cmd = '/usr/bin/convert -thumbnail 300x300 {} {}'.format(pngfile, thumbfile)
os.system(cmd)
if verbose:
print('Wrote {}'.format(pngfile))
In [29]:
indx = np.sum(sample['FRACMASKED'] < fraccut, axis=1) == 3
pngsample = sample[indx]
print('Identified {} groups with <2% pixels masked in grz.'.format(len(pngsample)))
In [30]:
pngsample
Out[30]:
In [31]:
%time make_png(pngsample, verbose=False)
To test the webpage before release, do
In [58]:
reject = ['pgc1667789-pgc1667493-pgc1666803-pgc1666735',
'pgc156324-pgc156326-pgc1151284',
'pgc200252-pgc1591421-ngc3040',
'pgc008096-ugc01617-ugc01618-pgc008116-ugc01620'
]
toss = np.zeros(len(pngsample), dtype=bool)
for ii, gg in enumerate(pngsample['GALAXY']):
for rej in np.atleast_1d(reject):
toss[ii] = rej in gg.lower().replace(',', '-')
if toss[ii]:
break
print('Rejecting {} groups.'.format(np.sum(toss)))
pngkeep = pngsample[~toss]
if np.sum(toss) > 0:
pngrej = pngsample[toss]
else:
pngrej = []
In [62]:
nperrow = 5
nrow = np.ceil(len(pngkeep) / nperrow).astype('int')
pngsplit = np.array_split(pngkeep, nrow)
print('Splitting the sample into {} rows with {} mosaics per row.'.format(nrow, nperrow))
In [63]:
htmlfile = os.path.join(gallerydir, 'index.html')
baseurl = 'http://legacysurvey.org/viewer-dev'
In [66]:
with open(htmlfile, '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>Image Gallery of Galaxy Groups</h1>\n')
html.write("""<p>This gallery highlights some of the galaxy groups in the DR5 footprint. The
groups were identified using an input catalog of galaxies with an apparent angular diameter greater than 24 arcsec
and a friends-of-friends algorithm with a linking length of 2.5 arcminute. The color mosaics
were generated using <a class="reference external" href="http://www.stsci.edu/~dcoe/trilogy">trilogy.py</a>.
Each thumbnail links to a larger image while the galaxy names below the thumbnails link to the
<a href="http://legacysurvey.org/viewer">Sky Viewer</a>. For reference, the horizontal white bar in
each image represents 30 arcsec.</p>\n""")
html.write("""<p><a href="http://portal.nersc.gov/project/cosmo/data/legacysurvey/dr4/gallery/">
Visit the DR4 Gallery.</a></p>\n""")
html.write('<table><tbody>\n')
for pngrow in pngsplit:
html.write('<tr>\n')
for gal in pngrow:
galaxy = galaxyname(gal)
pngfile = os.path.join('png', '{}.png'.format(galaxy))
thumbfile = os.path.join('png', 'thumb-{}.png'.format(galaxy))
img = 'class="ls-gallery" src="{}" alt="{}"'.format(thumbfile, galaxy.upper())
html.write('<td class="ls-gallery"><a href="{}"><img {} /></a></td>\n'.format(pngfile, img))
html.write('</tr><tr>\n')
for gal in pngrow:
galaxy = galaxyname(gal, html=True)
href = '{}/?layer=decals-{}&ra={:.8f}&dec={:.8f}&zoom=13'.format(baseurl, dr, gal['RA'], gal['DEC'])
html.write('<td class="ls-gallery"><a href="{}">{}</a></td>\n'.format(href, galaxy.upper()))
html.write('</tr>\n')
html.write('</table></tbody>')
html.write('</body></html>\n')
In [60]:
htmlfile_reject = os.path.join(gallerydir, 'index-reject.html')
In [61]:
if len(pngrej) > 0:
with open(htmlfile_reject, 'w') as html:
html.write('<html><body>\n')
for gal in pngrej:
galaxy = galaxyname(gal)
pngfile = os.path.join('png', '{}.png'.format(galaxy))
img = 'class="ls-gallery" src="{}" alt="{}"'.format(pngfile, galaxy.upper())
html.write('<a href="{}/?ra={:.8f}&dec={:.8f}&layer=decals-{}"><img {}></a>\n'.format(
baseurl, gal['RA'], gal['DEC'], dr, img))
html.write('</body></html>\n')
In [ ]: