SA13 master catalogue¶
Preparation of Legacy Survey data¶
The catalogue comes from dmu0_LegacySurvey.
In the catalogue, we keep:
- The identifier (it's unique in the catalogue);
- The position;
- The stellarity;
- The aperture fluxes. Are these aperture corrected?
- The kron magnitude to be used as total magnitude (no “auto” magnitude is 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))
plt.style.use('ggplot')
from collections import OrderedDict
import os
from astropy import units as u
from astropy import visualization as vis
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, nb_plot_mag_ap_evol, \
nb_plot_mag_vs_apcor, remove_duplicates
from herschelhelp_internal.utils import astrometric_correction, mag_to_flux, aperture_correction, flux_to_mag
In [3]:
OUT_DIR = os.environ.get('TMP_DIR', "./data_tmp")
try:
os.makedirs(OUT_DIR)
except FileExistsError:
pass
RA_COL = "legacy_ra"
DEC_COL = "legacy_dec"
In [4]:
# Pritine LS catalogue
orig_legacy = Table.read("../../dmu0/dmu0_LegacySurvey/data/LegacySurvey-dr4_SA13.fits")
I - Aperture correction¶
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:
- The evolution of the magnitudes of the objects by plotting on the same plot aperture number vs the mean magnitude.
- The mean gain (loss when negative) of magnitude is each aperture compared to the previous (except for the first of course).
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.
In [5]:
bands = ["g", "r", "z"]
apertures = [0, 1, 2, 3, 4, 5, 6, 7]
aperture_sizes = [0.5, 0.75, 1.0, 1.5, 2.0, 3.5, 5.0, 7.0] #arcsec aperture sizes
flux = {}
flux_errors ={}
magnitudes = {}
flux_errors ={}
magnitude_errors = {}
stellarities = {}
flux_to_mag_vect = np.vectorize(flux_to_mag)
for band in bands:
flux[band] = np.transpose(np.array( orig_legacy["apflux_{}".format(band)], dtype=np.float ))
flux_errors[band] = np.transpose(np.array( orig_legacy["apflux_ivar_{}".format(band)], dtype=np.float ))
magnitudes[band], magnitude_errors[band] = flux_to_mag_vect(flux[band] * 3.631e-6 ,flux_errors[band] * 3.631e-6)
stellarities[band] = np.full(len(orig_legacy),0., dtype='float32')
stellarities[band][np.array( orig_legacy["type"]) == "PSF" ] = 1.
# Some sources have an infinite magnitude
mask = np.isinf(magnitudes[band])
magnitudes[band][mask] = np.nan
magnitude_errors[band][mask] = np.nan
mag_corr = {}
I.a - g band¶
In [6]:
nb_plot_mag_ap_evol(magnitudes['g'], stellarities['g'], labels=apertures)
We will use aperture 5 as target.
In [7]:
nb_plot_mag_vs_apcor(magnitudes['g'][4],
magnitudes['g'][5],
stellarities['g'])
We will use magnitudes between 17.0 and 18.5
In [8]:
# Aperture correction
mag_corr['g'], num, std = aperture_correction(
magnitudes['g'][4], magnitudes['g'][5],
stellarities['g'],
mag_min=17.0, mag_max=18.5)
print("Aperture correction for g band:")
print("Correction: {}".format(mag_corr['g']))
print("Number of source used: {}".format(num))
print("RMS: {}".format(std))
I.b - r band¶
In [9]:
nb_plot_mag_ap_evol(magnitudes['r'], stellarities['r'], labels=apertures)
We will use aperture 5 as target.
In [10]:
nb_plot_mag_vs_apcor(magnitudes['r'][4],
magnitudes['r'][5],
stellarities['r'])
We use magnitudes between 16.5 and 18.
In [11]:
# Aperture correction
mag_corr['r'], num, std = aperture_correction(
magnitudes['r'][4], magnitudes['r'][5],
stellarities['r'],
mag_min=16.5, mag_max=18.5)
print("Aperture correction for r band:")
print("Correction: {}".format(mag_corr['r']))
print("Number of source used: {}".format(num))
print("RMS: {}".format(std))
I.c - z band¶
In [12]:
nb_plot_mag_ap_evol(magnitudes['z'], stellarities['z'], labels=apertures)
We will use aperture 5 as target.
In [13]:
nb_plot_mag_vs_apcor(magnitudes['z'][4],
magnitudes['z'][4],
stellarities['z'])
We use magnitudes between 16.0 and 17.5.
In [14]:
# Aperture correction
mag_corr['z'], num, std = aperture_correction(
magnitudes['z'][4], magnitudes['z'][5],
stellarities['z'],
mag_min=16.0, mag_max=17.5)
print("Aperture correction for z band:")
print("Correction: {}".format(mag_corr['z']))
print("Number of source used: {}".format(num))
print("RMS: {}".format(std))
II - Stellarity¶
Legacy Survey does not provide a 0 to 1 stellarity so we replace items flagged as PSF accpording to the following table:
\begin{equation*} P(star) = \frac{ \prod_{i} P(star)_i }{ \prod_{i} P(star)_i + \prod_{i} P(galaxy)_i } \end{equation*}where $i$ is the band, and with using the same probabilities as UKDISS:
| HSC flag | UKIDSS flag | Meaning | P(star) | P(galaxy) | P(noise) | P(saturated) |
|---|---|---|---|---|---|---|
| -9 | Saturated | 0.0 | 0.0 | 5.0 | 95.0 | |
| -3 | Probable galaxy | 25.0 | 70.0 | 5.0 | 0.0 | |
| -2 | Probable star | 70.0 | 25.0 | 5.0 | 0.0 | |
| 0 | -1 | Star | 90.0 | 5.0 | 5.0 | 0.0 |
| 0 | Noise | 5.0 | 5.0 | 90.0 | 0.0 | |
| 1 | +1 | Galaxy | 5.0 | 90.0 | 5.0 | 0.0 |
In [15]:
stellarities['g'][np.isclose(stellarities['g'], 1.)] = 0.9
stellarities['g'][np.isclose(stellarities['g'], 0.)] = 0.05
In [16]:
orig_legacy.add_column(Column(data=stellarities['g'], name="stellarity")) #Stelarites computed earlier
II - Column selection¶
In [17]:
imported_columns = OrderedDict({
"objid": "legacy_id",
"ra": "legacy_ra",
"dec": "legacy_dec",
"flux_g": "f_bass_g",
"flux_ivar_g": "ferr_bass_g",
"apflux_g": "f_ap_bass_g",
"apflux_ivar_g": "ferr_ap_bass_g",
"flux_r": "f_bass_r",
"flux_ivar_r": "ferr_bass_r",
"apflux_r": "f_ap_bass_r",
"apflux_ivar_r": "ferr_ap_bass_r",
"flux_z": "f_bass_z",
"flux_ivar_z": "ferr_bass_z",
"apflux_z": "f_ap_bass_z",
"apflux_ivar_z": "ferr_ap_bass_z",
"stellarity": "legacy_stellarity"
})
catalogue = orig_legacy[list(imported_columns)]
for column in imported_columns:
catalogue[column].name = imported_columns[column]
epoch = 2017
# Clean table metadata
catalogue.meta = None
In [18]:
# Adding flux and band-flag columns
for col in catalogue.colnames:
if col.startswith('f_'):
errcol = "ferr{}".format(col[1:])
#First we take aperture 4 if it is an aperture flux
if 'ap' in col:
catalogue[col] = catalogue[col][:, 4]
catalogue[errcol] = catalogue[errcol][:, 4]
#Convert nanomaggies to uJy
# 1 nanomaggy = 1.e-9 maggy
# 1 maggy = 3631 Jy
# 1 nanomaggy = 3.631×10-6 Jy
catalogue[col] = catalogue[col] * 3.631 #* 1.e9
catalogue[errcol] = (1/np.sqrt(catalogue[errcol])) * 3.631 #* 1.e9
catalogue[col].unit = u.microjansky
catalogue[errcol].unit = u.microjansky
mag, magerror = flux_to_mag(np.array(catalogue[col])* 1.e-6, np.array(catalogue[errcol])* 1.e-6)
mag[mag == np.inf] = np.nan #The very low fluxes yield infinite mags
magerror[magerror == np.inf] = np.nan
if 'ap' in col:
mag += mag_corr[col[-1]]
catalogue[col],catalogue[errcol] = mag_to_flux(mag,magerror)
# Add magnitudes
catalogue.add_column(Column(mag , name="m{}".format(col[1:])))
catalogue.add_column(Column(magerror , name="m{}".format(errcol[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.
In [19]:
catalogue[:10].show_in_notebook()
Out[19]:
III - Removal of duplicated sources¶
We remove duplicated objects from the input catalogues.
In [20]:
SORT_COLS = [
'merr_ap_bass_g', 'merr_ap_bass_r', 'merr_ap_bass_z']
FLAG_NAME = 'legacy_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 [21]:
gaia = Table.read("../../dmu0/dmu0_GAIA/data/GAIA_SA13.fits")
gaia_coords = SkyCoord(gaia['ra'], gaia['dec'])
In [22]:
nb_astcor_diag_plot(catalogue[RA_COL], catalogue[DEC_COL],
gaia_coords.ra, gaia_coords.dec)
In [23]:
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 [24]:
catalogue[RA_COL].unit = u.deg
catalogue[DEC_COL].unit = u.deg
catalogue[RA_COL] = catalogue[RA_COL] + delta_ra.to(u.deg)
catalogue[DEC_COL] = catalogue[DEC_COL] + delta_dec.to(u.deg)
In [25]:
nb_astcor_diag_plot(catalogue[RA_COL], catalogue[DEC_COL],
gaia_coords.ra, gaia_coords.dec)
IV - Flagging Gaia objects¶
In [26]:
catalogue.add_column(
gaia_flag_column(SkyCoord(catalogue[RA_COL], catalogue[DEC_COL]), epoch, gaia)
)
In [27]:
GAIA_FLAG_NAME = "legacy_flag_gaia"
catalogue['flag_gaia'].name = GAIA_FLAG_NAME
print("{} sources flagged.".format(np.sum(catalogue[GAIA_FLAG_NAME] > 0)))
V - Saving to disk¶
In [28]:
catalogue.write("{}/LegacySurvey.fits".format(OUT_DIR), overwrite=True)