from herschelhelp_internal import git_version print("This notebook was run with herschelhelp_internal version: \n{}".format(git_version()))
%matplotlib inline
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.coordinates import SkyCoord from astropy.table import Column, Table, join 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, flux_to_mag
OUT_DIR = os.environ.get('TMP_DIR', "./data_tmp") try: os.makedirs(OUT_DIR) except FileExistsError: pass
RA_COL = "acs_ra" DEC_COL = "acs_dec"
b_catalogue = Table.read("../../dmu0/dmu0_GOODS-ACS/data/h_goods_nb_r2.0z_cat.fits") #[list(imported_columns)] i_catalogue = Table.read("../../dmu0/dmu0_GOODS-ACS/data/h_goods_ni_r2.0z_cat.fits") #[list(imported_columns)] v_catalogue = Table.read("../../dmu0/dmu0_GOODS-ACS/data/h_goods_nv_r2.0z_cat.fits") #[list(imported_columns)] z_catalogue = Table.read("../../dmu0/dmu0_GOODS-ACS/data/h_goods_nz_r2.0z_cat.fits") #[list(imported_columns)]
for column in b_catalogue.colnames: if not b_catalogue[column].name == 'ID_IAU': bcatalogue[column].name = 'b' + b_catalogue[column].name
for column in i_catalogue.colnames: if not i_catalogue[column].name == 'ID_IAU': icatalogue[column].name = 'i' + i_catalogue[column].name
for column in v_catalogue.colnames: if not v_catalogue[column].name == 'ID_IAU': vcatalogue[column].name = 'v' + v_catalogue[column].name
for column in z_catalogue.colnames: if not z_catalogue[column].name == 'ID_IAU': zcatalogue[column].name = 'z' + z_catalogue[column].name
epoch = 2012 #Year of publication
catalogue = b_catalogue
catalogue = join(b_catalogue,i_catalogue, keys='ID_IAU') catalogue = join(catalogue,v_catalogue, keys='ID_IAU') catalogue = join(catalogue,z_catalogue, keys='ID_IAU')
catalogue[:10].show_in_notebook()
imported_columns = OrderedDict({ 'ID_IAU': "acs_id", 'b_ALPHA_J2000': "acs_ra", 'b_DELTA_J2000': "acs_dec", 'b_MAG_AUTO': "m_acs_b", 'b_MAGERR_AUTO': "merr_acs_b", 'b_FLUX_AUTO': "f_acs_b", 'b_FLUXERR_AUTO': "ferr_acs_b", 'e_F140W': "ferr_candels_f140w",
})
catalogue = catalogue[list(imported_columns)]
for column in imported_columns: catalogue[column].name = imported_columns[column]
"Adding flux and band-flag columns" for col in catalogue.colnames: if col.startswith('f_'):
errcol = "ferr{}".format(col[1:])
#Calculate mags, errors including the fact that fluxes are in units of 0.3631 uJy
mag, error = flux_to_mag(np.array(catalogue[col]) * 0.3631e-6, np.array(catalogue[errcol] * 0.3631e-6))
# magnitudes are added
catalogue.add_column(Column(mag, name="m{}".format(col[1:])))
catalogue.add_column(Column(error, name="m{}".format(errcol[1:])))
#Correct flux units to uJy
catalogue[col] = catalogue[col] * 0.3631
catalogue[col].unit = u.microjansky
catalogue[errcol] = catalogue[errcol] * 0.3631
catalogue[errcol].unit = u.microjansky
if ('125' in col) or ('814' in col) or ('606' in col) :
# We add nan filled aperture photometry for consistency
catalogue.add_column(Column(np.full(len(catalogue), np.nan), name="m_ap{}".format(col[1:])))
catalogue.add_column(Column(np.full(len(catalogue), np.nan), name="merr_ap{}".format(col[1:])))
catalogue.add_column(Column(np.full(len(catalogue), np.nan), name="f_ap{}".format(col[1:])))
catalogue.add_column(Column(np.full(len(catalogue), np.nan), name="ferr_ap{}".format(col[1:])))
# Band-flag column
if "ap" not in col:
catalogue.add_column(Column(np.zeros(len(catalogue), dtype=bool), name="flag{}".format(col[1:])))
TODO: Set to True the flag columns for fluxes that should not be used for SED fitting.
catalogue[:10].show_in_notebook()
We remove duplicated objects from the input catalogues.
SORT_COLS = ['merr_candels_f140w', 'merr_candels_f160w', 'merr_candels_f606w', 'merr_candels_f814w', 'merr_candels_f125w'] FLAG_NAME = 'candels_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])))
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.
gaia = Table.read("../../dmu0/dmu0_GAIA/data/GAIA_XMM-LSS.fits") gaia_coords = SkyCoord(gaia['ra'], gaia['dec'])
nb_astcor_diag_plot(catalogue[RA_COL], catalogue[DEC_COL], gaia_coords.ra, gaia_coords.dec)
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))
catalogue[RA_COL] += delta_ra.to(u.deg) catalogue[DEC_COL] += delta_dec.to(u.deg)
nb_astcor_diag_plot(catalogue[RA_COL], catalogue[DEC_COL], gaia_coords.ra, gaia_coords.dec)
catalogue.add_column( gaia_flag_column(SkyCoord(catalogue[RA_COL], catalogue[DEC_COL]), epoch, gaia) )
GAIA_FLAG_NAME = "candels_flag_gaia"
catalogue['flag_gaia'].name = GAIA_FLAG_NAME print("{} sources flagged.".format(np.sum(catalogue[GAIA_FLAG_NAME] > 0)))
catalogue.write("{}/CANDELS.fits".format(OUT_DIR), overwrite=True)