Make a full overview of number counts by band by field¶
The reviewer requested that we produce number counts for each band on each field
In [1]:
%matplotlib inline
#%config InlineBackend.figure_format = 'svg'
import matplotlib as mpl
mpl.use('pdf')
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import numpy as np
#plt.rc('figure', figsize=(10, 6))
from matplotlib_venn import venn3
import herschelhelp
from herschelhelp.utils import clean_table
from astropy.table import Table, vstack
import pyvo as vo
from pymoc import MOC
import time
import yaml
import warnings
warnings.filterwarnings('ignore')
In [4]:
#Then we establish the VO connection to our database
service = vo.dal.TAPService("https://herschel-vos.phys.sussex.ac.uk/__system__/tap/run/tap"
)
In [5]:
fields = yaml.load(open('../../../dmu2/meta_main.yml', 'r'))['fields']
In [6]:
bands = [
'mmt_g',
'omegacam_g',
'suprime_g',
'megacam_g',
'wfc_g',
'gpc1_g',
'decam_g',
'90prime_g',
'sdss_g',
'isaac_k',
'moircs_k',
'ukidss_k',
'newfirm_k',
'wircs_k',
'hawki_k',
'wircam_ks',
'vista_ks',
'moircs_ks',
'omega2000_ks',
'tifkam_ks'
]
mag_tables = {}
for band in bands:
mag_tables[band] = {}
for f in fields:
mag_tables[band].update({f['name']: None})
In [80]:
start = time.time()
for band in bands:
print(band)
for f in fields:
try:
mag_tables[band].update({f['name'] : Table.read('./data/{}_{}.fits'.format(band, f['name']))})
print("loaded {} from file ({} objects)".format(band, len(mag_tables[band][f['name']])))
continue
except FileNotFoundError:
print("Querying VOX for {} mags on {}".format(band, f['name']))
query = """
SELECT
m_{}
FROM herschelhelp.main
WHERE herschelhelp.main.m_{} IS NOT NULL
AND herschelhelp.main.field='{}'""".format(
band, band,
f['name'].replace('Lockman-SWIRE','Lockman SWIRE' ).replace('HATLAS-NGP','NGP' ))
job = service.submit_job(query, maxrec=100000000)
job.run()
job_url = job.url
job_result = vo.dal.tap.AsyncTAPJob(job_url)
start_time = time.time()
wait = 10.
while job.phase == 'EXECUTING':
#print('Job still running after {} seconds.'.format(round(time.time() - start_time)))
time.sleep(wait)
#wait *=2
print(job.phase)
result = job_result.fetch_result()
mag_tables[band].update({ f['name']: result.table})
if len(mag_tables[band][f['name']]) != 0:
print("Band {} field {} done in {} seconds with {} objects".format(
band,
f['name'],
round(time.time() - start_time),
len(mag_tables[band][f['name']])
))
print("Total time: {} seconds".format(round(time.time() - start)))
In [37]:
mag_tables
Out[37]:
In [26]:
#service.maxrec=100000000
print(service.hardlimit)
Herschel Stripe 82 hits the hard limit on rows so we must get it from two queries
In [28]:
query = """
SELECT
m_{}
FROM herschelhelp.main
WHERE herschelhelp.main.m_{} IS NOT NULL
AND herschelhelp.main.field='{}'
AND herschelhelp.main.m_decam_g<25
""".format('decam_g', 'decam_g', 'Herschel-Stripe-82')
job = service.submit_job(query)
job.run()
job_url = job.url
job_result = vo.dal.tap.AsyncTAPJob(job_url)
start_time = time.time()
wait = 10.
while job.phase == 'EXECUTING':
#print('Job still running after {} seconds.'.format(round(time.time() - start_time)))
time.sleep(wait)
#wait *=2
print(job.phase)
result = job_result.fetch_result()
In [29]:
len(result.table)
Out[29]:
In [30]:
query = """
SELECT
m_{}
FROM herschelhelp.main
WHERE herschelhelp.main.m_{} IS NOT NULL
AND herschelhelp.main.field='{}'
AND herschelhelp.main.m_decam_g>25
""".format('decam_g', 'decam_g', 'Herschel-Stripe-82')
job = service.submit_job(query)
job.run()
job_url = job.url
job_result = vo.dal.tap.AsyncTAPJob(job_url)
start_time = time.time()
wait = 10.
while job.phase == 'EXECUTING':
#print('Job still running after {} seconds.'.format(round(time.time() - start_time)))
time.sleep(wait)
#wait *=2
print(job.phase)
result2 = job_result.fetch_result()
In [31]:
len(result2.table)
Out[31]:
In [34]:
hs82 = vstack([clean_table(result.table), clean_table(result2.table)])
In [35]:
len(hs82)
Out[35]:
In [36]:
mag_tables['decam_g'].update({ 'Herschel-Stripe-82': hs82})
In [7]:
write = False
read = True
for band in bands:
for f in fields:
if write and (len(mag_tables[band][f['name']]) != 0):
clean_table(mag_tables[band][f['name']]).write(
'./data/{}_{}.fits'.format(band, f['name']), overwrite = True
)
print('Table cleaned and written to ./data/{}_{}.fits'.format(band, f['name']))
elif read:
mag_tables[band].update(
{f['name'] : Table.read('./data/{}_{}.fits'.format(band, f['name']))}
)
Get depth tables to calculate area¶
In [102]:
depth_query = """
SELECT
DISTINCT
hp_idx_o_10,"""
for band in bands:
depth_query += " ferr_{}_mean,".format(band)
depth_query = depth_query.strip(',')
depth_query +=""" FROM depth.main"""
job = service.submit_job(depth_query, maxrec=100000000)
job.run()
job_url = job.url
job_result = vo.dal.tap.AsyncTAPJob(job_url)
start_time = time.time()
wait = 10.
while job.phase == 'EXECUTING':
#print('Job still running after {} seconds.'.format(round(time.time() - start_time)))
time.sleep(wait)
#wait *=2
print('Job {} after {} seconds.'.format(print(job.phase), round(time.time() - start_time)))
result = job_result.fetch_result()
depth_result = result.table
In [104]:
depth_result = clean_table(depth_result)
depth_result.write('./data/depth_result.fits', overwrite = True)
In [9]:
depth_result = Table.read('./data/depth_result.fits')
In [105]:
depth_result
Out[105]:
Plot the histograms¶
In [12]:
bands_plotting = {
'mmt_g':['MMT $g$','b'],
'omegacam_g':['Omegacam $g$','g'],
'suprime_g':['HSC $g$','r'],
'megacam_g':['Megacam $g$','c'],
'wfc_g':['WFC $g$','m'],
'gpc1_g':['GPC1 $g$','y'],
'decam_g':['DECam $g$','k'],
'90prime_g':['90Prime $g$','darkcyan'],
'sdss_g':['SDSS $g$','olive'],
'isaac_k':['ISAAC $K$','hotpink'],
'moircs_k':['MOIRCS $K$','indigo'],
'ukidss_k':['UKIDSS $K$','midnightblue'],
'newfirm_k':['Newfirm $K$','violet'],
'wircs_k':['WIRKS $K$','darkorchid'],
'hawki_k':['HAWKI $K$','royalblue'],
'wircam_ks':['WIRCam $Ks$','chartreuse'],
'vista_ks':['VISTA $Ks$','rebeccapurple'],
'moircs_ks':['MOIRCS $Ks$','gold'],
'omega2000_ks':['Omega2000 $Ks$','coral'],
'tifkam_ks':['TIFKAM $Ks$','gray']
}
In [48]:
h = np.histogram(mag_tables['decam_g']['Herschel-Stripe-82']['m_decam_g'], bins = 100)
bin_width = (np.abs(h[1][5] - h[1][4]) )
area = 100
vals = plt.fill_between( h[1][:-1], h[0]/(bin_width*area), alpha=0.5)
plt.xlim(20.,27.)
plt.yscale('log')
#plt.ylim(1.e5,1.e7)
In [50]:
[b for b in mag_tables if b.endswith('g')]
Out[50]:
In [56]:
fig, ax = plt.subplots()
f = 'COSMOS'
area = 5
for band in [b for b in mag_tables if b.endswith('g')]:
mask = np.isfinite(mag_tables[band][f]['m_'+band])
mags = mag_tables[band][f][mask]['m_'+band]
if not np.sum(mask)==0:
#vz.hist(table[name][mask], bins='scott', label=label, alpha=.5)
h = np.histogram(mags, bins = 100)
bin_width = (np.abs(h[1][5] - h[1][4]) )
#ax.fill_between( h[1][:-1], h[0]/bin_width)#, alpha=0.4)
ax.plot( h[1][:-1], h[0]/bin_width , c=bands_plotting[band])#, alpha=0.4)
ax.legend(loc=1, fontsize=8)
plt.xlim(20.,27.)
plt.xlabel("Magnitude [mag]")
plt.yscale('log')
#plt.ylim(0.,0.4)
plt.ylabel('Number [dex$^{-1}$]')
plt.rc('font', family='serif', serif='Times')
plt.rc('text') #, usetex=True)
plt.rc('xtick', labelsize=12)
plt.rc('ytick', labelsize=12)
plt.rc('axes', labelsize=12)
#plt.savefig('./figs/numbers_g_en1.pdf', bbox_inches='tight')
#plt.savefig('./figs/numbers_g_en1.png', bbox_inches='tight')
In [112]:
for n, f in enumerate(fields):
print(f['name'])
In [10]:
areas = {}
for band in bands:
areas.update({band: {}})
for n, f in enumerate(fields):
f = f['name']
f_moc = MOC(filename='../../../dmu2/dmu2_field_coverages/{}_MOC.fits'.format(f))
band_moc = MOC(10,
depth_result[~np.isnan(depth_result['ferr_{}_mean'.format(band)])]['hp_idx_o_10']
)
area = band_moc.intersection( f_moc).area_sq_deg #.flattened(order=10)
areas[band].update({f: area})
np.save('./data/areas.npy', areas)
In [8]:
areas = np.load('./data/areas.npy').item()
In [9]:
areas
Out[9]:
In [15]:
dim = [4,6]
fig, axes = plt.subplots(dim[1], dim[0], sharex=True, sharey=True)
plt.rcParams.update({'font.size': 12})
#area_per_pixel = MOC(10, (1234)).area_sq_degrees
for n, f in enumerate(fields):
f = f['name']
x, y = np.floor_divide(n, dim[0]), np.remainder(n, dim[0])
for band in [b for b in mag_tables if b.endswith('g')]:
mask = np.isfinite(mag_tables[band][f]['m_'+band])
#mask &= (bands[band][0]['field'] == f)
area = areas[band][f]
#area=f_moc.area_sq_deg
mags = mag_tables[band][f][mask]['m_'+band]
if not np.sum(mask)==0:
#vz.hist(table[name][mask], bins='scott', label=label, alpha=.5)
h = np.histogram(mags, bins = 100)
bin_width = (np.abs(h[1][5] - h[1][4]) )
#ax.fill_between( h[1][:-1], h[0]/bin_width)#, alpha=0.4)
axes[x,y].plot( h[1][:-1], h[0]/(bin_width*area), c=bands_plotting[band][1])#, alpha=0.4)
axes[x,y].get_xaxis().set_tick_params(direction='out')
axes[x,y].xaxis.set_ticks_position('bottom')
axes[x,y].tick_params(axis='x', labelsize=8)
axes[x,y].set_xlim(18, 29)
axes[x,y].set_xticks([18,20,22,24,26,28])
axes[x,y].set_ylim(1.e2, 0.5e6)
axes[x,y].set_yscale('log')
axes[x,y].scatter([-99],[-99],
label=f,
c='w', s=0.0001)
axes[x,y].legend(frameon=False, loc=(-0.2, 0.8)) #, bbox_to_anchor=(0.1, 0.7, 0.2, 0.2)
for band in [b for b in mag_tables if b.endswith('g')]:
axes[dim[1]-1,dim[0]-1].scatter([-99],[-99],
label=bands_plotting[band][0],
c=bands_plotting[band][1], s=5.)
axes[dim[1]-1,dim[0]-1].legend( prop={'size': 6},ncol=2)
axes[dim[1]-1,dim[0]-1].tick_params(axis='x', labelsize=8)
#axes[dim[1]-1,dim[0]-1].set_xlabel('band')
fig.text(0.5, 0.07, '$g$ magnitude [mag]', ha='center')
fig.text(0.04, 0.5, 'Differential $g$ number counts [deg.$^{-2}$ dex$^{-1}$]', va='center', rotation='vertical')
fig.set_size_inches(10, 12)
fig.subplots_adjust(hspace=0, wspace=0)
plt.rc('axes', labelsize=12)
plt.savefig('./figs/numbers_g_allfields.pdf', bbox_inches='tight')
plt.savefig('./figs/numbers_g_allfields.png', bbox_inches='tight')
In [16]:
dim = [4,6]
fig, axes = plt.subplots(dim[1], dim[0], sharex=True, sharey=True)
plt.rcParams.update({'font.size': 12})
#area_per_pixel = MOC(10, (1234)).area_sq_degrees
for n, f in enumerate(fields):
f = f['name']
x, y = np.floor_divide(n, dim[0]), np.remainder(n, dim[0])
f_moc = MOC(filename='../../../dmu2/dmu2_field_coverages/{}_MOC.fits'.format(f))
for band in [b for b in mag_tables if not b.endswith('g')]:
mask = np.isfinite(mag_tables[band][f]['m_'+band])
#mask &= (bands[band][0]['field'] == f)
#band_moc = MOC(10,
# depth_result[~np.isnan(depth_result['ferr_{}_mean'.format(band)])]['hp_idx_o_10']
#)
#area = band_moc.intersection( f_moc) * area_per_pixel #.flattened(order=10)
#area=f_moc.area_sq_deg
area = areas[band][f]
mags = mag_tables[band][f][mask]['m_'+band]
if not np.sum(mask)==0:
#vz.hist(table[name][mask], bins='scott', label=label, alpha=.5)
h = np.histogram(mags, bins = 100)
bin_width = (np.abs(h[1][5] - h[1][4]) )
#ax.fill_between( h[1][:-1], h[0]/bin_width)#, alpha=0.4)
axes[x,y].plot( h[1][:-1], h[0]/(bin_width*area), c=bands_plotting[band][1])#, alpha=0.4)
axes[x,y].get_xaxis().set_tick_params(direction='out')
axes[x,y].xaxis.set_ticks_position('bottom')
axes[x,y].tick_params(axis='x', labelsize=8)
axes[x,y].set_xlim(15, 28)
axes[x,y].set_xticks([16,18,20,22,24,26, 28])
axes[x,y].set_ylim(1.e1, 0.5e6)
axes[x,y].set_yscale('log')
axes[x,y].scatter([-99],[-99],
label=f,
c='w', s=0.0001)
axes[x,y].legend(frameon=False, loc=(-0.2, 0.8))
for band in [b for b in mag_tables if not b.endswith('g')]:
axes[dim[1]-1,dim[0]-1].scatter([-99],[-99],
label=bands_plotting[band][0],
c=bands_plotting[band][1], s=5.)
axes[dim[1]-1,dim[0]-1].legend( prop={'size': 6},ncol=2)
axes[dim[1]-1,dim[0]-1].tick_params(axis='x', labelsize=8)
#axes[dim[1]-1,dim[0]-1].set_xlabel('band')
fig.text(0.5, 0.07, '$K$ or $Ks$ magnitude [mag]', ha='center')
fig.text(0.04, 0.5, 'Differential $K$ or $Ks$ number counts [deg.$^{-2}$ dex$^{-1}$]', va='center', rotation='vertical')
fig.set_size_inches(10, 12)
fig.subplots_adjust(hspace=0, wspace=0)
plt.rc('axes', labelsize=12)
plt.savefig('./figs/numbers_K_allfields.pdf', bbox_inches='tight')
plt.savefig('./figs/numbers_K_allfields.png', bbox_inches='tight')
Compare COSMOS counts¶
These are taken from Laigle et al 2016 (Figure 7)
https://iopscience.iop.org/article/10.3847/0067-0049/224/2/24/pdf
In [46]:
adeep = [
[19.250836363636363, 2329.5537063254733],
[19.756290909090907, 3540.714496578173],
[20.24850909090909, 5204.3157341950655],
[20.754254545454543, 7153.942654376493],
[21.246787878787877, 9430.733111377924],
[21.752654545454543, 12432.127473600984],
[22.258448484848486, 16805.620239098735],
[22.75100606060606, 21969.39927297683],
[23.25687272727273, 28961.308633637942],
[23.74940606060606, 38178.44026370508],
[24.2424, 42926.67721223526],
[24.74909090909091, 42568.75112370236],
[25.25672727272727, 30453.440440739156],
[25.751490909090908, 18581.93382418742]
]
a,b=zip(*adeep)
adeep= np.array([a,b])
In [47]:
aud = [
[19.224024242424242, 2449.57594791443],
[19.715951515151513, 3981.071705534969],
[20.221648484848483, 5564.863876983289],
[20.71410909090909, 7522.524941724341],
[21.219951515151514, 9999.999999999989],
[21.72572121212121, 13631.556706476138],
[22.218448484848484, 16805.620239098735],
[22.72431515151515, 22154.12212581185],
[23.230278787878788, 28242.8871620805],
[23.723030303030303, 34528.826464648126],
[24.215781818181817, 42213.80945169485],
[24.722036363636363, 48671.27545242507],
[25.21541818181818, 47863.0092322638],
[25.723078787878787, 33955.4187658821],
]
a,b=zip(*aud)
aud= np.array([a,b])
In [48]:
fontana = [
[24.00249696969697, 41512.78002752295],
[24.25561212121212, 44761.95959862869],
[24.495321212121212, 49493.19092475681],
[24.748242424242424, 57064.18728389801],
[25.001454545454543, 59503.90320456443],
[25.25461818181818, 63095.7344480193],
[25.494569696969698, 64161.236839294404],
[25.747854545454544, 65244.7324492195],
[26.001066666666667, 68034.19848864731],
]
a,b=zip(*fontana)
fontana= np.array([a,b])
In [49]:
bielby = [
[19.251151515151513, 2089.2961308540366],
[19.743151515151514, 3311.311214825908],
[20.248654545454546, 4949.319092475681],
[20.7544, 6803.419848864717],
[21.24681212121212, 9352.098899885981],
[21.752581818181817, 12748.366647836887],
[22.24509090909091, 16946.925063853138],
[22.75107878787879, 21424.420855233886],
[23.2438303030303, 26192.793448156273],
[23.750278787878788, 28242.8871620805],
]
a,b=zip(*bielby)
bielby= np.array([a,b])
In [50]:
mccracken = [
[19.344339393939393, 2196.9399272976784],
[19.836339393939394, 3481.9151350201055],
[20.328557575757575, 5117.889550210738],
[20.84763636363636, 7035.139723186507],
[21.340193939393938, 9196.791985117054],
[21.832727272727272, 12123.733007482104],
[22.338593939393938, 15982.192733572023],
[22.844630303030304, 19869.265737142672],
]
a,b=zip(*mccracken)
mccracken= np.array([a,b])
In [51]:
aihara = [
[19.224024242424242, 2449.57594791443],
[19.716290909090908, 3540.714496578173],
[20.22189090909091, 5117.889550210738],
[20.714254545454544, 7153.942654376493],
[21.22002424242424, 9751.937496801149],
[21.725890909090907, 12855.557319139032],
[22.218448484848484, 16805.620239098735],
[22.724266666666665, 22528.240444714964],
[23.230157575757573, 29450.38042071756],
[23.70962424242424, 35407.14496578173],
[24.215878787878786, 40823.39234476075],
[24.709478787878787, 37231.376317700866],
[25.21798787878788, 19703.594200718988],
]
a,b=zip(*aihara)
aihara= np.array([a,b])
In [52]:
ilbert = [
[19.250545454545453, 2575.781922652323],
[19.742836363636364, 3692.0938106277963],
[20.24850909090909, 5204.3157341950655],
[20.754327272727274, 6976.480162117927],
[21.24690909090909, 9044.064195957602],
[21.752751515151516, 12022.644346174131],
[22.24530909090909, 15716.782254997685],
[22.751272727272728, 20036.33027109271],
[23.244072727272727, 24088.967285183033],
[23.750472727272726, 26413.02739580969],
]
a,b=zip(*ilbert)
ilbert= np.array([a,b])
In [89]:
fig, ax = plt.subplots()
#plt.rcParams.update({'font.size': 12})
f='COSMOS'
for band in [b for b in mag_tables if not b.endswith('g')]:
mask = np.isfinite(mag_tables[band][f]['m_'+band])
area = areas[band][f]
mags = mag_tables[band][f][mask]['m_'+band]
if not np.sum(mask)==0:
#vz.hist(table[name][mask], bins='scott', label=label, alpha=.5)
h = np.histogram(mags, bins = 100)
bin_width = (np.abs(h[1][5] - h[1][4]) )
#ax.fill_between( h[1][:-1], h[0]/bin_width)#, alpha=0.4)
ax.plot( h[1][:-1], h[0]/(bin_width*area),
c=bands_plotting[band][1],
label=bands_plotting[band][0])#, alpha=0.4)
ax.scatter(adeep[0], 2*adeep[1], label='Laigle et al. 2016 $\mathcal{A}^{Deep}$',
c='gold',s=50, alpha=0.5, edgecolors= "black")
ax.scatter(aud[0], 2*aud[1], label='Laigle et al. 2016 $\mathcal{A}^{UD}$',
c='brown',s=50, alpha=0.5, edgecolors= "black")
ax.scatter(fontana[0], 2*fontana[1], label='Fontana et al. 2014',
c='blue',s=20, alpha=0.5, edgecolors= "black", marker='D')
ax.scatter(bielby[0], 2*bielby[1], label='Bielby et al. 2012',
c='white',s=20, alpha=0.5, edgecolors= "purple", marker='s')
ax.scatter(mccracken[0], 2*mccracken[1], label='McCracken et al. 2010',
c='white',s=20, alpha=0.5, edgecolors= "dodgerblue", marker='p')
ax.scatter(aihara[0], 2*aihara[1], label='Aihara et al. 2011',
c='darkcyan',s=20, alpha=0.5, edgecolors= "chartreuse", marker='v')
ax.scatter(ilbert[0], 2*ilbert[1], label='Ilbert et al. 2013',
c='white',s=20, alpha=0.5, edgecolors= "black", marker='h')
ax.get_xaxis().set_tick_params(direction='out')
ax.xaxis.set_ticks_position('bottom')
ax.tick_params(axis='x', labelsize=8)
ax.set_xlim(19, 26)
#ax.set_xlim(15, 28)
#ax.set_xticks([16,18,20,22,24,26, 28])
#ax.set_ylim(1.e3, 1.e5)
ax.set_ylim(2.e3, 0.2e6)
ax.set_yscale('log')
#ax.legend(loc='upper center', bbox_to_anchor=(1.5, 1.05), prop={'size': 6})
ax.legend(prop={'size': 6})
ax.set_xlabel("$K$ or $Ks$ [mag]")
ax.set_ylabel("N [deg.$^{-2}$ dex$^{-1}$]")
plt.rc('font', family='serif', serif='Times')
plt.rc('text') #, usetex=True)
plt.rc('xtick', labelsize=12)
plt.rc('ytick', labelsize=12)
plt.rc('axes', labelsize=12)
column_width_cm = 8.9
width_cm = 1.0 * column_width_cm
hieght_cm = width_cm
width_inches = width_cm/2.5
hieght_inches = hieght_cm/2.5
fig.set_size_inches(width_inches, hieght_inches)
plt.savefig('./figs/numbers_K_COSMOS_comparison.pdf', bbox_inches='tight')
plt.savefig('./figs/numbers_K_COSMOS_comparison.png', bbox_inches='tight')
