XMM-LSS master catalogue¶

Preparation of Spitzer Adaptation of the Red-sequence Cluster Survey (SpARCS) data¶

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:

  • The internal identifier (this one is only in HeDaM data);
  • The position;
  • The ugrzy magnitudes in the 8th aperture (11×0.186=2.046 arcsec).
  • The "auto" magnitudes.

The maps on the web page indicate they were observed in 2012 (or late 2011). Let's use 2012 as epoch.

In [1]:
from herschelhelp_internal import git_version
print("This notebook was run with herschelhelp_internal version: \n{}".format(git_version()))
This notebook was run with herschelhelp_internal version: 
33f5ec7 (Wed Dec 6 16:56:17 2017 +0000)
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
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
In [3]:
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"

I - Parametres for 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 [4]:
# We are using the aperture index 7 (0 base) that correspond to 11 pix * 0.186 arcsec/pix = 2.046 arcsec
AP_INDEX = 7
In [5]:
t = Table.read("../../dmu0/dmu0_SpARCS/data/SpARCS_HELP_XMM-LSS.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

mags_y = np.array(t['MAG_APER_y']).T
mags_y[mags_y == 99] = np.nan

del t
WARNING: UnitsWarning: '""' did not parse as fits unit: Invalid character at col 0 [astropy.units.core]

I.a r-band¶

In [6]:
nb_plot_mag_ap_evol(mags_r, stellarity)

We will use the 16th (aperture number above begin to 0) aperture as target.

In [7]:
nb_plot_mag_vs_apcor(mags_r[AP_INDEX], mags_r[15], stellarity)
/opt/herschelhelp_internal/herschelhelp_internal/utils.py:131: RuntimeWarning: invalid value encountered in greater_equal
  mask &= (mag >= mag_min)
/opt/herschelhelp_internal/herschelhelp_internal/utils.py:133: RuntimeWarning: invalid value encountered in less_equal
  mask &= (mag <= mag_max)

We use magnitudes between 17 and 17.9.

I.b u-band¶

In [8]:
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?

In [9]:
nb_plot_mag_vs_apcor(mags_u[AP_INDEX], mags_u[15], stellarity)
/opt/herschelhelp_internal/herschelhelp_internal/utils.py:131: RuntimeWarning: invalid value encountered in greater_equal
  mask &= (mag >= mag_min)
/opt/herschelhelp_internal/herschelhelp_internal/utils.py:133: RuntimeWarning: invalid value encountered in less_equal
  mask &= (mag <= mag_max)

We use magnitudes between 17 and 17.9.

I.c g-band¶

In [10]:
nb_plot_mag_ap_evol(mags_g, stellarity)

We will use the 16th (aperture number above begin to 0) aperture as target.

In [11]:
nb_plot_mag_vs_apcor(mags_g[AP_INDEX], mags_g[15], stellarity)
/opt/herschelhelp_internal/herschelhelp_internal/utils.py:131: RuntimeWarning: invalid value encountered in greater_equal
  mask &= (mag >= mag_min)
/opt/herschelhelp_internal/herschelhelp_internal/utils.py:133: RuntimeWarning: invalid value encountered in less_equal
  mask &= (mag <= mag_max)

We use magnitudes between 17.2 and 18.

I.d z-band¶

In [12]:
nb_plot_mag_ap_evol(mags_z, stellarity)

We will use the 16th (aperture number above begin to 0) aperture as target.

In [13]:
nb_plot_mag_vs_apcor(mags_z[AP_INDEX], mags_z[15], stellarity)
/opt/herschelhelp_internal/herschelhelp_internal/utils.py:131: RuntimeWarning: invalid value encountered in greater_equal
  mask &= (mag >= mag_min)
/opt/herschelhelp_internal/herschelhelp_internal/utils.py:133: RuntimeWarning: invalid value encountered in less_equal
  mask &= (mag <= mag_max)

We use magnitudes between 16 and 17.

I.e y-band¶

In [14]:
nb_plot_mag_ap_evol(mags_y, stellarity)

We will use the 16th (aperture number above begin to 0) aperture as target.

In [15]:
nb_plot_mag_vs_apcor(mags_y[AP_INDEX], mags_y[15], stellarity)
/opt/herschelhelp_internal/herschelhelp_internal/utils.py:131: RuntimeWarning: invalid value encountered in greater_equal
  mask &= (mag >= mag_min)
/opt/herschelhelp_internal/herschelhelp_internal/utils.py:133: RuntimeWarning: invalid value encountered in less_equal
  mask &= (mag <= mag_max)

We use mags between 16.5 and 17.5

II - Column selection¶

In [16]:
# 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,
    'y': 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.),
    'y': (16.5, 17.5)
}

