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.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, mag_to_flux
from herschelhelp_internal.masterlist import merge_catalogues, nb_merge_dist_plot
In [3]:
OUT_DIR = os.environ.get('TMP_DIR', "./data_tmp")
try:
os.makedirs(OUT_DIR)
except FileExistsError:
pass
RA_COL = "wirds_ra"
DEC_COL = "wirds_dec"
I - Column selection¶
In [4]:
#We have to import and combine the H, J and Ks catalogues separately.
#Fluxes are given in counts sowe compute them fresh from the magnitudes
epoch = 2007
imported_columns_h = OrderedDict({
'NUMBER': "wirds_h_id",
'ALPHA_J2000': "wirds_h_ra",
'DELTA_J2000': "wirds_h_dec",
'CLASS_STAR': "wirds_h_stellarity",
'MAG_APER': "m_ap_wirds_h",
'MAGERR_APER': "merr_ap_wirds_h",
'MAG_AUTO': "m_wirds_h",
'MAGERR_AUTO': "merr_wirds_h"
#'FLUX_APER': "f_ap_wirds_h",
#'FLUXERR_APER': "ferr_ap_wirds_h",
#'FLUX_AUTO': "f_wirds_h",
#'FLUXERR_AUTO': "ferr_wirds_h"
})
catalogue_h = Table.read("../../dmu0/dmu0_CFHT-WIRDS/data/COSMOS_H.fits")[list(imported_columns_h)]
for column in imported_columns_h:
catalogue_h[column].name = imported_columns_h[column]
for col in catalogue_h.colnames:
if col.startswith('m_'):
errcol = "merr{}".format(col[1:])
#catalogue_h[col].name = imported_columns_h[col]
flux, error = mag_to_flux(np.array(catalogue_h[col]), np.array(catalogue_h[errcol]))
# Fluxes are added in µJy
catalogue_h.add_column(Column(flux * 1.e6, name="f{}".format(col[1:])))
catalogue_h.add_column(Column(error * 1.e6, name="f{}".format(errcol[1:])))
# Band-flag column
if "ap" not in col:
catalogue_h.add_column(Column(np.zeros(len(catalogue_h), dtype=bool), name="flag{}".format(col[1:])))
# Clean table metadata
catalogue_h.meta = None
imported_columns_j = OrderedDict({
'NUMBER': "wirds_j_id",
'ALPHA_J2000': "wirds_j_ra",
'DELTA_J2000': "wirds_j_dec",
'CLASS_STAR': "wirds_j_stellarity",
'MAG_APER': "m_ap_wirds_j",
'MAGERR_APER': "merr_ap_wirds_j",
'MAG_AUTO': "m_wirds_j",
'MAGERR_AUTO': "merr_wirds_j"
#'FLUX_APER': "f_ap_wirds_j",
#'FLUXERR_APER': "ferr_ap_wirds_j",
#'FLUX_AUTO': "f_wirds_j",
#'FLUXERR_AUTO': "ferr_wirds_j"
})
catalogue_j = Table.read("../../dmu0/dmu0_CFHT-WIRDS/data/COSMOS_J.fits")[list(imported_columns_j)]
for column in imported_columns_j:
catalogue_j[column].name = imported_columns_j[column]
for col in catalogue_j.colnames:
if col.startswith('m_'):
errcol = "merr{}".format(col[1:])
#catalogue_h[col].name = imported_columns_h[col]
flux, error = mag_to_flux(np.array(catalogue_j[col]), np.array(catalogue_j[errcol]))
# Fluxes are added in µJy
catalogue_j.add_column(Column(flux * 1.e6, name="f{}".format(col[1:])))
catalogue_j.add_column(Column(error * 1.e6, name="f{}".format(errcol[1:])))
# Band-flag column
if "ap" not in col:
catalogue_j.add_column(Column(np.zeros(len(catalogue_j), dtype=bool), name="flag{}".format(col[1:])))
# Clean table metadata
catalogue_j.meta = None
imported_columns_ks = OrderedDict({
'NUMBER': "wirds_ks_id",
'ALPHA_J2000': "wirds_ks_ra",
'DELTA_J2000': "wirds_ks_dec",
'CLASS_STAR': "wirds_ks_stellarity",
'MAG_APER': "m_ap_wirds_ks",
'MAGERR_APER': "merr_ap_wirds_ks",
'MAG_AUTO': "m_wirds_ks",
'MAGERR_AUTO': "merr_wirds_ks"
#'FLUX_APER': "f_ap_wirds_ks",
#'FLUXERR_APER': "ferr_ap_wirds_ks",
#'FLUX_AUTO': "f_wirds_ks",
#'FLUXERR_AUTO': "ferr_wirds_ks"
})
catalogue_ks = Table.read("../../dmu0/dmu0_CFHT-WIRDS/data/COSMOS_Ks.fits")[list(imported_columns_ks)]
for column in imported_columns_ks:
catalogue_ks[column].name = imported_columns_ks[column]
for col in catalogue_ks.colnames:
if col.startswith('m_'):
errcol = "merr{}".format(col[1:])
#catalogue_h[col].name = imported_columns_h[col]
flux, error = mag_to_flux(np.array(catalogue_ks[col]), np.array(catalogue_ks[errcol]))
# Fluxes are added in µJy
catalogue_ks.add_column(Column(flux * 1.e6, name="f{}".format(col[1:])))
catalogue_ks.add_column(Column(error * 1.e6, name="f{}".format(errcol[1:])))
# Band-flag column
if "ap" not in col:
catalogue_ks.add_column(Column(np.zeros(len(catalogue_ks), dtype=bool), name="flag{}".format(col[1:])))
# Clean table metadata
catalogue_ks.meta = None
Merging different bands¶
CFHT-WIRDS has indivdual extractions from each band. We must therefore merge them as if they were individual catalogues (they have different
In [5]:
catalogue_ks[:10].show_in_notebook()
Out[5]:
In [6]:
nb_merge_dist_plot(
SkyCoord(catalogue_ks['wirds_ks_ra'], catalogue_ks['wirds_ks_dec']),
SkyCoord(catalogue_j['wirds_j_ra'], catalogue_j['wirds_j_dec'])
)
In [7]:
catalogue = catalogue_ks
catalogue_ks['wirds_ks_ra'].name = 'ra'
catalogue_ks['wirds_ks_dec'].name = 'dec'
# Given the graph above, we use 0.8 arc-second radius
catalogue = merge_catalogues(catalogue, catalogue_j, "wirds_j_ra", "wirds_j_dec", radius=0.8*u.arcsec)
In [8]:
nb_merge_dist_plot(
SkyCoord(catalogue['ra'], catalogue['dec']),
SkyCoord(catalogue_h['wirds_h_ra'], catalogue_h['wirds_h_dec'])
)
In [9]:
# Given the graph above, we use 0.8 arc-second radius
catalogue = merge_catalogues(catalogue, catalogue_h, "wirds_h_ra", "wirds_h_dec", radius=0.8*u.arcsec)
In [10]:
#rename radec colums
catalogue['ra'].name = 'wirds_ra'
catalogue['dec'].name = 'wirds_dec'
In [11]:
for col in catalogue.colnames:
if "m_" in col or "merr_" in col or "f_" in col or "ferr_" in col or "stellarity" in col:
catalogue[col].fill_value = np.nan
elif "flag" in col:
catalogue[col].fill_value = 0
elif "id" in col:
catalogue[col].fill_value = 999999
catalogue = catalogue.filled()
Generate internal id¶
Since every source has an independent id we combine them in 6 digit groups so that each iondividual id can be retrieved from the final integer
In [12]:
wirds_intid = catalogue['wirds_ks_id'] + catalogue['wirds_j_id']*10**6 + catalogue['wirds_h_id']*10**12
In [13]:
catalogue.add_column(Column(data=wirds_intid, name="wirds_intid"))
catalogue.remove_columns(['wirds_ks_id','wirds_j_id','wirds_h_id'])
In [14]:
stellarity_columns = [column for column in catalogue.colnames
if 'stellarity' in column]
print(", ".join(stellarity_columns))
In [15]:
# We create an masked array with all the stellarities and get the maximum value, as well as its
# origin. Some sources may not have an associated stellarity.
stellarity_array = np.array([catalogue[column] for column in stellarity_columns])
stellarity_array = np.ma.masked_array(stellarity_array, np.isnan(stellarity_array))
max_stellarity = np.max(stellarity_array, axis=0)
max_stellarity.fill_value = np.nan
catalogue.add_column(Column(data=max_stellarity.filled(), name="wirds_stellarity"))
catalogue.remove_columns(stellarity_columns)
catalogue['flag_merged'].name = 'wirds_flag_merged'
In [16]:
catalogue[:10].show_in_notebook()
Out[16]:
II - Removal of duplicated sources¶
We remove duplicated objects from the input catalogues.
In [17]:
SORT_COLS = ['merr_ap_wirds_ks',
'merr_ap_wirds_j',
'merr_ap_wirds_h']
FLAG_NAME = 'wirds_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 [18]:
gaia = Table.read("../../dmu0/dmu0_GAIA/data/GAIA_COSMOS.fits")
gaia_coords = SkyCoord(gaia['ra'], gaia['dec'])
In [19]:
nb_astcor_diag_plot(catalogue[RA_COL], catalogue[DEC_COL],
gaia_coords.ra, gaia_coords.dec)
In [20]:
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 [21]:
catalogue[RA_COL] = catalogue[RA_COL] + delta_ra.to(u.deg)
catalogue[DEC_COL] = catalogue[DEC_COL] + delta_dec.to(u.deg)
In [22]:
nb_astcor_diag_plot(catalogue[RA_COL], catalogue[DEC_COL],
gaia_coords.ra, gaia_coords.dec)
IV - Flagging Gaia objects¶
In [23]:
catalogue.add_column(
gaia_flag_column(SkyCoord(catalogue[RA_COL], catalogue[DEC_COL]), epoch, gaia)
)
In [24]:
GAIA_FLAG_NAME = "cfht-wirds_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 [25]:
catalogue.write("{}/CFHT-WIRDS.fits".format(OUT_DIR), overwrite=True)