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 [ ]: