Make depth summary figures¶
This notebook is to produce some figures for the ELAIS-N1 paper.
First we will, for every helpix 13 in HELP get the highest K or Ks depth, g depth, and irac 1 depth. We will then plot cumulative area histograms and generate the 25%, 50%, 75% areas to summarise the depths across all HELP.
In [1]:
from herschelhelp_internal import git_version
print("This notebook was run with herschelhelp_internal version: \n{}".format(git_version()))
In [2]:
import pyvo as vo
import glob
import time
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
import random
import herschelhelp as hh
from herschelhelp_internal.utils import flux_to_mag
from astropy.table import Table, Column, vstack, join, unique
from pymoc import MOC
import pandas as pd
import yaml
from pcigale.sed import SED
from pcigale.sed_modules import get_module
In [3]:
#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 [4]:
herschelhelp_python_loc = ('../../../../herschelhelp_python/')
filter_locs = glob.glob(herschelhelp_python_loc + 'database_builder/filters/*')
In [5]:
filters = [f.split('/')[-1].split('.')[0] for f in filter_locs]
filters.remove('readme')
bands = ['u', 'g', 'r', 'i', 'z', 'y',
'J', 'H', 'K', 'Ks',
'i1', 'i2', 'i3', 'i4']
In [6]:
filters
Out[6]:
In [7]:
irac_i1_query="""
SELECT DISTINCT hp_idx_O_10, ferr_ap_irac_i1_mean
FROM depth.main
WHERE ferr_ap_irac_i1_mean IS NOT NULL
"""
#Then we execute the query
#resultset = service.run_async(irac_i1_query)
job = service.submit_job(irac_i1_query)
job.run()
job_url = job.url
job_result = vo.dal.tap.AsyncTAPJob(job_url)
start_time = time.time()
while job.phase == 'EXECUTING':
print('Job still running after {} seconds.'.format(round(time.time() - start_time)))
time.sleep(10) #wait ten seconds and try again
table = job_result.fetch_result()
i1_table = table.table
In [8]:
i1_table
Out[8]:
In [9]:
plt.hist(np.log10(i1_table['ferr_ap_irac_i1_mean']) )
Out[9]:
In [10]:
irac_i1_moc = MOC(10, i1_table['hp_idx_o_10'])
irac_i1_moc.area_sq_deg
total_area = irac_i1_moc.area_sq_deg
total_area
Out[10]:
In [11]:
fig, ax = plt.subplots()
cells, fluxes = np.histogram(np.log10(i1_table['ferr_ap_irac_i1_mean']*1.e-6), bins=50)
ax.plot(fluxes[1:],
np.cumsum(cells)*total_area/cells.sum() ,
drawstyle='steps')
ax.set_xlabel('Log$_{10}$ ( Mean IRAC 1 flux error [Jy] )')
ax.set_ylabel('Cumulative area [deg.$^2$]')
ax.set_xlim([-6.5,-3.0])
#y_vals = ax.get_yticks()
#ax.set_yticklabels([n for n in y_vals])
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.4 * column_width_cm
hieght_cm = width_cm / 1.618
width_inches = width_cm/2.5
hieght_inches = hieght_cm/2.5
fig.set_size_inches(width_inches, hieght_inches)
plt.tight_layout()
plt.savefig('./figs/IRAC_i1_cumulative_area_depth.pdf')
plt.savefig('./figs/IRAC_i1_cumulative_area_depth.png')
In [12]:
i1_p25, i1_p50, i1_p75 = np.nanpercentile(i1_table['ferr_ap_irac_i1_mean'] * 3.,[25, 50, 75])
"""
The IRAC i1 band has coverage over {} square degrees with 3 sigma
depths at the 25th, 50th and 75th percentiles of {}, {}, {} respectively.""".format(total_area, i1_p25, i1_p50, i1_p75)
Out[12]:
In mags as per reviewer request
In [56]:
fig, ax = plt.subplots()
i1_pix_mag_depths, _ = flux_to_mag(i1_table['ferr_ap_irac_i1_mean']*1.e-6*5)
cells, fluxes = np.histogram(i1_pix_mag_depths, bins=1000)
#We want cumulative from faint to bright
cells =np.flip(cells)
fluxes = np.flip(fluxes)
ax.plot(fluxes[1:],
np.cumsum(cells)*total_area/cells.sum() ,
drawstyle='steps')
ax.set_xlabel('5$\sigma$ IRAC $i1$ depth [mag]')
ax.set_ylabel('Cumulative area [deg.$^2$]')
ax.set_xlim([23,16])
ax.set_ylim([0,300])
#y_vals = ax.get_yticks()
#ax.set_yticklabels([n for n in y_vals])
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/IRAC_i1_cumulative_area_depth_mag.pdf', bbox_inches='tight')
plt.savefig('./figs/IRAC_i1_cumulative_area_depth_mag.png', bbox_inches='tight')
In [14]:
mi1_p25, mi1_p50, mi1_p75 = np.nanpercentile(i1_pix_mag_depths ,[25, 50, 75])
"""
The IRAC i1 band has coverage over {} square degrees with 5 sigma
depths at the 25th, 50th and 75th percentiles of {}, {}, {} respectively.""".format(total_area, mi1_p25, mi1_p50, mi1_p75)
Out[14]:
g bands¶
There are multiple g type bands so we the the maximum depth of all of them.
In [15]:
g_filters = [f for f in filters if f.split('_')[1] == 'g']
spaced_list = ''
for g in g_filters:
spaced_list = spaced_list + ', ferr_ap_' + g + '_mean'
print(spaced_list)
In [16]:
g_bands_query="""
SELECT DISTINCT hp_idx_O_10{}
FROM depth.main
""".format(spaced_list)
g_bands_query
Out[16]:
In [17]:
#Then we execute the query
#resultset = service.run_async(irac_i1_query)
job = service.submit_job(g_bands_query)
job.run()
job_url = job.url
job_result = vo.dal.tap.AsyncTAPJob(job_url)
start_time = time.time()
while job.phase == 'EXECUTING':
print('Job still running after {} seconds.'.format(round(time.time() - start_time)))
time.sleep(10) #wait ten seconds and try again
print(job.phase)
table = job_result.fetch_result()
In [18]:
g_depth_cols = ["ferr_ap_{}_mean".format(band) for band in g_filters]
g_depth_cols
Out[18]:
In [19]:
g_table = Table(table.table)
for col in g_depth_cols:
g_table[col].fill_value = np.nan
g_table = g_table.filled()
In [20]:
np.nanmin([g_table[column] for column in g_depth_cols], axis=0)
Out[20]:
In [21]:
for col in g_depth_cols:
print(col, np.sum(np.log10(np.array(g_table[col])*1.e-6) < -9))
In [22]:
decam_low = (np.log10(np.array(g_table['ferr_ap_decam_g_mean'])*1.e-6) < -9)
prime_low = (np.log10(np.array(g_table['ferr_ap_90prime_g_mean'])*1.e-6) < -9)
In [23]:
decam_bad_moc = MOC(10, g_table['hp_idx_o_10'][decam_low])
prime_bad_moc = MOC(10, g_table['hp_idx_o_10'][decam_low])
fields = yaml.load(open('../../../dmu2/meta_main.yml', 'r'))
for field in fields['fields']:
field_moc = MOC(filename=field['region'].replace('dmu_products/', '../../../'))
decam_bad_area = decam_bad_moc.intersection(field_moc).area_sq_deg
prime_bad_area = prime_bad_moc.intersection(field_moc).area_sq_deg
print("{} has {} sq deg of bad DECam and {} sq degrees of bad 90Prime".format(field['name'],
decam_bad_area,
prime_bad_area))
In [24]:
g_table['ferr_ap_decam_g_mean'][decam_low] = g_table['ferr_ap_decam_g_mean'][decam_low] * 1.e6
g_table['ferr_ap_90prime_g_mean'][prime_low] = g_table['ferr_ap_90prime_g_mean'][prime_low] * 1.e6
In [25]:
g_table.add_column(Column(data=np.nanmin([g_table[column] for column in g_depth_cols], axis=0), name='ferr_g_min'))
In [26]:
g_table['hp_idx_o_10', 'ferr_g_min']
Out[26]:
In [27]:
g_moc = MOC(10, g_table['hp_idx_o_10'])
g_moc.area_sq_deg
total_area_g = g_moc.area_sq_deg
total_area_g
Out[27]:
In [28]:
g_table['hp_idx_o_10', 'ferr_g_min']
fig, ax = plt.subplots()
good = np.log10(np.array(g_table['ferr_g_min'])*1.e-6) > -20
good &= np.log10(np.array(g_table['ferr_g_min'])*1.e-6) <20
cells, fluxes = np.histogram(np.log10(np.array(g_table['ferr_g_min'][good ])*1.e-6), bins=100)
ax.plot(fluxes[1:],
np.cumsum(cells)*1270./cells.sum() ,
drawstyle='steps')
ax.set_xlabel('Log$_{10}$ ( Mean g flux error [Jy] )')
ax.set_ylabel('Cumulative area [deg.$^2$]')
ax.set_xlim([-8,-4.0])
#y_vals = ax.get_yticks()
#ax.set_yticklabels([n for n in y_vals])
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 / 1.618
width_inches = width_cm/2.5
hieght_inches = hieght_cm/2.5
fig.set_size_inches(width_inches, hieght_inches)
plt.tight_layout()
plt.savefig('./figs/g_cumulative_area_depth.pdf')
plt.savefig('./figs/g_cumulative_area_depth.png')
In [29]:
g_p25, g_p50, g_p75 = np.nanpercentile(g_table['ferr_g_min'][good ] * 3.,[25, 50, 75])
"""
The g band has coverage over {} square degrees with 3 sigma
depths at the 25th, 50th and 75th percentiles of {}, {}, {} respectively.""".format(1270, g_p25, g_p50, g_p75)
Out[29]:
In mags as per reviewer request
In [53]:
g_table['hp_idx_o_10', 'ferr_g_min']
fig, ax = plt.subplots()
good = np.log10(np.array(g_table['ferr_g_min'])*1.e-6) > -20
good &= np.log10(np.array(g_table['ferr_g_min'])*1.e-6) <20
g_pix_mag_depths, _ = flux_to_mag(np.array(g_table['ferr_g_min'][good ])*1.e-6*5)
cells, fluxes = np.histogram(g_pix_mag_depths, bins=1000)
cells =np.flip(cells)
fluxes = np.flip(fluxes)
ax.plot(fluxes[1:],
np.cumsum(cells)*1270./cells.sum() ,
drawstyle='steps')
ax.set_xlabel('5$\sigma$ $g$ depth [mag] ')
ax.set_ylabel('Cumulative area [deg.$^2$]')
ax.set_xlim([27,20])
ax.set_ylim([0,1300])
#y_vals = ax.get_yticks()
#ax.set_yticklabels([n for n in y_vals])
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/g_cumulative_area_depth_mag.pdf', bbox_inches='tight')
plt.savefig('./figs/g_cumulative_area_depth_mag.png', bbox_inches='tight')
In [31]:
mg_p25, mg_p50, mg_p75 = np.nanpercentile(g_pix_mag_depths,[25, 50, 75])
"""
The g band has coverage over {} square degrees with 3 sigma
depths at the 25th, 50th and 75th percentiles of {}, {}, {} respectively.""".format(1270, mg_p25, mg_p50, mg_p75)
Out[31]:
In [32]:
g_table[prime_low]
Out[32]:
K or Ks filters¶
the same for K
In [33]:
K_filters = [f for f in filters if f.split('_')[1].startswith('k')]
K_filters
Out[33]:
In [34]:
spaced_list = ''
for K in K_filters:
spaced_list = spaced_list + ', ferr_ap_' + K + '_mean'
print(spaced_list)
In [35]:
K_bands_query="""
SELECT DISTINCT hp_idx_O_10{}
FROM depth.main
""".format(spaced_list)
K_bands_query
Out[35]:
In [36]:
#Then we execute the query
#resultset = service.run_async(irac_i1_query)
job = service.submit_job(K_bands_query)
job.run()
job_url = job.url
job_result = vo.dal.tap.AsyncTAPJob(job_url)
start_time = time.time()
while job.phase == 'EXECUTING':
print('Job still running after {} seconds.'.format(round(time.time() - start_time)))
time.sleep(10) #wait ten seconds and try again
print(job.phase)
table = job_result.fetch_result()
In [37]:
depth_cols = ["ferr_ap_{}_mean".format(band) for band in K_filters]
depth_cols
Out[37]:
In [38]:
K_table = Table(table.table)
for col in depth_cols:
K_table[col].fill_value = np.nan
K_table = K_table.filled()
In [39]:
len(K_table)
Out[39]:
In [40]:
K_table.add_column(Column(data=np.nanmin([K_table[column] for column in depth_cols], axis=0), name='ferr_K_min'))
In [41]:
np.sum(np.isnan(K_table['ferr_K_min'] ))
Out[41]:
In [42]:
K_moc = MOC(10, K_table[~np.isnan(K_table['ferr_K_min'] )]['hp_idx_o_10'])
K_moc.area_sq_deg
total_area_K = K_moc.area_sq_deg
total_area_K
Out[42]:
In [43]:
K_table['hp_idx_o_10', 'ferr_K_min']
fig, ax = plt.subplots()
good = np.log10(np.array(K_table['ferr_K_min'])*1.e-6) > -20
good &= np.log10(np.array(K_table['ferr_K_min'])*1.e-6) <20
cells, fluxes = np.histogram(np.log10(np.array(K_table['ferr_K_min'][good ])*1.e-6), bins=100)
ax.plot(fluxes[1:],
np.cumsum(cells)*total_area_K/cells.sum() ,
drawstyle='steps')
ax.set_xlabel('Log$_{10}$ ( Mean K or Ks flux error [Jy] )')
ax.set_ylabel('Cumulative area [deg.$^2$]')
ax.set_xlim([-6.5,-4.5])
#y_vals = ax.get_yticks()
#ax.set_yticklabels([n for n in y_vals])
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.4 * column_width_cm
hieght_cm = width_cm / 1.618
width_inches = width_cm/2.5
hieght_inches = hieght_cm/2.5
fig.set_size_inches(width_inches, hieght_inches)
plt.tight_layout()
plt.savefig('./figs/K_cumulative_area_depth.pdf')
plt.savefig('./figs/K_cumulative_area_depth.png')
In [44]:
K_p25, K_p50, K_p75 = np.nanpercentile(K_table['ferr_K_min'][good ] * 3.,[25, 50, 75])
"""
The K band has coverage over {} square degrees with 3 sigma
depths at the 25th, 50th and 75th percentiles of {}, {}, {} respectively.""".format(total_area_K, K_p25, K_p50, K_p75)
Out[44]:
And in mags as per reviewer request
In [50]:
K_table['hp_idx_o_10', 'ferr_K_min']
fig, ax = plt.subplots()
good = np.log10(np.array(K_table['ferr_K_min'])*1.e-6) > -20
good &= np.log10(np.array(K_table['ferr_K_min'])*1.e-6) <20
#cells, fluxes = np.histogram(np.log10(np.array(K_table['ferr_K_min'][good ])*1.e-6), bins=100)
K_pix_mag_depths, _ = flux_to_mag(np.array(K_table['ferr_K_min'][good ])*1.e-6*5)
cells, fluxes = np.histogram(K_pix_mag_depths, bins=1000)
cells =np.flip(cells)
fluxes = np.flip(fluxes)
ax.plot(fluxes[1:],
np.cumsum(cells)*total_area_K/cells.sum() ,
drawstyle='steps')
ax.set_xlabel('5$\sigma$ $K$ or $Ks$ depth [mag]')
ax.set_ylabel('Cumulative area [deg.$^2$]')
ax.set_xlim([22,19])
ax.set_ylim([0,1200])
#y_vals = ax.get_yticks()
#ax.set_yticklabels([n for n in y_vals])
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/K_cumulative_area_depth_mag.pdf', bbox_inches='tight')
plt.savefig('./figs/K_cumulative_area_depth_mag.png', bbox_inches='tight')
In [49]:
K_p25, K_p50, K_p75 = np.nanpercentile(K_pix_mag_depths,[25, 50, 75])
("The K or Ks band has coverage over {} square degrees with 5 sigma" +
"depths at the 25th, 50th and 75th percentiles of {}, {}, {} respectively.").format(total_area_K, K_p25, K_p50, K_p75)
Out[49]:
Make a table of bands available on each field¶
In [50]:
bands = ['u', 'g', 'r', 'i', 'z', 'y',
'J', 'H', 'K', 'Ks',
'i1', 'i2', 'i3', 'i4']
In [51]:
spaced_list = ''
for m, band in enumerate(bands):
specific_bands = [f for f in filters if f.split('_')[1] == band.lower()]
for specific_band in specific_bands:
spaced_list = spaced_list + ', ferr_ap_' + specific_band + '_mean'
spaced_list = spaced_list.replace(' ferr_ap_lbc_u_mean,', '')
spaced_list = spaced_list.replace(' ferr_ap_lbc_y_mean,', '')
In [52]:
spaced_list
Out[52]:
In [74]:
query = """
SELECT DISTINCT hp_idx_O_10{}
FROM depth.main
""".format(spaced_list)
job = service.submit_job(query)
job.run()
job_url = job.url
job_result = vo.dal.tap.AsyncTAPJob(job_url)
start_time = time.time()
print(time.time())
while job.phase == 'EXECUTING':
print('Job still running after {} seconds.'.format(round(time.time() - start_time)))
time.sleep(50) #wait ten seconds and try again
print(job.phase)
print('Fetching result', time.time())
table = job_result.fetch_result()
In [75]:
depth_table = table.table
depth_table
Out[75]:
In [76]:
cosmos_coverage = Table.read('../../../dmu1/dmu1_ml_COSMOS/data/depths_cosmos_20190402.fits')
In [77]:
cosmos_coverage[:10].show_in_notebook()
Out[77]:
In [78]:
np.sum(~np.isnan(cosmos_coverage['ferr_irac_i1_mean']))
Out[78]:
In [79]:
plt.hist(cosmos_coverage[~np.isnan(cosmos_coverage['ferr_irac_i1_mean'])]['ferr_irac_i1_mean'], bins=100)
Out[79]:
In [80]:
#We use total here because apertur fluxes are absent
cosmos_coverage['ferr_irac_i1_mean'].name = 'ferr_ap_irac_i1_mean'
cosmos_coverage['ferr_irac_i2_mean'].name = 'ferr_ap_irac_i2_mean'
cosmos_coverage['ferr_irac_i3_mean'].name = 'ferr_ap_irac_i3_mean'
cosmos_coverage['ferr_irac_i4_mean'].name = 'ferr_ap_irac_i4_mean'
In [81]:
cosmos_coverage['hp_idx_O_10'].name='hp_idx_o_10'
In [82]:
len(cosmos_coverage)
Out[82]:
In [83]:
len(np.unique(cosmos_coverage['hp_idx_o_10']))
Out[83]:
In [84]:
cosmos_coverage.sort('hp_idx_o_10')
is_duplicate = cosmos_coverage['hp_idx_o_10'][1:] == cosmos_coverage['hp_idx_o_10'][:-1]
#np.append(is_duplicate,False)
cosmos_coverage = cosmos_coverage[~np.append(is_duplicate,False)]
In [85]:
len(cosmos_coverage)
Out[85]:
In [86]:
shared_cols = set(cosmos_coverage.colnames).intersection(set(depth_table.colnames))
In [87]:
shared_cols
Out[87]:
In [88]:
cosmos_pixels = np.in1d(depth_table['hp_idx_o_10'], np.unique(cosmos_coverage['hp_idx_o_10']))
In [89]:
np.sum(cosmos_pixels)
Out[89]:
In [90]:
len(depth_table)
Out[90]:
In [91]:
depth_table= depth_table[~cosmos_pixels]
len(depth_table)
Out[91]:
In [92]:
depth_table = vstack([depth_table, cosmos_coverage[list(shared_cols)]])
len(depth_table)
Out[92]:
In [93]:
for col in depth_table.colnames:
if depth_table[col].dtype == 'float64':
depth_table[col].fill_value = np.nan
if col.startswith('ferr'):
depth_table[col][depth_table[col] < 1.e-3] *= 1.e6
depth_table = depth_table.filled()
coverage = Table()
coverage.add_column(Column(data=depth_table['hp_idx_o_10'], name='hp_idx_o_10'))
for band in bands:
specific_bands = [f for f in filters if f.split('_')[1] == band.lower()]
print(specific_bands)
for specific_band in specific_bands:
if 'ferr_ap_{}_mean'.format(specific_band) not in depth_table.colnames:
specific_bands.remove(specific_band)
print(specific_bands)
lowest_depths = np.nanmin([depth_table[column] for column in
['ferr_ap_{}_mean'.format(s) for s in specific_bands]]
, axis=0)
coverage.add_column(Column(data=lowest_depths, name='ferr_ap_{}_mean_min'.format(band)))
In [94]:
coverage
Out[94]:
In [95]:
#for m, band in enumerate(bands):
# field_coverages.add_column(Column(data=np.full(len(field_coverages), np.nan), name='ferr_ap_{}_mean_median'.format(band)))
fields_table = Table()
for n, field in enumerate(fields['fields']):
field_moc = MOC(filename=field['region'].replace('dmu_products/', '../../../'))
field_cells= Table()
field_cells.add_column(Column(data=list(field_moc.flattened(order=10)), name='hp_idx_o_10'))
field_cells.add_column(Column(data=np.full(len(field_cells), field['name'] ), name='field'))
fields_table = vstack([fields_table, field_cells])
In [96]:
coverage = join(coverage, fields_table, join_type='left', keys='hp_idx_o_10')
In [97]:
coverage.write('./data/coverage_band_overview.fits', overwrite=True)
In [98]:
coverage = Table.read('./data/coverage_band_overview.fits')
In [99]:
coverage
Out[99]:
Make field summaries¶
In [100]:
order10_area = MOC(order=10, cells=[1]).area_sq_deg
order10_area
Out[100]:
In [101]:
#Make a list of all bands
#this was too slow so tried to run whole query at once
field_coverages = Table()
field_coverages.add_column(Column(data=[field['name'] for field in fields['fields']], name='field' ))
for m, band in enumerate(bands):
field_coverages.add_column(Column(data=np.full(len(field_coverages), np.nan),
name='{}_band_area'.format(band)))
field_coverages.add_column(Column(data=np.full(len(field_coverages), np.nan),
name='ferr_ap_{}_mean_meadian'.format(band)))
for n, field in enumerate(fields['fields']):
in_field = coverage['field'] == field['name']
for m, band in enumerate(bands):
has_band = in_field & ~np.isnan(coverage['ferr_ap_{}_mean_min'.format(band)])
area = np.sum(has_band) * order10_area
field_coverages['{}_band_area'.format(band)][field_coverages['field'] == field['name']] = area
if np.sum(has_band) != 0:
depth_median = np.nanmedian(coverage['ferr_ap_{}_mean_min'.format(band)][has_band])
depth_max = np.nanmax(coverage['ferr_ap_{}_mean_min'.format(band)][has_band])
else:
depth_median = np.nan
depth_max = np.nan
field_coverages['ferr_ap_{}_mean_meadian'.format(band)][field_coverages['field'] == field['name']] = depth_median
if band == 'i1':
print("Field: {}, {} band coverage: {}, median depth = {}, max depth = {}".format(field['name'],
band,
area,
depth_median,
depth_max))
In [102]:
field_coverages
Out[102]:
In [103]:
field_coverages['field', 'ferr_ap_i1_mean_meadian'][field_coverages['ferr_ap_i1_mean_meadian'] > 10]
Out[103]:
In [104]:
field_coverages.write('./data/field_coverages_band_overview.fits', overwrite=True)
In [105]:
field_coverages = Table.read('./data/field_coverages_band_overview.fits')
In [106]:
table = 'field '
for band in bands:
table += '& \\multicolumn{2}{c}{' + band + '}'
table += ' \\\\ \n'
table += '\\cline{1-29} \n'
for band in bands:
table += '& a & \\sigma '
table += '\\\\ \n'
In [107]:
len('AKARI-NEP&&&9.6&1.5&9.6&1.6&9.6&1.5&9.6&3.2&9.6&7.2&&&&&&&&&7.2&1.7&7.2&1.4&&&&\\'.split('&'))
Out[107]:
In [108]:
print(table)
In [109]:
for n, field in enumerate(field_coverages):
latex = ''
for element in field:
if isinstance(element , float):
if np.isnan(element):
element = ''
elif np.isclose(element, 0.0):
element = ''
else:
element = round(element,1)
latex += str(element) + '&'
print(latex[:-1] + '\\\\')
table += latex[:-1] + '\\\\ \n'
table += '\\hline'
In [110]:
print(table)
In [111]:
round(np.nan,2)
Out[111]:
In [112]:
np.sum(np.isnan(field_coverages['ferr_ap_u_mean_meadian']))
Out[112]:
Plot general overview of depths¶
In [113]:
def depths_sample(band):
mask = ~np.isnan(coverage['ferr_ap_{}_mean_min'.format(band)])
mask &= coverage['ferr_ap_{}_mean_min'.format(band)] > 0.
mask &= coverage['ferr_ap_{}_mean_min'.format(band)] <1.e3
area = (np.sum(mask)/len(coverage)) * 1270.
return np.log10(np.array(coverage['ferr_ap_{}_mean_min'.format(band)][mask]) *5.e-6 ), area
data = [ depths_sample(band)[0] for band in bands ]
areas = [ depths_sample(band)[1] for band in bands ]
In [114]:
len(data)
Out[114]:
In [115]:
# fake data
fs = 10 # fontsize
#pos = [1, 2, 4, 5, 7, 8]
#data = np.array([np.array(coverage['ferr_ap_{}_mean_min'.format(b)]) for b in bands]).T
cmap = mpl.cm.viridis
norm = mpl.colors.Normalize(vmin=np.min(areas), vmax=np.max(areas))
scalmap = mpl.cm.ScalarMappable( cmap=cmap, norm=norm)
colors = scalmap.to_rgba(areas) # The color is the angle
colors[:, 3] = 1
fig, ax = plt.subplots()
def set_axis_style(ax, labels):
ax.get_xaxis().set_tick_params(direction='out')
ax.xaxis.set_ticks_position('bottom')
ax.set_xticks(np.arange(1, len(labels) + 1))
ax.set_xticklabels(labels)
ax.set_xlim(0.25, len(labels) + 0.75)
ax.set_xlabel('band')
#ax.violinplot(np.array(coverage['ferr_ap_{}_mean_min'.format('g') ] ) )
#ax.set_title('Custom violinplot 1', fontsize=fs)
ax.set_ylabel('log10( 5$\sigma$ Depths [Jy] )')
set_axis_style(ax, ['$' + band + '$' for band in bands])
ax.set_ylim(-7, -3)
parts = ax.violinplot(data, showmeans=False, showmedians=False,
showextrema=False)
for n, part in enumerate(parts['bodies']):
part.set_facecolor(colors[n])
part.set_alpha(1)
cax, _ = mpl.colorbar.make_axes(ax)
n_ticks = 7
values = np.linspace(0,1200, n_ticks)
ticks = values/np.max(areas)
cbar = mpl.colorbar.ColorbarBase(cax, cmap=cmap, ticks = ticks)
cax.set_yticklabels([int(d) for d in values])
cbar.set_label('area [square degrees]')
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)
#fig.suptitle("Violin Plotting Examples")
#fig.subplots_adjust(hspace=0.4)
#plt.ylim(-10,10)
column_width_cm = 8.9
width_cm = 3.0 * column_width_cm
hieght_cm = width_cm / 1.9
width_inches = width_cm/2.5
hieght_inches = hieght_cm/2.5
fig.set_size_inches(width_inches, hieght_inches)
plt.savefig('./figs/band_depths_overviews.pdf', bbox_inches='tight')
plt.savefig('./figs/band_depths_overviews.png', bbox_inches='tight')
In [116]:
widths = [np.max(np.histogram(data[n], bins = 100, density = True)[0])*areas[n] for n in np.arange(len(data)) ]
widths /= np.max(widths)
print(widths)
In [117]:
# fake data
fs = 10 # fontsize
#pos = [1, 2, 4, 5, 7, 8]
#data = np.array([np.array(coverage['ferr_ap_{}_mean_min'.format(b)]) for b in bands]).T
cmap = mpl.cm.viridis
norm = mpl.colors.Normalize(vmin=np.min(areas), vmax=np.max(areas))
scalmap = mpl.cm.ScalarMappable( cmap=cmap, norm=norm)
colors = scalmap.to_rgba(areas) # The color is the angle
colors[:, 3] = 1
fig, ax = plt.subplots()
def set_axis_style(ax, labels):
ax.get_xaxis().set_tick_params(direction='out')
ax.xaxis.set_ticks_position('bottom')
ax.set_xticks(np.arange(1, len(labels) + 1))
ax.set_xticklabels(labels)
ax.set_xlim(0.25, len(labels) + 0.75)
ax.set_xlabel('band')
#ax.violinplot(np.array(coverage['ferr_ap_{}_mean_min'.format('g') ] ) )
#ax.set_title('Custom violinplot 1', fontsize=fs)
ax.set_ylabel('log10( 5$\sigma$ Depths [Jy] )')
set_axis_style(ax, ['$' + band + '$' for band in bands])
ax.set_ylim(-7, -3)
#areas/np.max(areas)
parts = ax.violinplot(data, widths=widths, showmeans=False, showmedians=False,
showextrema=False)
for n, part in enumerate(parts['bodies']):
part.set_facecolor(colors[n])
part.set_alpha(1)
cax, _ = mpl.colorbar.make_axes(ax)
n_ticks = 7
values = np.linspace(0,1200, n_ticks)
ticks = values/np.max(areas)
cbar = mpl.colorbar.ColorbarBase(cax, cmap=cmap, ticks = ticks)
cax.set_yticklabels([int(d) for d in values])
cbar.set_label('Area [square degrees]')
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)
#fig.suptitle("Violin Plotting Examples")
#fig.subplots_adjust(hspace=0.4)
#plt.ylim(-10,10)
column_width_cm = 8.9
width_cm = 3.0 * column_width_cm
hieght_cm = width_cm / 1.9
width_inches = width_cm/2.5
hieght_inches = hieght_cm/2.5
fig.set_size_inches(width_inches, hieght_inches)
plt.savefig('./figs/band_depths_overviews_areaweighted.pdf', bbox_inches='tight')
plt.savefig('./figs/band_depths_overviews_areaweighted.png', bbox_inches='tight')
In [118]:
# fake data
fs = 10 # fontsize
#pos = [1, 2, 4, 5, 7, 8]
#data = np.array([np.array(coverage['ferr_ap_{}_mean_min'.format(b)]) for b in bands]).T
cmap = mpl.cm.viridis
norm = mpl.colors.Normalize(vmin=np.min(areas), vmax=np.max(areas))
scalmap = mpl.cm.ScalarMappable( cmap=cmap, norm=norm)
colors = scalmap.to_rgba(areas) # The color is the angle
colors[:, 3] = 1
fig, ax = plt.subplots()
def set_axis_style(ax, labels):
ax.get_xaxis().set_tick_params(direction='out')
ax.xaxis.set_ticks_position('bottom')
ax.set_xticks(np.arange(1, len(labels) + 1))
ax.set_xticklabels(labels)
ax.set_xlim(0.25, len(labels) + 0.75)
ax.set_xlabel('band')
#ax.violinplot(np.array(coverage['ferr_ap_{}_mean_min'.format('g') ] ) )
#ax.set_title('Custom violinplot 1', fontsize=fs)
ax.set_ylabel('log10( 5$\sigma$ Depths [Jy] )')
set_axis_style(ax, ['$' + band + '$' for band in bands])
ax.set_ylim(-7, -3)
parts = ax.violinplot(data, widths=np.log10(areas)/np.log10(np.max(areas)), showmeans=False, showmedians=False,
showextrema=False)
for n, part in enumerate(parts['bodies']):
part.set_facecolor(colors[n])
part.set_alpha(1)
cax, _ = mpl.colorbar.make_axes(ax)
n_ticks = 7
values = np.linspace(0,1200, n_ticks)
ticks = values/np.max(areas)
cbar = mpl.colorbar.ColorbarBase(cax, cmap=cmap, ticks = ticks)
cax.set_yticklabels([int(d) for d in values])
cbar.set_label('area [square degrees]')
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)
#fig.suptitle("Violin Plotting Examples")
#fig.subplots_adjust(hspace=0.4)
#plt.ylim(-10,10)
column_width_cm = 8.9
width_cm = 3.0 * column_width_cm
hieght_cm = width_cm / 1.9
width_inches = width_cm/2.5
hieght_inches = hieght_cm/2.5
fig.set_size_inches(width_inches, hieght_inches)
plt.savefig('./figs/band_depths_overviews_logareaweighted.pdf', bbox_inches='tight')
plt.savefig('./figs/band_depths_overviews_logareaweighted.png', bbox_inches='tight')
Add SED to plot¶
In [119]:
#s.add_contribution("HELP_SED", orig_spec['WAVE'], orig_spec['LUMIN'])
In [120]:
cigale_filternames = {
'u': 'u_prime',
'g': 'MCam_g',
'r': 'MCam_r',
'i': 'MCam_i',
'z': 'MCam_z',
'y': 'WFCAM_Y',
'J': 'WFCAM_J',
'H': 'WFI_H',
'K': 'WFI_K',
'Ks': 'Ks_2mass',
'i1': 'IRAC1',
'i2': 'IRAC2',
'i3': 'IRAC3',
'i4': 'IRAC4'
}
In [121]:
gal1 = './data/HELP_J095946.083p021914.438_best_model.fits'
gal2 = './data/HELP_J100130.443p020929.494_best_model.fits'
gal3 = './data/HELP_J095809.302p013203.775_best_model.fits'
gal4 = './data/HELP_J095822.986p013145.336_best_model.figs'
gal5 = './data/HELP_J003412.527-441056.846_best_model.fits' #z=4
gal6 = './data/HELP_J095738.934+021508.530_best_model.fits' #From hedam - high fluxes
gal7 = './data/HELP_J095818.598+013057.910_best_model.fits'
gal8 = './data/HELP_J095809.488+013225.513_best_model.fits'
gal9 = './data/HELP_J100108.952+022730.528_best_model.fits'
orig_spec = Table.read(gal9)
s = SED()
# This is wrong because the best SED we get from CIGALE is redshifted (written by Yannick)
s.add_contribution("HELP_SED", orig_spec['wavelength'], orig_spec['L_lambda_total'])
#z=1
zs = np.arange(1, 5, 1)
gal_fluxes = np.full([len(zs), len(bands)], np.nan)
for n, z in enumerate(zs):
print('z = {}:'.format(z))
sed = s.copy()
mod = get_module("redshifting", redshift=z)
mod.process(sed)
for m, band in enumerate(cigale_filternames):
flux = sed.compute_fnu(cigale_filternames[band])
print("{} band: {} mJy".format(band,flux))
gal_fluxes[n, m] = flux
In [122]:
fig, ax = plt.subplots()
for n, z in enumerate(zs):
print('z = {}:'.format(z))
sed = s.copy()
mod = get_module("redshifting", redshift=z)
mod.process(sed)
ax.plot(np.log10(sed.wavelength_grid), np.log10(sed.fnu * 1.e-3), label=str(z))
#ax.plot(np.log10(orig_spec['wavelength']), np.log10(orig_spec['Fnu']))
plt.legend()
ax.set_ylim(-7, 0)
#ax.set_xlim(0, 14.75)
Out[122]:
In [123]:
gal_fluxes
Out[123]:
In [124]:
pos = [3561., #u (SDSS)
4866., #g (GPC1)
6215., #r (GPC1)
7545., #i (GPC1)
8680., #z (GPC1)
9633., #y (GPC1)
12510., #J (UKIRT)
16377., #H (UKIRT)
22081., #K (UKIRT)
21496., #Ks (WIRCam)
36000., #i1
45000., #i2
56000., #i3
80000. #i4
]
fwhms = [
[3048., 4028.], #u (SDSS)
[3943., 5593.], #g (GPC1)
[5386., 7036.], #r (GPC1)
[6778., 8304.], #i (GPC1)
[8028., 9346.], #z (GPC1)
[9100., 10838.], #y (GPC1)
[11690., 13280], #J (UKIRT)
[14920., 17840.], #H (UKIRT)
[20290., 23800.], #K (UKIRT)
[19578., 23431.], #Ks (WIRCam)
[31296, 39614 ], #i1
[39173, 50561], #i2
[48983, 65089], #i3
[62994, 95876] #i4
]
In [125]:
fig, ax = plt.subplots()
line_styles = [':', '-.', '--', '-']
colours = ['b', 'g', 'r', 'k']
for n, z in enumerate(zs):
sed = s.copy()
mod = get_module("redshifting", redshift=z)
mod.process(sed)
ax.plot(sed.wavelength_grid*10,
np.log10(sed.fnu * 1.e-3),
#c='k',
c= colours[n],
#linestyle = line_styles[n],
)
for m, band in enumerate(cigale_filternames):
if m == 0:
lab = 'z = {}'.format(z)
else:
lab=None
ax.plot([fwhms[m][0], fwhms[m][1]], [np.log10(gal_fluxes[n, m] )-3,
np.log10(gal_fluxes[n, m] )-3],
# c='k',
c= colours[n],
#linestyle = line_styles[n],
label=lab
)
plt.legend()
plt.xscale('log')
ax.set_ylim(-7, -3)
ax.set_xlim(3000., 100000)
Out[125]:
In [126]:
# fake data
fs = 10 # fontsize
# Angstroms
#data = np.array([np.array(coverage['ferr_ap_{}_mean_min'.format(b)]) for b in bands]).T
cmap = mpl.cm.viridis
norm = mpl.colors.Normalize(vmin=np.min(areas), vmax=np.max(areas))
scalmap = mpl.cm.ScalarMappable( cmap=cmap, norm=norm)
colors = scalmap.to_rgba(areas) # The color is the angle
colors[:, 3] = 1
fig, ax = plt.subplots()
def set_axis_style(ax, labels, positions=False):
ax.get_xaxis().set_tick_params(direction='out')
ax.xaxis.set_ticks_position('bottom')
if not positions:
ax.set_xticks(np.arange(1, len(labels) + 1))
ax.set_xticklabels(labels)
ax.set_xlim(0.25, len(labels) + 0.75)
else:
ax.set_xticks(positions)
ax.set_xticklabels(labels)
ax.set_xlim(3000., 100000)
plt.xscale('log')
ax.set_xlabel('band')
#ax.violinplot(np.array(coverage['ferr_ap_{}_mean_min'.format('g') ] ) )
#ax.set_title('Custom violinplot 1', fontsize=fs)
ax.set_ylabel('log10( 5$\sigma$ Depths [Jy] )')
set_axis_style(ax, ['$' + band + '$' for band in bands])
ax.set_ylim(-7, -3)
#widths
log_widths = np.ones(len(pos)) * (pos) * .2 * widths
#areas/np.max(areas)
parts = ax.violinplot(data,
#positions=pos,
widths=widths,
showmeans=False,
showmedians=False, #widths=widths,
showextrema=False)
for n, part in enumerate(parts['bodies']):
part.set_facecolor(colors[n])
part.set_alpha(1.)
cax, _ = mpl.colorbar.make_axes(ax)
n_ticks = 7
values = np.linspace(0,1200, n_ticks)
ticks = values/np.max(areas)
cbar = mpl.colorbar.ColorbarBase(cax, cmap=cmap, ticks = ticks)
cax.set_yticklabels([int(d) for d in values])
cbar.set_label('Area [square degrees]')
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)
#fig.suptitle("Violin Plotting Examples")
#fig.subplots_adjust(hspace=0.4)
#plt.ylim(-10,10)
column_width_cm = 8.9
width_cm = 3.0 * column_width_cm
hieght_cm = width_cm / 1.9
width_inches = width_cm/2.5
hieght_inches = hieght_cm/2.5
fig.set_size_inches(width_inches, hieght_inches)
line_styles = [':', '-.', '--', '-']
for n, z in enumerate(zs):
for m, band in enumerate(cigale_filternames):
if m == 0:
lab = 'z = {}'.format(z)
else:
lab=None
ax.plot([m + 0.5, m + 1.5], [np.log10(gal_fluxes[n, m] *1.e-3 ),
np.log10(gal_fluxes[n, m] *1.e-3 )],
c='k',
linestyle = line_styles[n],
label=lab
)
ax.legend(loc=2)
plt.savefig('./figs/band_depths_overviews_areaweighted_with_seds.pdf', bbox_inches='tight')
plt.savefig('./figs/band_depths_overviews_areaweighted_with_seds.png', bbox_inches='tight')
In [127]:
# fake data
fs = 10 # fontsize
# Angstroms
#data = np.array([np.array(coverage['ferr_ap_{}_mean_min'.format(b)]) for b in bands]).T
cmap = mpl.cm.viridis
norm = mpl.colors.Normalize(vmin=np.min(areas), vmax=np.max(areas))
scalmap = mpl.cm.ScalarMappable( cmap=cmap, norm=norm)
colors = scalmap.to_rgba(areas) # The color is the angle
colors[:, 3] = 1
fig, ax = plt.subplots()
ax2 = ax.twiny()
#ax2.set_xlim(3000., 100000)
#ax2.set_xticks(pos)
#ax2.set_xticklabels(['$' + band + '$' for band in bands])
#ax2.set_xlabel('band')
#ax2.set_xscale('log')
ax2.set_xlim(np.log10(3000.), np.log10(100000))
ax2.set_xticks(np.log10(pos))
ax2.set_xticklabels(['${}$'.format(band).replace('$Ks$', '').replace('K', 'K/Ks') for band in bands])
ax2.set_xlabel('band')
#ax2.set_xscale('log')
line_styles = [':', '-.', '--', '-']
colours = ['y', 'b', 'g', 'r', 'k']
for n, z in enumerate([0.25, 1, 2, 3, 4]):
sed = s.copy()
mod = get_module("redshifting", redshift=z)
mod.process(sed)
ax.plot(sed.wavelength_grid*10,
np.log10(sed.fnu * 1.e-3),
#c='k',
c= colours[n],
#linestyle = line_styles[n],
label= 'z = {}'.format(z),
alpha=0.5
)
for m, band in enumerate(cigale_filternames):
continue
if m == 0:
lab = 'z = {}'.format(z)
else:
lab=None
ax.plot([fwhms[m][0], fwhms[m][1]], [np.log10(gal_fluxes[n, m] )-3,
np.log10(gal_fluxes[n, m] )-3],
# c='k',
c= colours[n],
#linestyle = line_styles[n],
label=lab,
alpha=0.7
)
ax.set_ylabel('log10( 5$\sigma$ Depths [Jy] )')
#ax.set_xticks(pos)
#ax.set_xticklabels(['${}$'.format(band.replace('Ks', '').replace('K', 'K/Ks')) for band in bands])
ax.set_xlim(3000., 100000)
ax.set_xscale('log')
ax.set_xlabel('Wavelength [nm]')
ax.set_ylim(-7, -3)
ax.legend(loc=2)
#widths
log_widths = np.ones(len(pos)) * (pos) * .2 * widths
#areas/np.max(areas)
parts = ax.violinplot(data,
positions=pos,
widths=log_widths,
showmeans=False,
showmedians=False, #widths=widths,
showextrema=False)
for n, part in enumerate(parts['bodies']):
part.set_facecolor(colors[n])
part.set_alpha(.9)
cax, _ = mpl.colorbar.make_axes([ax, ax2])
n_ticks = 7
values = np.linspace(0,1200, n_ticks)
ticks = values/np.max(areas)
cbar = mpl.colorbar.ColorbarBase(cax, cmap=cmap, ticks = ticks)
cax.set_yticklabels([int(d) for d in values])
cbar.set_label('Area [square degrees]')
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)
#fig.suptitle("Violin Plotting Examples")
#fig.subplots_adjust(hspace=0.4)
#plt.ylim(-10,10)
column_width_cm = 8.9
width_cm = 3.0 * column_width_cm
hieght_cm = width_cm / 1.9
width_inches = width_cm/2.5
hieght_inches = hieght_cm/2.5
fig.set_size_inches(width_inches, hieght_inches)
plt.savefig('./figs/band_depths_overviews_areaweighted_with_seds_wave.pdf', bbox_inches='tight')
plt.savefig('./figs/band_depths_overviews_areaweighted_with_seds_wave.png', bbox_inches='tight')
In [128]:
# fake data
fs = 10 # fontsize
# Angstroms
#data = np.array([np.array(coverage['ferr_ap_{}_mean_min'.format(b)]) for b in bands]).T
cmap = mpl.cm.viridis
norm = mpl.colors.Normalize(vmin=np.min(areas), vmax=np.max(areas))
scalmap = mpl.cm.ScalarMappable( cmap=cmap, norm=norm)
colors = scalmap.to_rgba(areas) # The color is the angle
colors[:, 3] = 1
fig, ax = plt.subplots()
ax2 = ax.twiny()
#ax2.set_xlim(3000., 100000)
#ax2.set_xticks(pos)
#ax2.set_xticklabels(['$' + band + '$' for band in bands])
#ax2.set_xlabel('band')
#ax2.set_xscale('log')
ax2.set_xlim(np.log10(3000.), np.log10(100000))
ax2.set_xticks(np.log10(pos))
ax2.set_xticklabels(['${}$'.format(band).replace('$Ks$', '').replace('K', 'K/Ks') for band in bands])
ax2.set_xlabel('band')
#ax2.set_xscale('log')
line_styles = [':', '-.', '--', '-']
colours = ['y', 'b', 'g', 'r', 'k']
for n, z in enumerate([0.25, 1, 2, 3, 4]):
sed = s.copy()
mod = get_module("redshifting", redshift=z)
mod.process(sed)
ax.plot(sed.wavelength_grid*10,
np.log10(sed.fnu * 1.e-3),
c='k',
#c= colours[n],
#linestyle = line_styles[n],
label= 'z = {}'.format(z),
linewidth=1.0,
alpha=1.
)
for m, band in enumerate(cigale_filternames):
continue
if m == 0:
lab = 'z = {}'.format(z)
else:
lab=None
ax.plot([fwhms[m][0], fwhms[m][1]], [np.log10(gal_fluxes[n, m] )-3,
np.log10(gal_fluxes[n, m] )-3],
# c='k',
c= colours[n],
#linestyle = line_styles[n],
label=lab,
alpha=1.0
)
ax.set_ylabel('log10( 5$\sigma$ Depths [Jy] )')
#ax.set_xticks(pos)
#ax.set_xticklabels(['${}$'.format(band.replace('Ks', '').replace('K', 'K/Ks')) for band in bands])
ax.set_xlim(3000., 100000)
ax.set_xscale('log')
ax.set_xlabel('Wavelength [nm]')
ax.set_ylim(-7, -3)
#ax.legend(loc=2)
#widths
log_widths = np.ones(len(pos)) * (pos) * .2 * widths
#areas/np.max(areas)
parts = ax.violinplot(data,
positions=pos,
widths=log_widths,
showmeans=False,
showmedians=False, #widths=widths,
showextrema=False)
for n, part in enumerate(parts['bodies']):
part.set_facecolor(colors[n])
part.set_alpha(1.0)
cax, _ = mpl.colorbar.make_axes([ax, ax2])
n_ticks = 7
values = np.linspace(0,1200, n_ticks)
ticks = values/np.max(areas)
cbar = mpl.colorbar.ColorbarBase(cax, cmap=cmap, ticks = ticks)
cax.set_yticklabels([int(d) for d in values])
cbar.set_label('Area [square degrees]')
plt.figtext(0.49, 0.83, 'z = 0.25')
plt.figtext(0.71, 0.645, 'z = 1')
plt.figtext(0.71, 0.51, 'z = 2')
plt.figtext(0.71, 0.425, 'z = 3')
plt.figtext(0.71, 0.365, 'z = 4')
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)
#fig.suptitle("Violin Plotting Examples")
#fig.subplots_adjust(hspace=0.4)
#plt.ylim(-10,10)
column_width_cm = 8.9
width_cm = 3.0 * column_width_cm
hieght_cm = width_cm / 1.9
width_inches = width_cm/2.5
hieght_inches = hieght_cm/2.5
fig.set_size_inches(width_inches, hieght_inches)
plt.savefig('./figs/band_depths_overviews_areaweighted_with_black_seds_wave.pdf', bbox_inches='tight')
plt.savefig('./figs/band_depths_overviews_areaweighted_with_black_seds_wave.png', bbox_inches='tight')
In [129]:
# fake data
fs = 10 # fontsize
# Angstroms
#data = np.array([np.array(coverage['ferr_ap_{}_mean_min'.format(b)]) for b in bands]).T
cmap = mpl.cm.viridis
norm = mpl.colors.Normalize(vmin=np.min(areas), vmax=np.max(areas))
scalmap = mpl.cm.ScalarMappable( cmap=cmap, norm=norm)
colors = scalmap.to_rgba(areas) # The color is the angle
colors[:, 3] = 1
fig, ax = plt.subplots()
ax2 = ax.twiny()
#ax2.set_xlim(3000., 100000)
#ax2.set_xticks(pos)
#ax2.set_xticklabels(['$' + band + '$' for band in bands])
#ax2.set_xlabel('band')
#ax2.set_xscale('log')
ax2.set_xlim(np.log10(300.), np.log10(10000))
ax2.set_xticks(np.log10(np.array(pos)/10))
ax2.set_xticklabels(['${}$'.format(band).replace('$Ks$', '').replace('K', 'K/Ks') for band in bands])
ax2.set_xlabel('band')
#ax2.set_xscale('log')
line_styles = [':', '-.', '--', '-']
colours = ['y', 'b', 'g', 'r', 'k']
for n, z in enumerate([0.25, 1, 2, 3, 4]):
sed = s.copy()
mod = get_module("redshifting", redshift=z)
mod.process(sed)
ax.plot(sed.wavelength_grid, #*10 factor ten turns angstrom to nm fix unit error
#np.log10(sed.fnu * 1.e-3), #This is for flux in Jy
flux_to_mag(sed.fnu * 1.e-3)[0], #This is for flux in mag
c='k',
#c= colours[n],
#linestyle = line_styles[n],
label= 'z = {}'.format(z),
linewidth=1.0,
alpha=1.
)
for m, band in enumerate(cigale_filternames):
continue
if m == 0:
lab = 'z = {}'.format(z)
else:
lab=None
ax.plot([fwhms[m][0], fwhms[m][1]], [np.log10(gal_fluxes[n, m] )-3,
np.log10(gal_fluxes[n, m] )-3],
# c='k',
c= colours[n],
#linestyle = line_styles[n],
label=lab,
alpha=1.0
)
ax.set_ylabel('5$\sigma$ depth [mag]')
#ax.set_xticks(pos)
#ax.set_xticklabels(['${}$'.format(band.replace('Ks', '').replace('K', 'K/Ks')) for band in bands])
ax.set_xlim(300., 10000)
ax.set_xscale('log')
ax.set_xlabel('Wavelength [nm]')
ax.set_ylim(27, 16)
#ax.legend(loc=2)
#widths
log_widths = np.ones(len(pos)) * (pos) * .2 * widths
#areas/np.max(areas)
parts = ax.violinplot([flux_to_mag(10**d)[0] for d in data], #flux_to_mag assumes Jy
positions=np.array(pos)/10, #fix unit error
widths=log_widths/10,
showmeans=False,
showmedians=False, #widths=widths,
showextrema=False)
for n, part in enumerate(parts['bodies']):
part.set_facecolor(colors[n])
part.set_alpha(1.0)
cax, _ = mpl.colorbar.make_axes([ax, ax2])
n_ticks = 7
values = np.linspace(0,1200, n_ticks)
ticks = values/np.max(areas)
cbar = mpl.colorbar.ColorbarBase(cax, cmap=cmap, ticks = ticks)
cax.set_yticklabels([int(d) for d in values])
cbar.set_label('Area [deg.$^2$]')
plt.figtext(0.49, 0.81, 'z = 0.25')
plt.figtext(0.70, 0.645, 'z = 1')
plt.figtext(0.70, 0.52, 'z = 2')
plt.figtext(0.70, 0.435, 'z = 3')
plt.figtext(0.70, 0.385, 'z = 4')
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)
#fig.suptitle("Violin Plotting Examples")
#fig.subplots_adjust(hspace=0.4)
#plt.ylim(-10,10)
column_width_cm = 8.9
width_cm = 3.0 * column_width_cm
hieght_cm = width_cm / 1.9
width_inches = width_cm/2.5
hieght_inches = hieght_cm/2.5
fig.set_size_inches(width_inches, hieght_inches)
plt.savefig('./figs/band_depths_overviews_areaweighted_with_black_seds_wave_mag.pdf', bbox_inches='tight')
plt.savefig('./figs/band_depths_overviews_areaweighted_with_black_seds_wave_mag.png', bbox_inches='tight')
In [130]:
['${}$'.format(band).replace('$Ks$', '').replace('K', 'K/Ks') for band in bands]
Out[130]:
In [131]:
data[0].shape
Out[131]:
In [132]:
fig, ax = plt.subplots()
#cmap = mpl.cm.get_cmap('viridis', 256)
#norm = mpl.colors.Normalize(vmin=np.min(areas), vmax=np.min(areas))
cmap = mpl.cm.viridis
norm = mpl.colors.Normalize(vmin=np.min(areas), vmax=np.max(areas))
scalmap = mpl.cm.ScalarMappable( cmap=cmap, norm=norm)
colors = scalmap.to_rgba(areas) # The color is the angle
colors[:, 3] = 1
#cb = mpl.colorbar.ColorbarBase(ax, cmap=cmap, norm=norm)
#cmap = mpl.colors.ListedColormap('viridis')
#colors = cmap(np.array(areas))
im = plt.scatter(np.arange(len(areas)), areas, c=colors)
cax, _ = mpl.colorbar.make_axes(ax)
cbar = mpl.colorbar.ColorbarBase(cax, cmap=cmap)
cbar.set_label('area [square degrees]')
#cax.set_yticklabels(['$-\pi$', '0', '$\pi$'])
#cb = mpl.colorbar.ColorbarBase(ax, cmap=cmap, norm=norm)
# Optionally add a colorbar
#cax, _ = mpl.colorbar.make_axes(ax)
#cbar = mpl.colorbar.ColorbarBase(cax, cmap=cmap)
In [133]:
colors
Out[133]:
In [134]:
cmap = mpl.cm.hot
norm = mpl.colors.Normalize(vmin=np.min(areas), vmax=np.max(areas))
scalmap = mpl.cm.ScalarMappable( cmap=cmap, norm=norm)
colors = scalmap.to_rgba(areas) # The color is the angle
#colors[:, 3] = 0.5
print(scalmap.to_rgba(500),
scalmap.to_rgba(1000),
scalmap.to_rgba(100))
colors
Out[134]:
In [135]:
mask = ~np.isnan(coverage['ferr_ap_{}_mean_min'.format('i1')])
mask |= coverage['ferr_ap_{}_mean_min'.format('i1')] > 0.
mask |= coverage['ferr_ap_{}_mean_min'.format('i1')] <1.e3
plt.hist(np.log10(np.array(coverage['ferr_ap_{}_mean_min'.format('i1')][mask ]) *1.e-6), bins=50)
Out[135]:
In [136]:
areas
Out[136]:
In [137]:
np.sum(np.log10(coverage['ferr_ap_{}_mean_min'.format('i1')] ) >3)
Out[137]:
In [138]:
np.sum(coverage['ferr_ap_{}_mean_min'.format('i1')] >10.)
Out[138]:
All field comparison¶
In [139]:
def depths_sample_field(band, field):
mask = ~np.isnan(coverage['ferr_ap_{}_mean_min'.format(band)])
mask &= coverage['ferr_ap_{}_mean_min'.format(band)] > 0.
mask &= coverage['ferr_ap_{}_mean_min'.format(band)] <1.e3
mask &= coverage['field'] == field
area = (np.sum(mask)/len(coverage[coverage['field'] == field])) #* 1270. #Return fractional areas
pixel_depths = np.log10(np.array(coverage['ferr_ap_{}_mean_min'.format(band)][mask]) *5.e-6 )
if np.sum(mask)== 0:
pixel_depths = np.full(100, -9.)
return pixel_depths, area
f = 'ELAIS-N1'
data_field = [ depths_sample_field(band, f)[0] for band in bands ]
areas_field = [ depths_sample_field(band, f)[1] for band in bands ]
In [140]:
areas_field
Out[140]:
In [141]:
dim = [4,6]
fig, axes = plt.subplots(dim[1], dim[0], sharex=True, sharey=True)
plt.rcParams.update({'font.size': 12})
cmap = mpl.cm.viridis
norm = mpl.colors.Normalize(vmin=0, vmax=1)
scalmap = mpl.cm.ScalarMappable( cmap=cmap, norm=norm)
for n, f in enumerate(np.unique(coverage['field'])):
x, y = np.floor_divide(n, dim[0]), np.remainder(n, dim[0])
#axes[x,y].scatter( np.arange(20), np.arange(20))
#axes[x,y].legend()
#####MAKE FIELD DATA#####
data_field = [ depths_sample_field(band, f)[0] for band in bands ]
areas_field = [ depths_sample_field(band, f)[1] for band in bands ]
data_field_mag = [flux_to_mag(10**d)[0] for d in data_field]
colors = scalmap.to_rgba(areas_field) # The color is the angle
colors[:, 3] = 1
def set_axis_style(ax, labels):
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_xticks(np.arange(1, len(labels) + 1))
axes[x,y].set_xticklabels(labels)
axes[x,y].set_xlim(0.25, len(labels) + 0.75)
#axes[x,y].set_xlabel('band')
#axes[x,y].set_ylabel('5$\sigma$ depth [mag]')
#ax.violinplot(np.array(coverage['ferr_ap_{}_mean_min'.format('g') ] ) )
#ax.set_title('Custom violinplot 1', fontsize=fs)
#axes[x,y].set_ylabel('log10( 5$\sigma$ Depths [Jy] )')
set_axis_style(axes[x,y], ['$' + band + '$' for band in bands])
axes[x,y].set_ylim(28, 17)
parts = axes[x,y].violinplot(data_field_mag, widths=0.75, showmeans=False, showmedians=False,
showextrema=False)
#This is just to add the field name
axes[x,y].scatter([-99],[-99],
label=f.replace('SGP', 'HATLAS-SGP').replace('NGP', 'HATLAS-NGP'),
c='w', s=0.0001)
axes[x,y].legend(frameon=False, loc=4)
for n, part in enumerate(parts['bodies']):
part.set_facecolor(colors[n])
part.set_alpha(1)
#cax, _ = mpl.colorbar.make_axes(axes[dim[1]-1,dim[0]-1], panchor=(-8.8, 10.6)) #5 3 depend on 4 x 6
#This bit moves the color bar
cax = inset_axes(axes[dim[1]-1,dim[0]-1],
width="5%", # width = 5% of parent_bbox width
height="85%", # height : 50%
loc=3,
bbox_to_anchor=(0.05, 0.08, 1, 1),
bbox_transform=axes[dim[1]-1,dim[0]-1].transAxes,
borderpad=0,
)
ticks = [0, 0.5, 1]
cbar = mpl.colorbar.ColorbarBase(cax, cmap=cmap, ticks = ticks)
cax.set_yticklabels(ticks)
cbar.set_label('coverage')
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, 'Bands', ha='center')
fig.text(0.04, 0.5, '5$\sigma$ depth [mag]', 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/fields_depths_comparison_grid.pdf', bbox_inches='tight')
plt.savefig('./figs/fields_depths_comparison_grid.png', bbox_inches='tight')
Summarise the filters¶
In [142]:
for band in bands:
specific_bands = [f for f in filters if f.split('_')[1] == band.lower()]
print(specific_bands)
for specific_band in specific_bands:
if 'ferr_ap_{}_mean'.format(specific_band) not in coverage.colnames:
specific_bands.remove(specific_band)
#print(specific_bands)
In [143]:
filter_info = Table.read(herschelhelp_python_loc + 'documentation/filters.csv')
In [144]:
has_broad_typical = np.full(len(filter_info), False)
for band in bands:
has_broad_typical |= filter_info['filter'] == band
filter_info[has_broad_typical]
Out[144]:
In [145]:
has_optical = np.full(len(filter_info), False)
optical_bands = ['u', 'g', 'r', 'i', 'z', 'y']
has_near_infrared = np.full(len(filter_info), False)
near_infrared_bands = ['J', 'H', 'K', 'Ks']
for band in bands:
if band in optical_bands:
has_optical |= filter_info['filter'] == band
if band in near_infrared_bands:
has_near_infrared |= filter_info['filter'] == band
In [146]:
def get_bands(camera):
band_list = []
for b in filter_info['filter'][filter_info['camera'] == camera]:
#print(b)
band_list.append( b)
band_list_ordered_string = ''
for band in bands:
if band in band_list:
band_list_ordered_string += band
return '$' + band_list_ordered_string + '$'
print(get_bands('HyperSuprimeCam'))
In [147]:
optical_sentence = ''
for camera in np.unique(filter_info['camera'][has_optical]):
if camera == 'Standard camera response (camera unknown)':
continue
row = filter_info[filter_info['camera'] == camera][0]
optical_sentence += 'From {} on {} we have the optical bands {}. '.format(row['camera'].replace('The', 'the'),
row['telescope'].replace('The', 'the'),
get_bands(camera))
print(optical_sentence)
In [148]:
near_infrared_sentence = ''
for camera in np.unique(filter_info['camera'][has_near_infrared]):
if camera == 'Standard camera response (camera unknown)':
continue
row = filter_info[filter_info['camera'] == camera][0]
near_infrared_sentence += 'From {} on {} we have the optical bands {}. '.format(row['camera'].replace('The', 'the'),
row['telescope'].replace('The', 'the'),
get_bands(camera))
print(near_infrared_sentence)
Count the number of flagged objects¶
For each band lets count the number of objects that have a flag