<jupyter_start><jupyter_text>Herschel Stripe 82 master catalogue Preparation of VISTA/CFHT Telescope / VIRCAM/WIRCAM Camera dataVISTA/CFHT Telescope / VIRCAM/WIRCAM Camera VICS82 catalogue: the catalogue comes from `dmu0_VICS82`.The catalogue is described here: https://arxiv.org/pdf/1705.05451.pdfIn the catalogue, we keep:- The identifier (it's unique in the catalogue);- The position;- The stellarity;- The magnitude for each band.- The auto magnitude to be used as total 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
import numpy as np

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, ebv

from herschelhelp.database import get_filters
OUT_DIR =  os.environ.get('TMP_DIR', "./data_tmp")
try:
    os.makedirs(OUT_DIR)
except FileExistsError:
    pass

RA_COL = "vics82_ra"
DEC_COL = "vics82_dec"<jupyter_output><empty_output><jupyter_text>I - Column selection<jupyter_code>imported_columns = OrderedDict({
        'VICS82ID': "vics82_id",
        'ALPHA_J2000': "vics82_ra",
        'DELTA_J2000': "vics82_dec",
        'CLASS_STAR':  "vics82_stellarity",
        'MAG_AUTO': "m_vics82_ks", 
        'MAG_APER': "m_ap_vics82_ks", 
        'JMAG_AUTO': "m_vics82_j", 
        'JMAG_APER': "m_ap_vics82_j", 
        'IMAGE': 'image_name',
})

#No error column?

catalogue = Table.read("../../dmu0/dmu0_VICS82/data/VICS82_FULL_SDSS_FEB2017_K22_HELP-coverage_intIDs.fits")[list(imported_columns)]
for column in imported_columns:
    catalogue[column].name = imported_columns[column]

epoch = 2011

# Clean table metadata
catalogue.meta = None

#replace list of aperture magnitudes with appropriate choice from list
# Apertures (1", 1.5", 2", 2.5" and 3"
catalogue["m_ap_vics82_j"] = catalogue["m_ap_vics82_j"][:, 2]
catalogue["m_ap_vics82_ks"] = catalogue["m_ap_vics82_ks"][:, 2]
#catalogue["m_ap_vics82_j"]  = [el[0] for el in catalogue["m_ap_vics82_j"] ]
#catalogue["m_ap_vics82_j"]  = [el[0] for el in catalogue["m_ap_vics82_j"] ]
# Adding flux and band-flag columns
for col in catalogue.colnames:
    if col.startswith('m_'):
        
        # Add err col with all nan values because errors absent
        errcol = "merr{}".format(col[1:])
        nan_values = np.empty(len(catalogue))
        nan_values[:] = np.NAN
        catalogue.add_column(Column(nan_values, name=errcol))
        
        # Some object have a magnitude to 0 or 99., we suppose this means missing value
        mask = ((catalogue[col] <= 0) |
                (catalogue[errcol] <= 0) |
                (catalogue[col] > 50.))
        catalogue[col][mask] = np.nan
        catalogue[errcol][mask] = 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:])))
        
