<jupyter_start><jupyter_text>XMM-LSS master catalogue Preparation of the Subaru/XMM-Newton Deep Survey (SXDS) dataThe catalogue is in `dmu0_SXDS`.In the catalogue, we keep:- The position;- The stellarity;- The aperture magnitude B, V, R, i, z (2 arcsec).- The total magnitude B, V, R, i, z (Kron like aperture magnitude).<jupyter_code>from herschelhelp_internal import git_version
print("This notebook was run with herschelhelp_internal version: \n{}".format(git_version()))
import datetime
print("This notebook was executed on: \n{}".format(datetime.datetime.now()))
%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, vstack
import numpy as np
#from numpy.core.defchararray import add

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
OUT_DIR =  os.environ.get('TMP_DIR', "./data_tmp")
try:
    os.makedirs(OUT_DIR)
except FileExistsError:
    pass

RA_COL = "sxds_ra"
DEC_COL = "sxds_dec"<jupyter_output><empty_output><jupyter_text>I - Column selectionFor each band we have 5 independent and overlapping catalogues. We must first stack the catalogues and remove duplicates then merge the bands together.<jupyter_code>#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 # TODO: check this

bands = ["B", "V", "R", "i", "z"]
cat_numbers = [ 1, 2, 3, 4, 5 ]

def imported_columns(band): 
    return OrderedDict({
        'NUMBER': "sxds_{}_id".format(band.lower()),
        'ra': "sxds_{}_ra".format(band.lower()),
        'dec': "sxds_{}_dec".format(band.lower()),
        'CLASS_STAR_{}'.format(band):  "sxds_{}_stellarity".format(band.lower()),
        'twoApertureMag_{}'.format(band): "m_ap_sxds_{}".format(band.lower()),
        'Err_twoApertureMag_{}'.format(band): "merr_ap_sxds_{}".format(band.lower()),
        'MAG_AUTO_{}'.format(band): "m_sxds_{}".format(band.lower()),
        'Err_MAG_AUTO_{}'.format(band): "merr_sxds_{}".format(band.lower()) 
    })


catalogue_B = vstack([ Table.read("../../dmu0/dmu0_SXDS/data/sxdsB{}_dr1.fits".format(1))[
                                   list(imported_columns("B"))],
                      Table.read("../../dmu0/dmu0_SXDS/data/sxdsB{}_dr1.fits".format(2))[
                                   list(imported_columns("B"))],
                      Table.read("../../dmu0/dmu0_SXDS/data/sxdsB{}_dr1.fits".format(3))[
                                   list(imported_columns("B"))],
                      Table.read("../../dmu0/dmu0_SXDS/data/sxdsB{}_dr1.fits".format(4))[
                                   list(imported_columns("B"))],
                      Table.read("../../dmu0/dmu0_SXDS/data/sxdsB{}_dr1.fits".format(5))[
                                   list(imported_columns("B"))]
                     ])
                      

for column in imported_columns("B"):
    catalogue_B[column].name = imported_columns("B")[column]
SORT_COLS = ['merr_ap_sxds_b']
FLAG_NAME = 'sxds_flag_cleaned_b'

nb_orig_sources = len(catalogue_B)

catalogue_B = remove_duplicates(catalogue_B, 'sxds_b_ra', 'sxds_b_dec',  sort_col=SORT_COLS, flag_name=FLAG_NAME)

nb_sources = len(catalogue_B)

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_B[FLAG_NAME])))
catalogue_V = vstack([ Table.read("../../dmu0/dmu0_SXDS/data/sxdsV{}_dr1.fits".format(1))[
                                   list(imported_columns("V"))],
                      Table.read("../../dmu0/dmu0_SXDS/data/sxdsV{}_dr1.fits".format(2))[
                                   list(imported_columns("V"))],
                      Table.read("../../dmu0/dmu0_SXDS/data/sxdsV{}_dr1.fits".format(3))[
                                   list(imported_columns("V"))],
                      Table.read("../../dmu0/dmu0_SXDS/data/sxdsV{}_dr1.fits".format(4))[
                                   list(imported_columns("V"))],
                      Table.read("../../dmu0/dmu0_SXDS/data/sxdsV{}_dr1.fits".format(5))[
                                   list(imported_columns("V"))]
                     ])
                      