epoch = 2012

sparcs_tmp = Table.read("../../dmu0/dmu0_SpARCS/data/SpARCS_HELP_XMM-LSS.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',
    'MAG_APER_y', 'MAGERR_APER_y', 'MAG_AUTO_y', 'MAGERR_AUTO_y']


# 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', 'y']:
    
    # 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_megacam_{}".format(band)
    ))
    catalogue.add_column(Column(
        data = magerr_aper.data,
        name = "merr_ap_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_megacam_{}".format(band)
    ))
    catalogue.add_column(Column(
        data = fluxerr_aper * 1.e6,
        name = "ferr_ap_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_megacam_{}".format(band)
    ))
    catalogue.add_column(Column(
        data = magerr_auto,
        name = "merr_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_megacam_{}".format(band)
    ))
    catalogue.add_column(Column(
        data = fluxerr * 1.e6,
        name = "ferr_megacam_{}".format(band)
    ))
    
    # Band-flag column
    catalogue.add_column(Column(np.zeros(len(catalogue), dtype=bool), 
                             name="flag_megacam_{}".format(band)
    ))
        
# TODO: Set to True the flag columns for fluxes that should not be used for SED fitting.

del sparcs_tmp
WARNING: UnitsWarning: '""' did not parse as fits unit: Invalid character at col 0 [astropy.units.core]
/opt/anaconda3/envs/herschelhelp_internal/lib/python3.6/site-packages/astropy/table/column.py:1096: MaskedArrayFutureWarning: setting an item on a masked array which has a shared mask will not copy the mask and also change the original mask array in the future.
Check the NumPy 1.11 release notes for more information.
  ma.MaskedArray.__setitem__(self, index, value)
/opt/herschelhelp_internal/herschelhelp_internal/utils.py:131: RuntimeWarning: invalid value encountered in greater_equal
  mask &= (mag >= mag_min)
/opt/herschelhelp_internal/herschelhelp_internal/utils.py:133: RuntimeWarning: invalid value encountered in less_equal
  mask &= (mag <= mag_max)
Aperture correction for SpARCS band u:
Correction: -0.21997451782226562
Number of source used: 325
RMS: 0.06430433467975358

Aperture correction for SpARCS band g:
Correction: -0.23264789581298828
Number of source used: 487
RMS: 0.06903040173673632

Aperture correction for SpARCS band r:
Correction: -0.2197418212890625
Number of source used: 821
RMS: 0.10640508916122655

Aperture correction for SpARCS band z:
Correction: -0.2187061309814453
Number of source used: 1058
RMS: 0.06586610521457599

Aperture correction for SpARCS band y:
Correction: -0.2814064025878906
Number of source used: 521
RMS: 0.055379872136121035

