SV Fields

This example is based on a plot initially developed by Anand Raichoor.


In [ ]:
%matplotlib inline
import os
from pkg_resources import resource_filename
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
import astropy.units as u
from astropy.coordinates import SkyCoord, HeliocentricTrueEcliptic, ICRS, SkyOffsetFrame, Longitude
from astropy.io import fits
import healpy as hp
from desiutil.plots import init_sky

In [ ]:
def get_footprint(survey):
    """Define sky areas for various surveys.
    
    Parameters
    ----------
    survey : str
        The name of the survey.
    
    Returns
    -------
    dict
        A dictionary containing various sub-regions of `survey`.
    """
    mydict = {}
    if (survey=='hscdr1'):
        mydict['xmm']     = '29,40,-7,-3'
        mydict['gama09h'] = '130,138,-1,3'
        mydict['cosmos']  = '148,152,0,4'
        mydict['wide12h'] = '177,183,-2,2'
        mydict['gama15h'] = '213,221,-2,2'
        mydict['vvds']    = '331,342,-1,3'
        mydict['deep2-3'] = '350,354,-2,2'
        mydict['hectomap']= '242,249,42,45'
        mydict['aegis']   = '213,217,51,54'
        mydict['elais-n1']= '240,246,53,57'
    elif (survey=='hscdr2'):
        mydict['w05_1']   = '330,4,-2,3'
        mydict['w05_2']   = '330,345,3,7'
        mydict['w01']     = '15,23,-2,3'
        mydict['w02']     = '28,40,-8,0'
        mydict['w03']     = '126,162,-3,6'
        mydict['w04']     = '165,228,-3,6'
        mydict['w06']     = '198,252,41.5,45'
        mydict['w07']     = '212,216,51,54'
    elif (survey=='hscdr2_dud'):
        mydict['xmm']     = '33,38,-7,-3'
        mydict['cosmos']  = '148,152,0,4'
        mydict['elais']   = '239,246,53,57'
        mydict['deep2_3'] = '350,354,-2,2'
    elif (survey=='deep2'):
        mydict['deep2_1'] = '213,217,51,54'
        mydict['deep2_2'] = '251,254,34,36'
        mydict['deep2_3'] = '351,354,-1,1'
        mydict['deep2_4'] = '36,39,0,1'
    elif (survey=='kids'):
        mydict['s1']      = '0,53.5,-35.6,-25.7'
        mydict['s2']      = '329.5,360,-35.6,-25.7'
        mydict['n1']      = '156,225,-5,4'
        mydict['n2']      = '225,238,-3,4'
    elif (survey=='ebosselg'):
        mydict['eboss21'] = '315,360,-2,2'
        mydict['eboss22'] = '0,45,-5,5'
        mydict['eboss23a']= '126,143,16,29'
        mydict['eboss23b']= '136,157,13,27'
        mydict['eboss25a']= '131,166,29,32.5'
        mydict['eboss25b']= '142.5,166,27,29'
        mydict['eboss25c']= '157,166,23,27'
    elif (survey=='dr8b'):
        mydict['s82']     = '0,45,-1.25,1.25'
        mydict['hsc_sgc'] = '30,40,-6.5,-1.25'
        mydict['hsc_ngc'] = '177.5,182.5,-1,1'
        mydict['edr']     = '240,245,5,12'
        mydict['hsc_north']='240,250,42,45'
        mydict['egs']     = '213,216.5,52,54'
    elif (survey=='cosmos'):
        mydict['cosmos']  = '149,151.5,1,3.5'
    elif (survey=='vvdswide'):
        mydict['10']     = '150.4,152.0,1.0,2.8'
        mydict['14']     = '208.9,211.0,4.1,5.6'
        mydict['22']     = '333.4,335.1,-0.6,1.3'
    elif (survey=='vvdsdeep'):
        mydict['02']     = '36.2,37.1,-4.9,-4.0'
        mydict['03']     = '52.9,53.3,-28.0,-27.6'
    elif (survey=='cfhtls'):
        mydict['w1']     = '30,39,-11.5,-3.5'
        mydict['w2']     = '131,137,-7.0,-0.5'
        mydict['w3']     = '208,221,51,58'
        mydict['w4']     = '329.5,336,-1.5,5'
    elif (survey=='vipers'):
        mydict['w1']     = '30.1,38.8,-6.0,-4.15'
        mydict['w4']      = '330,335.4,0.85,2.4'
    return mydict


def get_poly_mw(string, ax):
    """Create polygon-compatible arrays from `string`.
    """
    ramin,ramax,decmin,decmax = [float(x) for x in string.split(',')] 
    raminmw,decminmw          = ax.projection_ra(np.array([ramin])), ax.projection_dec(np.array([decmin]))
    raminmw,decminmw          = raminmw[0],decminmw[0]
    ramaxmw,decmaxmw          = ax.projection_ra(np.array([ramax])), ax.projection_dec(np.array([decmax]))
    ramaxmw,decmaxmw          = ramaxmw[0],decmaxmw[0]
    tmpx    = [raminmw, ramaxmw, ramaxmw, raminmw, raminmw]
    tmpy    = [decminmw,decminmw,decmaxmw,decmaxmw,decminmw]
    tmppoly = np.zeros((len(tmpx), 2))
    tmppoly[:,0] = tmpx
    tmppoly[:,1] = tmpy
    return tmppoly

