This catalogue comes from dmu0_SpARCS
. Alexandru Tudorica confirmed that the magnitudes are AB ones and are not aperture corrected.
In the catalogue, we keep:
Is there y band data?
The maps on the web page indicate they were observed in 2012 (or late 2011). Let's use 2012 as epoch.
from herschelhelp_internal import git_version
print("This notebook was run with herschelhelp_internal version: \n{}".format(git_version()))
%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
import numpy as np
from herschelhelp_internal.flagging import gaia_flag_column
from herschelhelp_internal.masterlist import nb_astcor_diag_plot, remove_duplicates, nb_plot_mag_ap_evol, nb_plot_mag_vs_apcor
from herschelhelp_internal.utils import aperture_correction, astrometric_correction, mag_to_flux
OUT_DIR = os.environ.get('TMP_DIR', "./data_tmp")
try:
os.makedirs(OUT_DIR)
except FileExistsError:
pass
RA_COL = "sparcs_ra"
DEC_COL = "sparcs_dec"
To compute aperture correction we need to dertermine two parametres: the target aperture and the range of magnitudes for the stars that will be used to compute the correction.
Target aperture: To determine the target aperture, we simulate a curve of growth using the provided apertures and draw two figures:
As target aperture, we should use the smallest (i.e. less noisy) aperture for which most of the flux is captures.
Magnitude range: To know what limits in aperture to use when doing the aperture correction, we plot for each magnitude bin the correction that is computed and its RMS. We should then use the wide limits (to use more stars) where the correction is stable and with few dispersion.
# We are using the aperture index 7 (0 base) that correspond to 11 pix * 0.186 arcsec/pix = 2.046 arcsec
AP_INDEX = 7
t = Table.read("../../dmu0/dmu0_SpARCS/data/SpARCS_HELP_ELAIS-N2.fits")
stellarity = t['CLASS_STAR']
mags_r = np.array(t['MAG_APER_r']).T
mags_r[mags_r == 99] = np.nan
mags_u = np.array(t['MAG_APER_u']).T
mags_u[mags_u == 99] = np.nan
mags_g = np.array(t['MAG_APER_g']).T
mags_g[mags_g == 99] = np.nan
mags_z = np.array(t['MAG_APER_z']).T
mags_z[mags_z == 99] = np.nan
del t
nb_plot_mag_ap_evol(mags_r, stellarity)
We will use the 16th (aperture number above begin to 0) aperture as target.
nb_plot_mag_vs_apcor(mags_r[AP_INDEX], mags_r[15], stellarity)
We use magnitudes between 17 and 17.9.
nb_plot_mag_ap_evol(mags_u, stellarity)
We will use the 16th (aperture number above begin to 0) aperture as target. Should we use the 12nd because of the increasing magnitude?
nb_plot_mag_vs_apcor(mags_u[AP_INDEX], mags_u[15], stellarity)
We use magnitudes between 17 and 17.9.
nb_plot_mag_ap_evol(mags_g, stellarity)
We will use the 16th (aperture number above begin to 0) aperture as target.
nb_plot_mag_vs_apcor(mags_g[AP_INDEX], mags_g[15], stellarity)
We use magnitudes between 17.2 and 18.
nb_plot_mag_ap_evol(mags_z, stellarity)
We will use the 16th (aperture number above begin to 0) aperture as target.
nb_plot_mag_vs_apcor(mags_z[AP_INDEX], mags_z[15], stellarity)
We use magnitudes between 16 and 17.
# Index of the target aperture when doing aperture correction
# (see 'sparcs_aperture_correction' notebook).
AP_TARG_INDEX = {
'u': 15,
'g': 15,
'r': 15,
'z': 15
}
# Magnitude range for aperture correction (see 'sparcs_aperture_correction' notebook).
APCOR_MAG_LIMITS = {
'u': (17., 17.9),
'g': (17., 17.9),
'r': (17.2, 18.),
'z': (16., 17.)
}
epoch = 2012
sparcs_tmp = Table.read("../../dmu0/dmu0_SpARCS/data/SpARCS_HELP_ELAIS-N2.fits")[
'internal_id', 'ALPHA_J2000', 'DELTA_J2000', 'CLASS_STAR',
'MAG_APER_r', 'MAGERR_APER_r', 'MAG_AUTO_r', 'MAGERR_AUTO_r',
'MAG_APER_u', 'MAGERR_APER_u', 'MAG_AUTO_u', 'MAGERR_AUTO_u',
'MAG_APER_g', 'MAGERR_APER_g', 'MAG_AUTO_g', 'MAGERR_AUTO_g',
'MAG_APER_z', 'MAGERR_APER_z', 'MAG_AUTO_z', 'MAGERR_AUTO_z']
# Clean table metadata
sparcs_tmp.meta = None
catalogue = Table(
data = sparcs_tmp['internal_id', 'ALPHA_J2000', 'DELTA_J2000', 'CLASS_STAR'],
names = ['sparcs_intid', 'sparcs_ra', 'sparcs_dec', 'sparcs_stellarity'])
for band in ['u', 'g', 'r', 'z']:
# Aperture magnitudes
mag_aper = sparcs_tmp["MAG_APER_{}".format(band)][:, AP_INDEX]
mag_aper_target = sparcs_tmp["MAG_APER_{}".format(band)][:, AP_TARG_INDEX[band]]
magerr_aper = sparcs_tmp["MAGERR_APER_{}".format(band)][:, AP_INDEX]
# Set bad values (99.0) to NaN
mask = (mag_aper > 90) | (mag_aper_target > 90) | (magerr_aper > 90)
mag_aper[mask] = np.nan
mag_aper_target[mask] = np.nan
magerr_aper[mask] = np.nan
# Aperture correction
mag_diff, num, std = aperture_correction(
mag_aper, mag_aper_target, catalogue['sparcs_stellarity'],
mag_min=APCOR_MAG_LIMITS[band][0], mag_max=APCOR_MAG_LIMITS[band][1]
)
print("Aperture correction for SpARCS band {}:".format(band))
print("Correction: {}".format(mag_diff))
print("Number of source used: {}".format(num))
print("RMS: {}".format(std))
print("")
mag_aper += mag_diff
catalogue.add_column(Column(
data = mag_aper.data,
name = "m_ap_cfht_megacam_{}".format(band)
))
catalogue.add_column(Column(
data = magerr_aper.data,
name = "merr_ap_cfht_megacam_{}".format(band)
))
# Computing the aperture flux columns
flux_aper, fluxerr_aper = mag_to_flux(mag_aper.data, magerr_aper.data)
catalogue.add_column(Column(
data = flux_aper * 1.e6,
name = "f_ap_cfht_megacam_{}".format(band)
))
catalogue.add_column(Column(
data = fluxerr_aper * 1.e6,
name = "ferr_ap_cfht_megacam_{}".format(band)
))
# Auto magnitudes
mag_auto = sparcs_tmp["MAG_AUTO_{}".format(band)]
magerr_auto = sparcs_tmp["MAGERR_AUTO_{}".format(band)]
# Set bad values (99.0) to NaN
mask = (mag_auto > 90) | (magerr_auto > 90)
mag_auto[mask] = np.nan
magerr_auto[mask] = np.nan
catalogue.add_column(Column(
data = mag_auto,
name = "m_cfht_megacam_{}".format(band)
))
catalogue.add_column(Column(
data = magerr_auto,
name = "merr_cfht_megacam_{}".format(band)
))
# Computing the flux columns
flux, fluxerr = mag_to_flux(mag_auto, magerr_auto)
catalogue.add_column(Column(
data = flux * 1.e6,
name = "f_cfht_megacam_{}".format(band)
))
catalogue.add_column(Column(
data = fluxerr * 1.e6,
name = "ferr_cfht_megacam_{}".format(band)
))
# Band-flag column
catalogue.add_column(Column(np.zeros(len(catalogue), dtype=bool),
name="flag_cfht_megacam_{}".format(band)
))
# TODO: Set to True the flag columns for fluxes that should not be used for SED fitting.
del sparcs_tmp
catalogue[:10].show_in_notebook()
We remove duplicated objects from the input catalogues.
SORT_COLS = ['merr_ap_cfht_megacam_r', 'merr_ap_cfht_megacam_u',
'merr_ap_cfht_megacam_g', 'merr_ap_cfht_megacam_z']
FLAG_NAME = 'sparcs_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_ELAIS-N2.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 = "sparcs_flag_gaia"
catalogue['flag_gaia'].name = GAIA_FLAG_NAME
print("{} sources flagged.".format(np.sum(catalogue[GAIA_FLAG_NAME] > 0)))
catalogue.write("{}/SpARCS.fits".format(OUT_DIR), overwrite=True)