# TODO: Set to True the flag columns for fluxes that should not be used for SED fitting.
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 = ['m_ap_vics82_j', 'm_ap_vics82_ks']
FLAG_NAME = 'vics82_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])))<jupyter_output>The initial catalogue had 15318703 sources.
The cleaned catalogue has 13914761 sources (1403942 removed).
The cleaned catalogue has 1239734 sources flagged as having been cleaned<jupyter_text>Determination of camera usedUsing the service at https://apps.automeris.io/wpd/ I extracted positions from Figure 2 in the VICS82 paper.![VICS82](../../dmu0/dmu0_VICS82/regions.png)The positions from right to left are below in decimal hours in x and decimal degrees in y:positions= [    [-0.3135593220338988, -0.0710346572223628],    [0.8983050847457625, 0.06283835061978316],    [0.9915254237288131, 0.004249936756892225],    [1.3855932203389827, 0.008955223880596108],    [1.4533898305084745, 0.2485707057930675],    [1.4872881355932202, 0.24897546167467688],    [2.3644067796610164, 0.02064254996205328],]which can be roughly cut usingIs the objects in the grey region between ih and 23h?in_region1 = ((catalogue[RA_COL] > (360. + positions[0][0]*15.) )              | (catalogue[RA_COL] < positions[1][0]*15. )              | ( (catalogue[RA_COL] < positions[2][0]*15. ) & (catalogue[DEC_COL] < positions[2][1]*15. ))             )          Small aquare at 1.4 hrsin_region2 = (((catalogue[RA_COL] > positions[3][0]*15.) & (catalogue[RA_COL] < positions[4][0]*15.))              | ((catalogue[DEC_COL] > positions[4][1]*15.) & (catalogue[RA_COL] > positions[4][0]*15.) & (catalogue[RA_COL] < positions[5][0]*15.)))REgion at far left beyond 2.4 hrin_region3 = ((catalogue[RA_COL] > positions[6][0]*15.) & (catalogue[RA_COL] < 50.))cfht = in_region1 | in_region2 | in_region3However, <NAME> informed us that the camera can be found by taking teh image name which contains a flag for the camera used.<jupyter_code>cfht = np.char.array(catalogue['image_name'].astype(str)).startswith('C')



column_bases = ['f', 'ferr', 'm', 'merr', 'flag', 'f_ap', 'ferr_ap', 'm_ap', 'merr_ap']
bands = ['j', 'ks']

for band in bands:
    for col_base in column_bases:
        catalogue.add_column(Column(data=catalogue[col_base + '_vics82_' + band], 
                                    name = col_base + '_wircam_' + band))   
        if col_base == 'flag':
            catalogue[col_base + '_wircam_' + band][~cfht] = False
        else:
            catalogue[col_base + '_wircam_' + band][~cfht] = np.nan
        catalogue.add_column(Column(data=catalogue[col_base + '_vics82_' + band], 
                                    name = col_base + '_vista_' + band))
        if col_base == 'flag':
            catalogue[col_base + '_vista_' + band][cfht] = False
        else:
            catalogue[col_base + '_vista_' + band][cfht] = np.nan
        catalogue.remove_column(col_base + '_vics82_' + band)

print("There were {} CFHT fluxes and {} VISTA fluxes out of {} total.".format(np.sum(cfht), 
                                                                              np.sum(~cfht), 
                                                                              len(catalogue)))


catalogue.remove_column('image_name')
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_Herschel-Stripe-82.fits")
gaia_coords = SkyCoord(gaia['ra'], gaia['dec'])
nb_astcor_diag_plot(catalogue[RA_COL], catalogue[DEC_COL], 
                    gaia_coords.ra, gaia_coords.dec, near_ra0=True)
delta_ra, delta_dec =  astrometric_correction(
    SkyCoord(catalogue[RA_COL], catalogue[DEC_COL]),
    gaia_coords, near_ra0=True
)

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, near_ra0=True)<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 = "vics82_flag_gaia"

catalogue['flag_gaia'].name = GAIA_FLAG_NAME
print("{} sources flagged.".format(np.sum(catalogue[GAIA_FLAG_NAME] > 0)))<jupyter_output>194610 sources flagged.<jupyter_text>V - Flagging objects near bright stars  VI - Saving to disk<jupyter_code>catalogue.write("{}/VICS82.fits".format(OUT_DIR), overwrite=True)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1252320122284092, 'fraction_numerical': 0.05306256141500164, 'mean_word_length': 3.9891067538126364, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 35, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '169296', 'license': 'MIT', 'n_tokens_mistral': 3476, 'n_tokens_neox': 3137, 'n_words': 686, 'path': 'dmu1/dmu1_ml_Herschel-Stripe-82/1.3_VICS82.ipynb', 'repo_name': 'djbsmith/dmu_products', 'size': '9162'}