In [ ]:
#
# General survey areas.
#
surveys = ('hscdr2', 'kids', 'cfhtls', 'cosmos', 'ebosselg', 'vipers', 'vvdswide', 'vvdsdeep', 'deep2')
#
# SV fields
# https://desi.lbl.gov/trac/wiki/TargetSelectionWG/SVFields_for_DR9
#
svdict = {'01_s82': '30,40,-7,2',
          '02_egs': '210,220,50,55',
          '03_gama09': '129,141,-2,3',
          '04_gama12': '174,186,-3,2',
          '05_gama15': '212,224,-2,3',
          '06_overlap': '135,160,30,35',
          '07_refnorth': '215,230,41,46',
          '08_ages': '215,220,30,40',
          '09_sagittarius': '200,210,5,10',
          '10_highebv_n': '140,150,65,70',
          '11_highebv_s': '240,245,20,25',
          '12_highstardens_n': '273,283,40,45',
          '13_highstardens_s': '260,270,15,20',
          '14_hscw05': '330,340,-2,3'}

In [ ]:
#
# DES
#
DESfile = resource_filename('desiutil', 'data/DES_footprint.txt')
desra, desdec = np.loadtxt(DESfile, unpack=True)

In [ ]:
#
# DECaLS dr8 grz + dec>-30
#
hpdict = {}
dr8fn = os.path.join(os.environ['DESI_ROOT'], 'target', 'catalogs',
                     'dr8', '0.31.1', 'pixweight', 'pixweight-dr8-0.31.1.fits')
with fits.open(dr8fn) as hdulist:
    nside, nest = hdulist[1].header['HPXNSIDE'], hdulist[1].header['HPXNEST']
    npix = hp.nside2npix(nside)
    theta, phi = hp.pix2ang(nside, np.arange(npix, dtype=int), nest=nest)
    hpdict['ra'], hpdict['dec']  = 180./np.pi*phi, 90.-180./np.pi*theta
    for key in ('fracarea', 'stardens', 'ebv', 'psfsize_g', 'psfsize_r', 'psfsize_z',
                'galdepth_g', 'galdepth_r', 'galdepth_z', 'psfdepth_w1', 'psfdepth_w2'):
        if key == 'stardens':
            hpdict[key]       = np.log10(hdulist[1].data[key])
            hpdict[key+'lab'] = 'log10(stardens)'
        elif key.startswith('galdepth') or key.startswith('psfdepth'):
            hpdict[key]       = 22.5 - 2.5*np.log10(5./np.sqrt(hdulist[1].data[key]))
            hpdict[key+'lab'] = r'5$\sigma$ '+key
        else:
            hpdict[key]       = hdulist[1].data[key]
            hpdict[key+'lab'] = key
        if key.startswith('psfsize'):
            hpdict[key][hpdict[key]==0] = np.nan
#
# north/south 
#
c               = SkyCoord(hpdict['ra']*u.degree, hpdict['dec']*u.degree, frame='icrs')
hpdict['north'] = (hpdict['fracarea']>0) & (hpdict['dec']>32.375) & (c.galactic.b.value>0)
hpdict['south'] = ((hpdict['fracarea']>0) & (hpdict['dec']>-30)    & 
                   ((c.galactic.b.value<0) | ((c.galactic.b.value>0) & (hpdict['dec']<32.375))))
print('North: {0:.0f} deg2'.format(hpdict['north'].sum() * hp.nside2pixarea(nside, degrees=True)))
print('South: {0:.0f} deg2'.format(hpdict['south'].sum() * hp.nside2pixarea(nside, degrees=True)))
#
# SV
#
hpdict['sv']    = np.zeros(npix, dtype=bool)
tmparea         = 0
for key in sorted(svdict.keys()):
    ramin, ramax, decmin, decmax = [float(x) for x in svdict[key].split(',')]
    tmp = (hpdict['ra']>ramin) & (hpdict['ra']<ramax) & (hpdict['dec']>decmin) & (hpdict['dec']<decmax)
    hpdict['sv'][tmp] = True
    tmparea += tmp.sum() * hp.nside2pixarea(nside, degrees=True)
    print('{0}: {1:.0f} deg2'.format(key, tmp.sum() * hp.nside2pixarea(nside, degrees=True)))
print('SV: {0:.0f} deg2'.format(tmparea))

In [ ]:
ax = init_sky()
#
# dr8
#
ramw, decmw = ax.projection_ra(hpdict['ra']), ax.projection_dec(hpdict['dec'])
for s, a, l in zip([hpdict['north'], hpdict['south']],
                   [0.05, 0.6],
                   ['ls/dr8_north', 'ls/dr8_south (dec>-30)']):
    print('{0}: {1:.0f} deg2'.format(l, s.sum() * hp.nside2pixarea(nside, degrees=True)))
    foo = ax.scatter(ramw[s], decmw[s], s=1, c='y', zorder=0, alpha=a, rasterized=True, label=l)
#
# DES
#
foo = ax.plot(ax.projection_ra(desra), ax.projection_dec(desdec), c='k', lw=0.5, label='des')
#
# Other surveys
#
for survey in surveys:
    mydict = get_footprint(survey)
    count  = 0
    if survey == 'deep2':
        c, lw = 'k', 2
    else:
        c, lw = ax._get_lines.get_next_color(), 1
    for key in mydict:
        tmppoly = get_poly_mw(mydict[key], ax)
        if count == 0:
            label = survey
        else:
            label = None
        foo = ax.plot(tmppoly[:, 0], tmppoly[:, 1], c=c, lw=lw, label=label)
        count +=1
#
# SV fields
#
for key in svdict:
    tmppoly = get_poly_mw(svdict[key], ax)
    foo = ax.add_patch(Polygon(tmppoly, closed=True, fill=False, hatch='xxx', lw=2, color='k', zorder=10))
foo = ax.legend(ncol=2, loc=1)
foo.set_zorder(20)
# plt.savefig(outroot+'.skymap.png', bbox_inches='tight')

In [ ]: