<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Gallery for DR6
The purpose of this notebook is to build the gallery for the sixth Legacy Survey data release, DR6.  The theme of this gallery is...the NGC catalog!
For future reference
Step1:  Preliminaries
Define the data release and the various output directories.
Step2:  Read the Open NGC catalog created by Mattia Verga
Step3:  Select the desired object types.
Here we choose ...
Step4:  Require "big" objects, particularly the galaxies (to cut down the sample size).
Step5:  Convert coordinates in decimal degrees.
Step10:  Generate (find) the sample of objects in the DR6 footprint.
Step12:  Generate the color mosaics for each object.
Step14:  Add labels and a scale bar.
Step19:  Finally, assemble the webpage of good and rejected gallery images.
To test the webpage before release, do
  * rsync -auvP /global/cscratch1/sd/ioannis/dr6/gallery/png /global/project/projectdirs/cosmo/www/temp/ioannis/dr6/gallery/
  * rsync -auvP /global/cscratch1/sd/ioannis/dr6/gallery/*.html /global/project/projectdirs/cosmo/www/temp/ioannis/dr6/gallery/
<ASSISTANT_TASK:>
Python Code:
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
import astropy.units as u
from astropy.coordinates import SkyCoord
from astropy.table import Table, Column, vstack
from astropy.io import ascii
from PIL import Image, ImageDraw, ImageFont
#from astrometry.util.starutil_numpy import hmsstring2ra
from astrometry.util.util import Tan
from astrometry.util.fits import merge_tables
from legacypipe.survey import LegacySurveyData
from legacypipe.runbrick import run_brick
import multiprocessing
nproc = multiprocessing.cpu_count() // 2
%matplotlib inline
Explanation: Gallery for DR6
The purpose of this notebook is to build the gallery for the sixth Legacy Survey data release, DR6.  The theme of this gallery is...the NGC catalog!
For future reference: The notebook must be run from https://jupyter-dev.nersc.gov with the following (approximate) activation script:
```bash
!/bin/bash
version=$1                                                                                                          
connection_file=$2                                                                                                    
desiconda_version=20170818-1.1.12-img                                                                               
module use /global/common/${NERSC_HOST}/contrib/desi/desiconda/$desiconda_version/modulefiles                       
module load desiconda                                                                                                 
export LEGACY_SURVEY_DIR=/global/cscratch1/sd/dstn/dr6plus                                                          
export LEGACYPIPE_DIR=$SCRATCH/repos/legacypipe                                                                       
export PATH=$LEGACYPIPE_DIR/bin:${PATH}                                                                             
export PATH=$SCRATCH//repos/build/bin:$PATH                                                                         
export PYTHONPATH=$LEGACYPIPE_DIR/py:${PYTHONPATH}                                                                  
export PYTHONPATH=$SCRATCH/repos/build/lib/python3.5/site-packages:$PYTHONPATH                                        
module use $LEGACYPIPE_DIR/bin/modulefiles/cori                                                                     
module load dust                                                                                                      
exec python -m ipykernel -f $connection_file                                                                        
```
Imports and paths
End of explanation
dr = 'dr6'
PIXSCALE = 0.262
gallerydir = os.path.join( os.getenv('SCRATCH'), dr, 'gallery' )
galleryfile = os.path.join(gallerydir, 'gallery-{}.fits'.format(dr))
jpgdir = os.path.join(gallerydir, 'jpg')
if not os.path.isdir(jpgdir):
    os.mkdir(jpgdir)
pngdir = os.path.join(gallerydir, 'png')
if not os.path.isdir(pngdir):
    os.mkdir(pngdir)
Explanation: Preliminaries
Define the data release and the various output directories.
End of explanation
names = ('name', 'type', 'ra_hms', 'dec_dms', 'const', 'majax', 'minax', 
         'pa', 'bmag', 'vmag', 'jmag', 'hmag', 'kmag', 'sbrightn', 'hubble', 
         'cstarumag', 'cstarbmag', 'cstarvmag', 'messier', 'ngc', 'ic', 
         'cstarnames', 'identifiers', 'commonnames', 'nednotes', 'ongcnotes')
NGC = ascii.read(os.path.join(gallerydir, 'NGC.csv'), delimiter=';', names=names)
NGC.write(os.path.join(gallerydir, 'NGC.fits'), overwrite=True)
NGC
Explanation: Read the Open NGC catalog created by Mattia Verga:
https://github.com/mattiaverga/OpenNGC
bash
wget https://raw.githubusercontent.com/mattiaverga/OpenNGC/master/NGC.csv
Name: Object name composed by catalog + number
    NGC: New General Catalogue
    IC: Index Catalogue
Type: Object type
*: Star
**: Double star
*Ass: Association of stars
OCl: Open Cluster
GCl: Globular Cluster
Cl+N: Star cluster + Nebula
G: Galaxy
GPair: Galaxy Pair
GTrpl: Galaxy Triplet
GGroup: Group of galaxies
PN: Planetary Nebula
HII: HII Ionized region
DrkN: Dark Nebula
EmN: Emission Nebula
Neb: Nebula
RfN: Reflection Nebula
SNR: Supernova remnant
Nova: Nova star
NonEx: Nonexistent object
RA: Right Ascension in J2000 Epoch (HH:MM:SS.SS)
Dec: Declination in J2000 Epoch (+/-DD:MM:SS.SS)
Const: Constellation where the object is located
MajAx: Major axis, expressed in arcmin
MinAx: Minor axis, expressed in arcmin
PosAng: Major axis position angle (North Eastwards)
B-Mag: Apparent total magnitude in B filter
V-Mag: Apparent total magnitude in V filter
J-Mag: Apparent total magnitude in J filter
H-Mag: Apparent total magnitude in H filter
K-Mag: Apparent total magnitude in K filter
SurfBr (only Galaxies): Mean surface brigthness within 25 mag isophot (B-band), expressed in mag/arcsec2
Hubble (only Galaxies): Morphological type (for galaxies)
Cstar U-Mag (only Planetary Nebulae): Apparent magnitude of central star in U filter
Cstar B-Mag (only Planetary Nebulae): Apparent magnitude of central star in B filter
Cstar V-Mag (only Planetary Nebulae): Apparent magnitude of central star in V filter
M: cross reference Messier number
NGC: other NGC identification, if the object is listed twice in the catalog
IC: cross reference IC number, if the object is also listed with that identification
Cstar Names (only Planetary Nebulae): central star identifications
Identifiers: cross reference with other catalogs
Common names: Common names of the object if any
NED Notes: notes about object exported from NED
OpenNGC Notes: notes about the object data from OpenNGC catalog
End of explanation
majax = ma.getdata(NGC['majax']) # arcmin
objtype = np.char.strip(ma.getdata(NGC['type']))
keeptype = ('G', 'PN', 'OCl', 'GCl', 'Cl+N') # Cl gives us GCl, OCl, and Cl+N
#keeptype = ('G', 'GPair', 'GTrpl', 'GGroup', 'PN', 'Cl') # Cl gives us GCl, OCl, and Cl+N
keep = np.zeros(len(NGC), dtype=bool)
for otype in keeptype:
    print('Working on {}'.format(otype))
    ww = [otype == tt for tt in objtype]
    keep = np.logical_or(keep, ww)
Explanation: Select the desired object types.
Here we choose ...
End of explanation
galtoss = (objtype == 'G') * (majax < 3)
keep = np.logical_and(keep, (majax > 0.3) * (majax < 20))
keep = np.logical_and(keep, ~galtoss)
nobj = np.count_nonzero(keep)
print('Keeping {} / {} objects'.format(nobj, len(NGC)))
cat = NGC[keep]
cat
print(np.unique(ma.getdata(cat['type'])))
print(np.unique(ma.getdata(cat['hubble'])))
ww = (cat['type'] == 'G')
_ = plt.hist(cat['majax'][ww], bins=100)
Explanation: Require "big" objects, particularly the galaxies (to cut down the sample size).
End of explanation
coord = SkyCoord(ra=cat['ra_hms'], dec=cat['dec_dms'], unit=(u.hourangle, u.deg))
cat.add_column(Column(name='ra', unit='deg', length=nobj))
cat.add_column(Column(name='dec', unit='deg', length=nobj))
cat['ra'] = coord.ra.value
cat['dec'] = coord.dec.value
Explanation: Convert coordinates in decimal degrees.
End of explanation
survey = LegacySurveyData()
survey.output_dir = gallerydir
def get_name(cat, nice=False):
    name = np.atleast_1d(ma.getdata(cat['name']))
    mess = np.atleast_1d(ma.getdata(cat['messier']))
    comm = np.atleast_1d(ma.getdata(cat['commonnames']))
    outname = []
    if nice:
        hubble_type = np.empty_like(name)
        for nn, mm, cc in zip(name, mess, comm):
            oo = nn.strip().replace('NED01', '').upper()
            if mm != 0:
                oo = '{} = M{}'.format(oo, mm)
            if cc != 0:
                oo = '{} = {}'.format(oo, str(cc).replace(',', ' = '))
            outname.append(oo)
    else:
        for nn in name:
            outname.append(nn.strip().replace(' ', '_').lower())
            
    if len(outname) == 1:
        outname = outname[0]
    return outname
def simple_wcs(obj, diam):
    Build a simple WCS object for a single object.
    size = np.rint(diam / PIXSCALE).astype('int') # [pixels]
    wcs = Tan(obj['ra'], obj['dec'], size/2+0.5, size/2+0.5,
                 -PIXSCALE/3600.0, 0.0, 0.0, PIXSCALE/3600.0, 
                 float(size), float(size))
    return wcs
def _build_sample_one(args):
    Wrapper function for the multiprocessing.
    return build_sample_one(*args)
def build_sample_one(obj, factor=0.5, verbose=False):
    Wrapper function to find overlapping grz CCDs for a given object.
    
    
    name = get_name(obj)
    print('Working on {}...'.format(name))
    diam = factor * ma.getdata(obj['majax']) * 60.0 # diameter [arcsec]
    wcs = simple_wcs(obj, diam)
    try:
        ccds = survey.ccds_touching_wcs(wcs) # , ccdrad=2*diam/3600)
    except:
        return None
    
    if ccds:
        # Is there 3-band coverage?
        if 'g' in ccds.filter and 'r' in ccds.filter and 'z' in ccds.filter:
            if verbose:
                print('For {} (type={}) found {} CCDs, RA = {:.5f}, Dec = {:.5f}, Diameter={:.4f} arcmin'.format(
                        obj['name'], obj['type'], len(ccds), obj['ra'], obj['dec'], obj['majax']))
            return obj
    return None
def build_sample(cat, factor=1.0):
    Build the full sample with grz coverage in DR6.
    sampleargs = list()
    for cc in cat:
        sampleargs.append( (cc, factor, True) ) # the False refers to verbose=False
    if 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))
    # Remove non-matching objects and write out the sample
    outcat = vstack(list(filter(None, result)))
    print('Found {}/{} objects in the DR6 footprint.'.format(len(outcat), len(cat)))
    
    return outcat
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 DR6 footprint.'.format(len(sample), len(cat)))
print('Total time = {:.3f} seconds.'.format(time.time() - t0))
print('Writing {}'.format(galleryfile))
sample.write(galleryfile, overwrite=True)
sample
print(np.unique(sample['type']).data)
def qa_sample():
    fig, ax = plt.subplots()
    ax.scatter(cat['ra'], cat['dec'], alpha=0.5, s=10, label='Trimmed NGC Catalog')
    ax.scatter(sample['ra'], sample['dec'], s=20, label='Objects in DR6 Footprint')
    ax.set_xlabel('RA')
    ax.set_ylabel('Dec')
    ax.legend(loc='lower right')
qa_sample()
Explanation: Generate (find) the sample of objects in the DR6 footprint.
End of explanation
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
def get_factor(objtype):
    Scale factors for the mosaics.
    ref = dict(
        G = 2,
        GCl = 2,
        OCl = 2,
        PN = 4,
    )
    return ref[objtype]
def make_coadds_one(obj, scale=PIXSCALE, clobber=False):
    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:
        factor = get_factor(obj['type'])
        diam = factor * ma.getdata(obj['majax']) * 60.0 # diameter [arcsec]
        
        size = np.rint(diam / scale).astype('int') # [pixels]
        print('Generating mosaic for {} (type={}) with width={} pixels.'.format(name, obj['type'], size))
        
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            run_brick(None, survey, radec=(obj['ra'], obj['dec']), pixscale=scale, 
                      width=size, height=size, stages=['image_coadds'], splinesky=True,
                      early_coadds=True, pixPsf=True, hybridPsf=True, normalizePsf=True,
                      write_pickles=False, depth_cut=False, apodize=True, threads=nproc,
                      do_calibs=False, ceres=False)
            
        sys.stdout.flush()    
        brickname = custom_brickname(obj, prefix='custom-')
        _jpgfile = os.path.join(survey.output_dir, 'coadd', 'cus', brickname, 
                               'legacysurvey-{}-image.jpg'.format(brickname))
        shutil.copy(_jpgfile, jpgfile)
        shutil.rmtree(os.path.join(survey.output_dir, 'coadd'))
def make_coadds(sample, clobber=False):
    for obj in sample:
        make_coadds_one(obj, clobber=clobber)
#make_coadds_one(sample[111], clobber=True)
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):
        make_coadds(sample, clobber=False)
print('Total time = {:.3f} minutes.'.format((time.time() - t0) / 60))
Explanation: Generate the color mosaics for each object.
End of explanation
barlen = np.round(60.0 / PIXSCALE).astype('int')
fonttype = os.path.join(gallerydir, 'Georgia.ttf')
def _add_labels_one(args):
    Wrapper function for the multiprocessing.
    return add_labels_one(*args)
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))
        
    im = Image.open(jpgfile)
    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 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
    cmd = '/usr/bin/convert -thumbnail 300x300 {} {}'.format(pngfile, thumbfile)
    os.system(cmd)
def add_labels(sample):
    labelargs = list()
    for obj in sample:
        labelargs.append((obj, 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)
%time add_labels(sample)
Explanation: Add labels and a scale bar.
End of explanation
def get_type(hubble):
    Convert Hubble type to numerical type, for sorting purposes.
    numtype = {
        'E': 0,
        'E-S0': 1,
        
        'S0': 2,
        'S0-a': 3,
        
        'Sa': 4,
        'SBa': 4,
        'SABa': 4,
        'Sab': 5,
        'SBab': 5,
        
        'Sb': 6,
        'SABb': 6,
        'SBb': 6,
        'Sbc': 7,
        
        'Sc': 8,
        'SABc': 8,
        'SBc': 8,
        
        'Scd': 9,
        'SBcd': 9,
        
        'Sd': 10,
        
        'Sm': 11,
        'SBm': 11,
        
        'I': 12,
        'IAB': 12,
        'IB': 12,
        
        '0': -1
    }
    return np.array([numtype[hh] for hh in hubble])
reject = ['ngc3587', 'ngc6832', 'ngc5982', 'ngc2832', 'ngc2340', 'ngc5195',
          'ngc5308', 'ngc4346', 'ngc4036', 'ngc2681', 'ngc3718', 'ngc5377',
          'ngc2146', 'ngc3126', 'ngc2841', 'ngc2683', 'ngc4217', 'ngc4357',
          'ngc5055', 'ngc4100', 'ngc5879', 'ngc5297', 'ngc4605', 'ngc6015',
          'ngc4144', 'ngc3733', 'ngc3079', 'ngc3198', 'ngc3430', 'ngc3877',
          'ngc4062', 'ngc4631', 'ngc4656_ned01', 'ngc4395']
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 = []
htmlfile = os.path.join(gallerydir, 'index.html')
htmlfile_reject = os.path.join(gallerydir, 'index-reject.html')
baseurl = 'http://legacysurvey.org/viewer-dev'
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', '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 {}></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=decals-{}&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')
objtype = ma.getdata(pngkeep['type'])
hubbletype = get_type(ma.getdata(pngkeep['hubble']))
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>DR6 Image Gallery</h1>\n')
    html.write(<p class="ls-gallery">This gallery highlights the exquisite image quality and diversity 
    of objects observed by the Legacy Survey, including planetary nebulae, globular clusters, and 
    large, nearby galaxies.  Each thumbnail links to a larger image while the object name below each 
    thumbnail links to the 
    <a href="http://legacysurvey.org/viewer">Sky Viewer</a>.  For reference, the horizontal white bar in 
    the lower-right corner of each image represents one arcminute.</p>\n)
    html.write(<p>We gratefully acknowledge the <a href="https://github.com/mattiaverga/OpenNGC" target="_blank">
    OpenNGC</a> catalog created by Mattia Verga, which was used to generate this sample.</p>\n)
    html.write(<p>For more eye candy, please visit the gallery of galaxy groups highlighted in the 
    <a href="http://portal.nersc.gov/project/cosmo/data/legacysurvey/dr5/gallery/">DR5 Gallery.</a></p>\n)
    
    # Split by object type
    
    html.write('<h2>Planetary Nebulae, Open Clusters, and Globular Clusters</h2>\n')
    these = np.logical_or( np.logical_or(objtype == 'PN', objtype == 'OCl'), objtype == 'GCl' )
    srt = np.argsort(objtype[these])[::-1]
    html_rows(pngkeep[these][srt])
    html.write('<br />\n')
    
    html.write('<h2>Spheroidal & Elliptical Galaxies</h2>\n')
    these = (objtype == 'G') * (hubbletype <= 2)
    srt = np.argsort(hubbletype[these])
    html_rows(pngkeep[these][srt])
    html.write('<h2>Early-Type Disk Galaxies</h2>\n')
    these = (objtype == 'G') * (hubbletype >= 3) * (hubbletype <= 6)
    srt = np.argsort(hubbletype[these])
    html_rows(pngkeep[these][srt])
    html.write('<h2>Late-Type Disk Galaxies</h2>\n')
    these = (objtype == 'G') * (hubbletype >= 7) * (hubbletype <= 10)
    srt = np.argsort(hubbletype[these])
    html_rows(pngkeep[these][srt])
    
    html.write('<h2>Irregular Galaxies</h2>\n')
    these = (objtype == 'G') * (hubbletype >= 11)
    srt = np.argsort(hubbletype[these])
    html_rows(pngkeep[these][srt])
    
    html.write('</body></html>\n')
if len(pngrej) > 0:
    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>DR6 Image Gallery - Rejected</h1>\n')
        html_rows(pngrej)
        html.write('</body></html>\n')
Explanation: Finally, assemble the webpage of good and rejected gallery images.
To test the webpage before release, do
  * rsync -auvP /global/cscratch1/sd/ioannis/dr6/gallery/png /global/project/projectdirs/cosmo/www/temp/ioannis/dr6/gallery/
  * rsync -auvP /global/cscratch1/sd/ioannis/dr6/gallery/*.html /global/project/projectdirs/cosmo/www/temp/ioannis/dr6/gallery/
End of explanation
<END_TASK>