for column in imported_columns("V"):
    catalogue_V[column].name = imported_columns("V")[column]
SORT_COLS = ['merr_ap_sxds_v']
FLAG_NAME = 'sxds_flag_cleaned_v'

nb_orig_sources = len(catalogue_V)

catalogue_V = remove_duplicates(catalogue_V, 'sxds_v_ra', 'sxds_v_dec',  sort_col=SORT_COLS, flag_name=FLAG_NAME)

nb_sources = len(catalogue_V)

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_V[FLAG_NAME])))
catalogue_R = vstack([ Table.read("../../dmu0/dmu0_SXDS/data/sxdsR{}_dr1.fits".format(1))[
                                   list(imported_columns("R"))],
                      Table.read("../../dmu0/dmu0_SXDS/data/sxdsR{}_dr1.fits".format(2))[
                                   list(imported_columns("R"))],
                      Table.read("../../dmu0/dmu0_SXDS/data/sxdsR{}_dr1.fits".format(3))[
                                   list(imported_columns("R"))],
                      Table.read("../../dmu0/dmu0_SXDS/data/sxdsR{}_dr1.fits".format(4))[
                                   list(imported_columns("R"))],
                      Table.read("../../dmu0/dmu0_SXDS/data/sxdsR{}_dr1.fits".format(5))[
                                   list(imported_columns("R"))]
                     ])
                      

for column in imported_columns("R"):
    catalogue_R[column].name = imported_columns("R")[column]
SORT_COLS = ['merr_ap_sxds_r']
FLAG_NAME = 'sxds_flag_cleaned_r'

nb_orig_sources = len(catalogue_R)

catalogue_R = remove_duplicates(catalogue_R, 'sxds_r_ra', 'sxds_r_dec',  sort_col=SORT_COLS, flag_name=FLAG_NAME)

nb_sources = len(catalogue_R)

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_R[FLAG_NAME])))
catalogue_i = vstack([ Table.read("../../dmu0/dmu0_SXDS/data/sxdsi{}_dr1.fits".format(1))[
                                   list(imported_columns("i"))],
                      Table.read("../../dmu0/dmu0_SXDS/data/sxdsi{}_dr1.fits".format(2))[
                                   list(imported_columns("i"))],
                      Table.read("../../dmu0/dmu0_SXDS/data/sxdsi{}_dr1.fits".format(3))[
                                   list(imported_columns("i"))],
                      Table.read("../../dmu0/dmu0_SXDS/data/sxdsi{}_dr1.fits".format(4))[
                                   list(imported_columns("i"))],
                      Table.read("../../dmu0/dmu0_SXDS/data/sxdsi{}_dr1.fits".format(5))[
                                   list(imported_columns("i"))]
                     ])
                      

for column in imported_columns("i"):
    catalogue_i[column].name = imported_columns("i")[column]
SORT_COLS = ['merr_ap_sxds_i']
FLAG_NAME = 'sxds_flag_cleaned_i'

nb_orig_sources = len(catalogue_i)

catalogue_i = remove_duplicates(catalogue_i, 'sxds_i_ra', 'sxds_i_dec',  sort_col=SORT_COLS, flag_name=FLAG_NAME)

nb_sources = len(catalogue_i)

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_i[FLAG_NAME])))
catalogue_z = vstack([ Table.read("../../dmu0/dmu0_SXDS/data/sxdsz{}_dr1.fits".format(1))[
                                   list(imported_columns("z"))],
                      Table.read("../../dmu0/dmu0_SXDS/data/sxdsz{}_dr1.fits".format(2))[
                                   list(imported_columns("z"))],
                      Table.read("../../dmu0/dmu0_SXDS/data/sxdsz{}_dr1.fits".format(3))[
                                   list(imported_columns("z"))],
                      Table.read("../../dmu0/dmu0_SXDS/data/sxdsz{}_dr1.fits".format(4))[
                                   list(imported_columns("z"))],
                      Table.read("../../dmu0/dmu0_SXDS/data/sxdsz{}_dr1.fits".format(5))[
                                   list(imported_columns("z"))]
                     ])
                      

