The catalogue comes from dmu0_ESIS-VOICE
.
In the catalogue, we keep:
We don't know when the maps have been observed. We will use the year of the reference paper.
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 = "voice_ra"
DEC_COL = "voice_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 4 (1 base) that correspond 2 arcsec
# Apertures 1 to 20 correspond to these diameters (in arc-seconds): 0.5, 1.0,
# 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 15.0, 20.0,
# 30.0, 40.0, and 50.0.
AP_INDEX = 4 - 1 #Our aperture nos start at 0
t = Table.read("../../dmu0/dmu0_ESIS-VOICE/data/esis_b2vr_cat_03_HELP-coverage.fits")
bands = {'b99':1, 'b123':2, 'v':3, 'r':4}
apertures = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]
aperture_sizes = [0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 15.0, 20.0, 30.0, 40.0, 50.0] #arcsec aperture sizes
flux = {}
flux_errors ={}
magnitudes = {}
flux_errors ={}
magnitude_errors = {}
stellarities = {}
for col in t.colnames:
if col.startswith('MAG'):
t[col][np.isclose(t[col], 99.)] = np.nan
stellarity = t['CLASS_STAR_2']
for band in bands:
magnitudes[band] = np.array([t[f"MAG_APER_{ap}_{bands[band]}"] for ap in apertures])
#magnitudes[band] = t['MAG_APER_1_{}'.format(bands[band])]
#for ap in apertures:
# magnitudes[band] = np.stack((magnitudes[band],t['MAG_APER_{}_{}'.format(ap, bands[band])]), axis=0)
#mags_b99 = np.array(t['MAG_APER_r']).T
#mags_r[mags_r == 99] = np.nan
#mags_b123 = np.array(t['MAG_APER_u']).T
#mags_u[mags_u == 99] = np.nan
#mags_v = np.array(t['MAG_APER_g']).T
#mags_g[mags_g == 99] = np.nan
#mags_r = np.array(t['MAG_APER_z']).T
#mags_z[mags_z == 99] = np.nan
#del t
nb_plot_mag_ap_evol(magnitudes['b99'], stellarity)
We will use the 13th (12+1) (aperture number above begin at 0) aperture as target.
nb_plot_mag_vs_apcor(magnitudes['b99'][AP_INDEX], magnitudes['b99'][12], stellarity)
We use magnitudes between 14.3 and 14.6.
nb_plot_mag_ap_evol(magnitudes['b123'], stellarity)
We will use the 13th (12+1) (aperture number above begin at 0) aperture as target.
nb_plot_mag_vs_apcor(magnitudes['b123'][AP_INDEX], magnitudes['b123'][12], stellarity)
We use magnitudes between 14.0 and 16.0.
nb_plot_mag_ap_evol(magnitudes['v'], stellarity)
We will use the 13th (12+1) (aperture number above begin at 0) aperture as target.
nb_plot_mag_vs_apcor(magnitudes['v'][AP_INDEX], magnitudes['v'][12], stellarity)
We use magnitudes between 14.0 and 16.0.
nb_plot_mag_ap_evol(magnitudes['r'], stellarity)
We will use the 13th (12+1) (aperture number above begin at 0) aperture as target.
nb_plot_mag_vs_apcor(magnitudes['r'][AP_INDEX], magnitudes['r'][12], stellarity)
We use magnitudes between 14.0 and 16.0.
#Aperture 6 is 3 arcsec, Aperture 4 is 2 arcsec. We are now using 2 arcsec!
#There appear to be four bands in the catalogue B99/B123/V89/R162 (1/2/3/4)
#Multiple class star - for each band
imported_columns = OrderedDict({
'ID':'voice_id',
'ALPHA_J2000':'voice_ra',
'DELTA_J2000':'voice_dec',
'CLASS_STAR_2':'voice_stellarity',
'MAG_APER_4_1':'m_ap_voice_b99',
'MAGERR_APER_4_1':'merr_ap_voice_b99',
'MAG_AUTO_1':'m_voice_b99',
'MAGERR_AUTO_1':'merr_voice_b99',
'MAG_APER_4_2':'m_ap_voice_b123',
'MAGERR_APER_4_2':'merr_ap_voice_b123',
'MAG_AUTO_2':'m_voice_b123',
'MAGERR_AUTO_2':'merr_voice_b123',
'MAG_APER_4_3':'m_ap_voice_v',
'MAGERR_APER_4_3':'merr_ap_voice_v',
'MAG_AUTO_3':'m_voice_v',
'MAGERR_AUTO_3':'merr_voice_v',
'MAG_APER_4_4':'m_ap_voice_r',
'MAGERR_APER_4_4':'merr_ap_voice_r',
'MAG_AUTO_4':'m_voice_r',
'MAGERR_AUTO_4':'merr_voice_r'
})
catalogue = Table.read("../../dmu0/dmu0_ESIS-VOICE/data/esis_b2vr_cat_03_HELP-coverage.fits")[list(imported_columns)]
for column in imported_columns:
catalogue[column].name = imported_columns[column]
epoch = 2016
# Clean table metadata
catalogue.meta = None
# Index of the target aperture when doing aperture correction
# (see 'sparcs_aperture_correction' notebook).
AP_TARG_INDEX = {
'b99': 12,
'b123': 12,
'v': 12,
'r': 12
}
# Magnitude range for aperture correction.
APCOR_MAG_LIMITS = {
'b99': (14.0, 16.0),
'b123': (14.0, 16.0),
'v': (14.0, 16.0),
'r': (14.0, 16.0)
}
for band in ['b99', 'b123', 'v', 'r']:
# Aperture magnitudes
mag_aper = magnitudes[band][AP_INDEX]
mag_aper_target = magnitudes[band][ AP_TARG_INDEX[band]]
magerr_aper = catalogue["merr_ap_voice_{}".format(band)]
# 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['voice_stellarity'],
mag_min=APCOR_MAG_LIMITS[band][0], mag_max=APCOR_MAG_LIMITS[band][1]
)
print("Aperture correction for VOICE 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["m_ap_voice_{}".format(band)] = mag_aper.data
# 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_voice_{}".format(band)
))
catalogue.add_column(Column(
data = fluxerr_aper * 1.e6,
name = "ferr_ap_voice_{}".format(band)
))
# Auto magnitudes
# Set bad values (99.0) to NaN
mask = (catalogue["m_voice_{}".format(band)] > 90) | (catalogue["merr_voice_{}".format(band)] > 90)
catalogue["m_voice_{}".format(band)][mask] = np.nan
catalogue["merr_voice_{}".format(band)][mask] = np.nan
# Computing the flux columns
flux, fluxerr = mag_to_flux(catalogue["m_voice_{}".format(band)],
catalogue["merr_voice_{}".format(band)])
catalogue.add_column(Column(
data = flux * 1.e6,
name = "f_voice_{}".format(band)
))
catalogue.add_column(Column(
data = fluxerr * 1.e6,
name = "ferr_voice_{}".format(band)
))
# Band-flag column
catalogue.add_column(Column(np.zeros(len(catalogue), dtype=bool),
name="flag_voice_{}".format(band)
))
# 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_ap_voice_b99', 'merr_ap_voice_b123', 'merr_ap_voice_v', 'merr_ap_voice_r']
FLAG_NAME = 'voice_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-S1.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 = "voice_flag_gaia"
catalogue['flag_gaia'].name = GAIA_FLAG_NAME
print("{} sources flagged.".format(np.sum(catalogue[GAIA_FLAG_NAME] > 0)))
catalogue.write("{}/ESIS-VOICE.fits".format(OUT_DIR), overwrite=True)