CDFS SWIRE master catalogue¶
Preparation of COMBO data¶
Classifying Objects by Medium-Band Observations - a spectrophotometric 17-filter survey (COMBO-17). COMBO catalogue: the catalogue comes from dmu0_COMBO-17.
In the catalogue, we keep:
- The identifier (it's unique in the catalogue);
- The position;
- The stellarity;
- The total magnitude (aperture magnitudes are not provided.
We don't know when the maps have been observed. We will use the year of the reference paper.
In [1]:
from herschelhelp_internal import git_version
print("This notebook was run with herschelhelp_internal version: \n{}".format(git_version()))
In [2]:
%matplotlib inline
#%config InlineBackend.figure_format = 'svg'
import matplotlib.pyplot as plt
plt.rc('figure', figsize=(10, 6))
from collections import OrderedDict
import os
from astropy import units as u
from astropy import constants as const
from astropy.coordinates import SkyCoord
from astropy.table import Column, Table
import numpy as np
from herschelhelp_internal.flagging import gaia_flag_column
from herschelhelp_internal.masterlist import nb_astcor_diag_plot, remove_duplicates
from herschelhelp_internal.utils import astrometric_correction, mag_to_flux, flux_to_mag
In [3]:
OUT_DIR = os.environ.get('TMP_DIR', "./data_tmp")
try:
os.makedirs(OUT_DIR)
except FileExistsError:
pass
RA_COL = "combo_ra"
DEC_COL = "combo_dec"
I - Column selection¶
In [4]:
imported_columns = OrderedDict({
'Seq':'combo_id',
'ra':'combo_ra',
'dec':'combo_dec',
#'dl':'combo_dl',
'stellarity':'combo_stellarity',
'Rmag':'m_combo_r', #The catalogue is R selected
'e_Rmag':'merr_combo_r',
#'Ap_Rmag':'m_ap_combo_r',
#'UjMag':'m_combo_uj', #These bands are derived absolute magnitudes
#'e_UjMag':'merr_combo_uj',
#'BjMag':'m_combo_bj',
#'e_BjMag':'merr_combo_bj',
#'VjMag':'m_combo_vj',
#'e_VjMag':'merr_combo_vj',
#'usMag':'m_combo_us',
#'e_usMag':'merr_combo_us',
#'gsMag':'m_combo_gs',
#'e_gsMag':'merr_combo_gs',
#'rsMag':'m_combo_rs',
#'e_rsMag':'merr_combo_rs',
#'UbMag':'m_combo_ub',
#'e_UbMag':'merr_combo_ub',
#'BbMag':'m_combo_bb',
#'e_BbMag':'merr_combo_bb',
#'VbMag':'m_combo_vb',
#'e_VbMag':'merr_combo_vb',
#'S280Mag':'m_combo_s280',
#'e_S280Mag':'merr_combo_s280',
##'S145Mag':'m_combo_s145',
#'e_S145Mag':'merr_combo_s145',
'W420F_E':'f_ap_combo_420', #The following values are given as fluxes
'e_W420F_E':'ferr_ap_combo_420',
'W462F_E':'f_ap_combo_462',
'e_W462F_E':'ferr_ap_combo_462',
'W485F_D':'f_ap_combo_485',
'e_W485F_D':'ferr_ap_combo_485',
'W518F_E':'f_ap_combo_518',
'e_W518F_E':'ferr_ap_combo_518',
'W571F_S':'f_ap_combo_571', #Combined flux from two runs
'e_W571F_S':'ferr_ap_combo_571',
'W604F_E':'f_ap_combo_604',
'e_W604F_E':'ferr_ap_combo_604',
'W646F_D':'f_ap_combo_646',
'e_W646F_D':'ferr_ap_combo_646',
'W696F_E':'f_ap_combo_696',
'e_W696F_E':'ferr_ap_combo_696',
'W753F_E':'f_ap_combo_753',
'e_W753F_E':'ferr_ap_combo_753',
'W815F_S':'f_ap_combo_815',
'e_W815F_S':'ferr_ap_combo_815',
'W856F_D':'f_ap_combo_856',
'e_W856F_D':'ferr_ap_combo_856',
'W914F_D':'f_ap_combo_914', #Two runs but no combined - taking first
'e_W914F_D':'ferr_ap_combo_914',
'UF_S':'f_ap_combo_u',
'e_UF_S':'ferr_ap_combo_u',
'BF_S':'f_ap_combo_b',
'e_BF_S':'ferr_ap_combo_b',
'VF_D':'f_ap_combo_v',
'e_VF_D':'ferr_ap_combo_v',
'RF_S':'f_ap_combo_r',
'e_RF_S':'ferr_ap_combo_r',
'IF_D':'f_ap_combo_i',
'e_IF_D':'ferr_ap_combo_i'
})
catalogue = Table.read("../../dmu0/dmu0_COMBO-17/data/table3.fits")[list(imported_columns)]
for column in imported_columns:
catalogue[column].name = imported_columns[column]
epoch = 2000 #table says 1999 to 2001
# Clean table metadata
catalogue.meta = None
In [5]:
catalogue[:10].show_in_notebook()
Out[5]:
Unit conversion¶
Using mid wavelength¶
The flux is presented in $\textrm{photons} .\textrm{s}^{-1} . \textrm{m}^{-2} .\textrm{nm}^{-1}$. We wish to convert these to micro Jansky; $10^{-32} \textrm{ W} . \textrm{m}^{-2} . \textrm{Hz}^{-1} $.
To convert $\textrm{photons} . \textrm{s}^{-1} $ to $\textrm{W}$ we must multiply by the average photon energy $h c / \lambda$. We presume that the COMBO mid point wavelength was used.
To convert $\textrm{nm}^{-1}$ to $\textrm{Hz}^{-1}$ we must differenciate:
$c = \nu \lambda $
$\nu = c / \lambda$
$\frac{d \nu}{d \lambda} = - c /\lambda^{2}$
$d \lambda = - (\lambda^{2} / c )\times d \nu$
The net result of this is to multiply by $\lambda^2 / c$.
Combining these two unit conversions leads to overall multiplying by $ h \lambda$:
$(\lambda^2 / c ) \times (h c / \lambda) = h \lambda$
In [6]:
#Example conversion from photon s^-1 m^-2 nm^-1 to Jy
flux_lambda = 0.0064185 * (u.m **-2) * (u.s ** -1) *( u.nm ** -1)
wavelength = 420 * u.nm
flux_nu = flux_lambda * const.h * wavelength
print('flux_lambda:', flux_lambda)
#print('f_lambda:', flux_lambda.decompose())
print('flux_nu:', flux_nu)
print('flux_nu in Jy:',flux_nu.to(u.Jy))
Using Vega flux¶
The paper provides fluxes of Vega in each band allowing conversion to Vega mag. It then provides the difference between Vega mag and AB mag per band to convert to AB. We can therefore calculate the AB magnitude and convert that back to flux to get the flux in Jy.
In [7]:
#Set wavelengths for unit conversion.
#All from http://cds.aanda.org/component/article?access=bibcode&bibcode=2004A%252526A...421..913W
#Wavelengths in the table column headings are marginally different to those in the paper.
# Band name \lambda Vega AB Fphot (Vega 10^8 photons s^-1 m^-2 nm^-1)
wavelengths = {
'f_ap_combo_420': [418, -0.19, 1.571],
'ferr_ap_combo_420':[418, -0.19, 1.571],
'f_ap_combo_462': [462, -0.18, 1.412],
'ferr_ap_combo_462':[462, -0.18, 1.412],
'f_ap_combo_485': [486, -0.06, 1.207],
'ferr_ap_combo_485':[486, -0.06, 1.207],
'f_ap_combo_518': [519, -0.06, 1.125],
'ferr_ap_combo_518':[519, -0.06, 1.125],
'f_ap_combo_571': [572, 0.04, 0.932],
'ferr_ap_combo_571':[572, 0.04, 0.932],
'f_ap_combo_604': [605, 0.10, 0.832],
'ferr_ap_combo_604':[605, 0.10, 0.832],
'f_ap_combo_646': [645, 0.22, 0.703],
'ferr_ap_combo_646':[645, 0.22, 0.703],
'f_ap_combo_696': [696, 0.27, 0.621],
'ferr_ap_combo_696':[696, 0.27, 0.621],
'f_ap_combo_753': [753, 0.36, 0.525],
'ferr_ap_combo_753':[753, 0.36, 0.525],
'f_ap_combo_815': [816, 0.45, 0.442],
'ferr_ap_combo_815':[816, 0.45, 0.442],
'f_ap_combo_856': [857, 0.56, 0.386],
'ferr_ap_combo_856':[857, 0.56, 0.386],
'f_ap_combo_914': [914, 0.50, 0.380],
'ferr_ap_combo_914':[914, 0.50, 0.380],
'f_ap_combo_u': [365, 0.77, 0.737],
'ferr_ap_combo_u': [365, 0.77, 0.737],
'f_ap_combo_b': [458, -0.13, 1.371],
'ferr_ap_combo_b': [458, -0.13, 1.371],
'f_ap_combo_v': [538, -0.02, 1.055],
'ferr_ap_combo_v': [538, -0.02, 1.055],
'f_ap_combo_r': [648, 0.19, 0.725],
'ferr_ap_combo_r': [648, 0.19, 0.725],
'f_ap_combo_i': [857, 0.49, 0.412],
'ferr_ap_combo_i': [857, 0.49, 0.412]
}
In [8]:
#Example conversion from photon s^-1 m^-2 nm^-1 to Jy
flux_lambda = 0.0064185 #* (u.m **-2) * (u.s ** -1) *( u.nm ** -1)
mag_vega = -2.5 *np.log10(flux_lambda/(wavelengths['f_ap_combo_420'][2]*1.e8))
print(mag_vega)
mag_AB = mag_vega + wavelengths['f_ap_combo_420'][1]
print(mag_AB)
flux_converted = mag_to_flux(mag_AB)
print(flux_converted)
#This is different to the value calculated using mid point wavelength as we expect
#because it should take account of the filter response better.
In [9]:
#Replace 0.0 with NaN values
for col in catalogue.colnames:
catalogue[col].unit = None
if col.startswith('m'): # | col.endswith('ra') | col.endswith('dec'):
catalogue[col][np.where(catalogue[col] == 0.0)] = np.nan
# Add magnitude, fix flux units and add band-flag columns
nancol = np.zeros(len(catalogue))
nancol.fill(np.nan)
for col in catalogue.colnames:
if col.startswith('f_'):
errcol = "ferr{}".format(col[1:])
#Replace mask with nan
catalogue[col].fill_value = np.nan
catalogue[errcol].fill_value = np.nan
#Calculate magnitudes using Vega fluxes from wavelengths (values from paper)
magnitude = -2.5*np.log10(catalogue[col].filled()/(wavelengths[col][2]*1.e8)) + wavelengths[col][1]
magnitude_error = 2.5/np.log(10)*(catalogue[errcol].filled()/catalogue[col].filled() )
catalogue.add_column(Column(nancol,
name="m{}".format(col[1:])))
catalogue.add_column(Column(nancol,
name="m{}".format(errcol[1:])))
# Add the AB magnitudes
catalogue["m{}".format(col[1:])] = magnitude
catalogue["m{}".format(errcol[1:])] = magnitude_error
flux_new, flux_new_error = mag_to_flux(magnitude, magnitude_error)
catalogue[col] = flux_new * 1.e6 # uJy
catalogue[errcol] = flux_new_error * 1.e6 # uJy
#We add NAN filled total columns because no total fluxes are present
if not col == 'f_ap_combo_r':
catalogue.add_column(Column(nancol,
name="m{}".format(col[4:])))
catalogue.add_column(Column(nancol,
name="merr{}".format(col[4:])))
catalogue.add_column(Column(nancol,
name="f{}".format(col[4:])))
catalogue.add_column(Column(nancol,
name="ferr{}".format(col[4:])))
#
# Band-flag column
catalogue.add_column(Column(np.zeros(len(catalogue), dtype=bool), name="flag{}".format(col[4:])))
#Add total fluxes for R band
f_combo_r, ferr_combo_r = mag_to_flux(catalogue['m_combo_r'],catalogue['merr_combo_r'])
catalogue['f_combo_r'] = f_combo_r
catalogue['ferr_combo_r'] = ferr_combo_r
In [10]:
catalogue[:10].show_in_notebook()
Out[10]:
II - Removal of duplicated sources¶
We remove duplicated objects from the input catalogues.
In [11]:
SORT_COLS = ['merr_combo_u',
'merr_combo_b',
'merr_combo_v',
'merr_combo_r',
'merr_combo_i']
FLAG_NAME = 'combo_flag_cleaned'
nb_orig_sources = len(catalogue)
catalogue = remove_duplicates(catalogue, RA_COL, DEC_COL, sort_col=SORT_COLS,flag_name=FLAG_NAME)
nb_sources = len(catalogue)
print("The initial catalogue had {} sources.".format(nb_orig_sources))
print("The cleaned catalogue has {} sources ({} removed).".format(nb_sources, nb_orig_sources - nb_sources))
print("The cleaned catalogue has {} sources flagged as having been cleaned".format(np.sum(catalogue[FLAG_NAME])))
III - Astrometry correction¶
We match the astrometry to the Gaia one. We limit the Gaia catalogue to sources with a g band flux between the 30th and the 70th percentile. Some quick tests show that this give the lower dispersion in the results.
In [12]:
gaia = Table.read("../../dmu0/dmu0_GAIA/data/GAIA_CDFS-SWIRE.fits")
gaia_coords = SkyCoord(gaia['ra'], gaia['dec'])
In [13]:
nb_astcor_diag_plot(catalogue[RA_COL], catalogue[DEC_COL],
gaia_coords.ra, gaia_coords.dec)
In [14]:
delta_ra, delta_dec = astrometric_correction(
SkyCoord(catalogue[RA_COL], catalogue[DEC_COL]),
gaia_coords
)
print("RA correction: {}".format(delta_ra))
print("Dec correction: {}".format(delta_dec))
In [15]:
catalogue[RA_COL] += delta_ra.to(u.deg)
catalogue[DEC_COL] += delta_dec.to(u.deg)
In [16]:
nb_astcor_diag_plot(catalogue[RA_COL], catalogue[DEC_COL],
gaia_coords.ra, gaia_coords.dec)
IV - Flagging Gaia objects¶
In [17]:
catalogue.add_column(
gaia_flag_column(SkyCoord(catalogue[RA_COL], catalogue[DEC_COL]), epoch, gaia)
)
In [18]:
GAIA_FLAG_NAME = "combo_flag_gaia"
catalogue['flag_gaia'].name = GAIA_FLAG_NAME
print("{} sources flagged.".format(np.sum(catalogue[GAIA_FLAG_NAME] > 0)))
V - Flagging objects near bright stars¶
VI - Saving to disk¶
In [19]:
catalogue.write("{}/COMBO.fits".format(OUT_DIR), overwrite=True)