for column in imported_columns("z"):
    catalogue_z[column].name = imported_columns("z")[column]
SORT_COLS = ['merr_ap_sxds_z']
FLAG_NAME = 'sxds_flag_cleaned_z'

nb_orig_sources = len(catalogue_z)

catalogue_z = remove_duplicates(catalogue_z, 'sxds_z_ra', 'sxds_z_dec',  sort_col=SORT_COLS, flag_name=FLAG_NAME)

nb_sources = len(catalogue_z)

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_z[FLAG_NAME])))<jupyter_output>The initial catalogue had 842590 sources.
The cleaned catalogue has 818902 sources (23688 removed).
The cleaned catalogue has 23589 sources flagged as having been cleaned<jupyter_text>Merging different bandsSXDS has indivdual extractions from each band. We must therefore merge them as if they were individual catalogues (they have different<jupyter_code>nb_merge_dist_plot(
    SkyCoord(catalogue_B['sxds_b_ra'], catalogue_B['sxds_b_dec']),
    SkyCoord(catalogue_V['sxds_v_ra'], catalogue_V['sxds_v_dec'])
)
catalogue = catalogue_B
catalogue_B['sxds_b_ra'].name = 'ra'
catalogue_B['sxds_b_dec'].name = 'dec'
# Given the graph above, we use 0.8 arc-second radius
catalogue = merge_catalogues(catalogue, catalogue_V, "sxds_v_ra", "sxds_v_dec", radius=0.8*u.arcsec)
del catalogue_B
del catalogue_V
nb_merge_dist_plot(
    SkyCoord(catalogue['ra'], catalogue['dec']),
    SkyCoord(catalogue_R['sxds_r_ra'], catalogue_R['sxds_r_dec'])
)
catalogue = merge_catalogues(catalogue, catalogue_R, "sxds_r_ra", "sxds_r_dec", radius=0.8*u.arcsec)
del catalogue_R
nb_merge_dist_plot(
    SkyCoord(catalogue['ra'], catalogue['dec']),
    SkyCoord(catalogue_i['sxds_i_ra'], catalogue_i['sxds_i_dec'])
)
catalogue = merge_catalogues(catalogue, catalogue_i, "sxds_i_ra", "sxds_i_dec", radius=0.8*u.arcsec)
del catalogue_i
nb_merge_dist_plot(
    SkyCoord(catalogue['ra'], catalogue['dec']),
    SkyCoord(catalogue_z['sxds_z_ra'], catalogue_z['sxds_z_dec'])
)
catalogue = merge_catalogues(catalogue, catalogue_z, "sxds_z_ra", "sxds_z_dec", radius=0.8*u.arcsec)
del catalogue_z
#rename radec colums
catalogue['ra'].name = 'sxds_ra'
catalogue['dec'].name = 'sxds_dec'<jupyter_output><empty_output><jupyter_text>Fill masked values and add fluxes and nans<jupyter_code>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 = ""
        
catalogue =  catalogue.filled()

for col in catalogue.colnames:
    if col.startswith('m_'):
        
        errcol = "merr{}".format(col[1:])
        # catalogue_h[col].name = imported_columns_h[col]
        
        # Turn 99.0 and 99.03 to nans
        catalogue[col][catalogue[col] > 90.] = np.nan
        catalogue[errcol][catalogue[errcol] > 90.] = np.nan 
    
        flux, error = mag_to_flux(np.array(catalogue[col]), np.array(catalogue[errcol]))
        
        # Fluxes are added in µJy
        catalogue.add_column(Column(flux * 1.e6, name="f{}".format(col[1:])))
        catalogue.add_column(Column(error * 1.e6, name="f{}".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:])))
        
# Clean table metadata
catalogue.meta = None
    