In [17]:
catalogue[:10].show_in_notebook()
Out[17]:
<Table masked=True length=10>
idxsparcs_intidsparcs_rasparcs_decsparcs_stellaritym_ap_megacam_umerr_ap_megacam_uf_ap_megacam_uferr_ap_megacam_um_megacam_umerr_megacam_uf_megacam_uferr_megacam_uflag_megacam_um_ap_megacam_gmerr_ap_megacam_gf_ap_megacam_gferr_ap_megacam_gm_megacam_gmerr_megacam_gf_megacam_gferr_megacam_gflag_megacam_gm_ap_megacam_rmerr_ap_megacam_rf_ap_megacam_rferr_ap_megacam_rm_megacam_rmerr_megacam_rf_megacam_rferr_megacam_rflag_megacam_rm_ap_megacam_zmerr_ap_megacam_zf_ap_megacam_zferr_ap_megacam_zm_megacam_zmerr_megacam_zf_megacam_zferr_megacam_zflag_megacam_zm_ap_megacam_ymerr_ap_megacam_yf_ap_megacam_yferr_ap_megacam_ym_megacam_ymerr_megacam_yf_megacam_yferr_megacam_yflag_megacam_y
degdegmagmagmagmagmagmagmagmagmagmag
0133.5163489916-3.774445153420.982108nannannannannannannannanFalsenannannannannannannannanFalse21.45940.008622499.467510.075187321.69980.008516277.587080.0595113FalsenannannannannannannannanFalsenannannannannannannannanFalse
1233.5164313574-3.607826040190.754893nannannannannannannannanFalsenannannannannannannannanFalse23.94290.04086120.9612330.036175624.02110.03693170.8944540.0304252FalsenannannannannannannannanFalsenannannannannannannannanFalse
2333.516451123-3.536405194060.948239nannannannannannannannanFalsenannannannannannannannanFalse23.4750.02843891.479170.038744223.55930.02900411.368640.0365617FalsenannannannannannannannanFalsenannannannannannannannanFalse
3433.5164923954-3.802985252240.796541nannannannannannannannanFalsenannannannannannannannanFalse24.05310.05078090.8684840.040619824.35960.04349050.654890.0262324FalsenannannannannannannannanFalsenannannannannannannannanFalse
4533.5165013484-3.670289099340.954353nannannannannannannannanFalsenannannannannannannannanFalse23.48080.03237851.471210.043873923.69790.03046431.204540.0337978FalsenannannannannannannannanFalsenannannannannannannannanFalse
5633.5165249482-3.39290102880.659968nannannannannannannannanFalsenannannannannannannannanFalse24.67440.07232910.4900710.032647324.90830.05489680.3950720.0199756FalsenannannannannannannannanFalsenannannannannannannannanFalse
6733.5165384578-3.404277470410.837242nannannannannannannannanFalsenannannannannannannannanFalse24.09220.04603020.8377630.035517324.26950.03927040.7115730.0257371FalsenannannannannannannannanFalsenannannannannannannannanFalse
7833.5165723264-3.263240542380.707974nannannannannannannannanFalsenannannannannannannannanFalse23.91390.0393870.9872580.035814524.11060.0384150.8237160.0291443FalsenannannannannannannannanFalsenannannannannannannannanFalse
8933.5166030275-3.634263113770.969087nannannannannannannannanFalsenannannannannannannannanFalse22.76160.0197672.853350.051948422.98310.01957162.326880.0419445FalsenannannannannannannannanFalsenannannannannannannannanFalse
91033.5166238733-3.359313572810.777242nannannannannannannannanFalsenannannannannannannannanFalse23.57110.03255491.353860.040594423.70860.03864081.192740.0424491FalsenannannannannannannannanFalsenannannannannannannannanFalse

II - Removal of duplicated sources¶

We remove duplicated objects from the input catalogues.

In [18]:
SORT_COLS = ['merr_ap_megacam_r', 'merr_ap_megacam_u',
             'merr_ap_megacam_g', 'merr_ap_megacam_z',
             'merr_ap_megacam_y']
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])))
/opt/anaconda3/envs/herschelhelp_internal/lib/python3.6/site-packages/astropy/table/column.py:1096: MaskedArrayFutureWarning: setting an item on a masked array which has a shared mask will not copy the mask and also change the original mask array in the future.
Check the NumPy 1.11 release notes for more information.
  ma.MaskedArray.__setitem__(self, index, value)
The initial catalogue had 447768 sources.
The cleaned catalogue has 447768 sources (0 removed).
The cleaned catalogue has 0 sources flagged as having been cleaned

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 [19]:
gaia = Table.read("../../dmu0/dmu0_GAIA/data/GAIA_XMM-LSS.fits")
gaia_coords = SkyCoord(gaia['ra'], gaia['dec'])
In [20]:
nb_astcor_diag_plot(catalogue[RA_COL], catalogue[DEC_COL], 
                    gaia_coords.ra, gaia_coords.dec)
In [21]:
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))
RA correction: 0.1243786283438908 arcsec
Dec correction: -0.10356365148096458 arcsec
In [22]:
catalogue[RA_COL] +=  delta_ra.to(u.deg)
catalogue[DEC_COL] += delta_dec.to(u.deg)
In [23]:
nb_astcor_diag_plot(catalogue[RA_COL], catalogue[DEC_COL], 
                    gaia_coords.ra, gaia_coords.dec)

IV - Flagging Gaia objects¶

In [24]:
catalogue.add_column(
    gaia_flag_column(SkyCoord(catalogue[RA_COL], catalogue[DEC_COL]), epoch, gaia)
)
In [25]:
GAIA_FLAG_NAME = "sparcs_flag_gaia"

catalogue['flag_gaia'].name = GAIA_FLAG_NAME
print("{} sources flagged.".format(np.sum(catalogue[GAIA_FLAG_NAME] > 0)))
6126 sources flagged.

V - Saving to disk¶

In [26]:
#We rename megacam to sparcs so that when we merge all the
#megacam fluxes from different surveys we have access to the right column name
for col in catalogue.colnames:
    if 'megacam' in col:
        catalogue[col].name = col.replace('megacam', 'sparcs')
In [27]:
catalogue.write("{}/SpARCS.fits".format(OUT_DIR), overwrite=True)