catalogue[:10].show_in_notebook()<jupyter_output><empty_output><jupyter_text>Combine stellarities<jupyter_code>stellarity_columns = [column for column in catalogue.colnames
                      if 'stellarity' in column]

print(", ".join(stellarity_columns))
# 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="sxds_stellarity"))

catalogue.remove_columns(stellarity_columns)
catalogue[:10].show_in_notebook()<jupyter_output><empty_output><jupyter_text>II - Removal of duplicated sources We remove duplicated objects from the input catalogues.<jupyter_code>SORT_COLS = ['merr_ap_sxds_b',
            'merr_ap_sxds_v',
            'merr_ap_sxds_r',
            'merr_ap_sxds_i',
            'merr_ap_sxds_z',]
FLAG_NAME = 'sxds_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])))
catalogue['sxds_flag_cleaned'] =  (catalogue['sxds_flag_cleaned'] | 
                                   catalogue['sxds_flag_cleaned_b'] |  
                                   catalogue['sxds_flag_cleaned_v'] | 
                                   catalogue['sxds_flag_cleaned_r'] | 
                                   catalogue['sxds_flag_cleaned_i'] |  
                                   catalogue['sxds_flag_cleaned_z'])
catalogue.remove_columns(['sxds_flag_cleaned_b',
                         'sxds_flag_cleaned_v',
                         'sxds_flag_cleaned_r',
                         'sxds_flag_cleaned_i',
                         'sxds_flag_cleaned_z'])
catalogue['flag_merged'].name = 'sxds_flag_merged'
catalogue[:10].show_in_notebook()<jupyter_output><empty_output><jupyter_text>III - Astrometry correctionWe 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.<jupyter_code>gaia = Table.read("../../dmu0/dmu0_GAIA/data/GAIA_XMM-LSS.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] = catalogue[RA_COL] + delta_ra.to(u.deg)
catalogue[DEC_COL] = catalogue[DEC_COL] + delta_dec.to(u.deg)
nb_astcor_diag_plot(catalogue[RA_COL], catalogue[DEC_COL], 
                    gaia_coords.ra, gaia_coords.dec)<jupyter_output><empty_output><jupyter_text>IV - Flagging Gaia objects<jupyter_code>catalogue.add_column(
    gaia_flag_column(SkyCoord(catalogue[RA_COL], catalogue[DEC_COL]), epoch, gaia)
)
GAIA_FLAG_NAME = "sxds_flag_gaia"

catalogue['flag_gaia'].name = GAIA_FLAG_NAME
print("{} sources flagged.".format(np.sum(catalogue[GAIA_FLAG_NAME] > 0)))
catalogue.add_column(Column(data=np.arange(len(catalogue)), name="sxds_intid"))
for col in catalogue.colnames:
    if col.endswith('_sxds_b'):
        catalogue.rename_column(col, col.replace('_sxds_b', '_sxds-suprime_b'))
    elif col.endswith('_sxds_v'):
        catalogue.rename_column(col, col.replace('_sxds_v', '_sxds-suprime_v'))
    elif col.endswith('_sxds_r'):
        catalogue.rename_column(col, col.replace('_sxds_r', '_sxds-suprime_rc'))
    elif col.endswith('_sxds_i'):
        catalogue.rename_column(col, col.replace('_sxds_i', '_sxds-suprime_ip'))
    elif col.endswith('_sxds_z'):
        catalogue.rename_column(col, col.replace('_sxds_z', '_sxds-suprime_zp'))<jupyter_output><empty_output><jupyter_text>V - Saving to disk<jupyter_code>catalogue.write("{}/SXDS.fits".format(OUT_DIR), overwrite=True)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13822031937051607, 'fraction_numerical': 0.010472112936820181, 'mean_word_length': 3.282705649157582, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 39, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 30, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '131244', 'license': 'MIT', 'n_tokens_mistral': 6209, 'n_tokens_neox': 5868, 'n_words': 962, 'path': 'dmu1/dmu1_ml_XMM-LSS/1.11_SXDS.ipynb', 'repo_name': 'djbsmith/dmu_products', 'size': '17284'}