Herschel Stripe 82 master catalogue¶

Preparation of DECam Legacy Survey data¶

This catalogue comes from dmu0_DECaLS.

In the catalogue, we keep:

  • The object_id as unique object identifier;
  • The position;
  • The u, g, r, i, z, Y aperture magnitude (2”);
  • The u, g, r, i, z, Y kron fluxes and magnitudes.

We check for all ugrizY then only take bands for which there are measurements

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))
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 = "decals_ra"
DEC_COL = "decals_dec"
In [4]:
# Pritine LS catalogue
orig_decals = Table.read("../../dmu0/dmu0_DECaLS/data/DECaLS_Herschel-Stripe-82.fits")
WARNING: UnitsWarning: '1/deg^2' did not parse as fits unit: Numeric factor not supported by FITS [astropy.units.core]
WARNING: UnitsWarning: 'nanomaggy' did not parse as fits unit: At col 0, Unit 'nanomaggy' not supported by the FITS standard.  [astropy.units.core]
WARNING: UnitsWarning: '1/nanomaggy^2' did not parse as fits unit: Numeric factor not supported by FITS [astropy.units.core]
WARNING: UnitsWarning: '1/arcsec^2' did not parse as fits unit: Numeric factor not supported by FITS [astropy.units.core]

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 = ["u", "g", "r", "i", "z", "y"]
band_index = {"u":0,"g":1, "r":2, "i":3, "z":4, "y":5}
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_decals["decam_apflux"][:,band_index[band]])) #np.transpose(np.array( orig_decals["decam_apflux"], dtype=np.float )) 
    flux_errors[band] = np.transpose(np.array(orig_decals["decam_apflux_ivar"][:,band_index[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_decals),0., dtype='float32')
    stellarities[band][np.array( orig_decals["type"]) == "PSF " ] = 1.
    stellarities[band][np.array( orig_decals["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 = {}
/opt/herschelhelp_internal/herschelhelp_internal/utils.py:76: RuntimeWarning: divide by zero encountered in log10
  magnitudes = 2.5 * (23 - np.log10(fluxes)) - 48.6
/opt/herschelhelp_internal/herschelhelp_internal/utils.py:80: RuntimeWarning: invalid value encountered in double_scalars
  errors = 2.5 / np.log(10) * errors_on_fluxes / fluxes
/opt/herschelhelp_internal/herschelhelp_internal/utils.py:76: RuntimeWarning: invalid value encountered in log10
  magnitudes = 2.5 * (23 - np.log10(fluxes)) - 48.6

1.a u band¶

In [6]:
nb_plot_mag_ap_evol(magnitudes['u'], stellarities['u'], labels=apertures)
/opt/anaconda3/envs/herschelhelp_internal/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:703: RuntimeWarning: Mean of empty slice
  warnings.warn("Mean of empty slice", RuntimeWarning)
/opt/anaconda3/envs/herschelhelp_internal/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:703: RuntimeWarning: Mean of empty slice
  warnings.warn("Mean of empty slice", RuntimeWarning)

u band is all nan

In [7]:
nb_plot_mag_vs_apcor(magnitudes['u'][4], 
                     magnitudes['u'][5], 
                     stellarities['u'])
/opt/anaconda3/envs/herschelhelp_internal/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1043: RuntimeWarning: All-NaN slice encountered
  warnings.warn("All-NaN slice encountered", RuntimeWarning)
In [8]:
# Aperture correction
mag_corr['u'] = np.nan
#mag_corr['u'], num, std = aperture_correction(
#    magnitudes['u'][4], magnitudes['u'][5], 
#    stellarities['u'],
#    mag_min=16.0, mag_max=19.0)
#print("Aperture correction for g band:")
#print("Correction: {}".format(mag_corr['g']))
#print("Number of source used: {}".format(num))
#print("RMS: {}".format(std))

I.a - g band¶

In [9]:
nb_plot_mag_ap_evol(magnitudes['g'], stellarities['g'], labels=apertures)

We will use aperture 5 as target.

In [10]:
nb_plot_mag_vs_apcor(magnitudes['g'][4], 
                     magnitudes['g'][5], 
                     stellarities['g'])

We will use magnitudes between 16.0 and 19.0

In [11]:
# Aperture correction

mag_corr['g'], num, std = aperture_correction(
    magnitudes['g'][4], magnitudes['g'][5], 
    stellarities['g'],
    mag_min=16.0, mag_max=19.0)
print("Aperture correction for g band:")
print("Correction: {}".format(mag_corr['g']))
print("Number of source used: {}".format(num)) 
print("RMS: {}".format(std))
Aperture correction for g band:
Correction: -0.0911514235846056
Number of source used: 151015
RMS: 0.02364389650630337

I.b - r band¶

In [12]:
nb_plot_mag_ap_evol(magnitudes['r'], stellarities['r'], labels=apertures)

We will use aperture 5 as target.

In [13]:
nb_plot_mag_vs_apcor(magnitudes['r'][4], 
                     magnitudes['r'][5], 
                     stellarities['r'])

We use magnitudes between 16.0 and 18.0.

In [14]:
# Aperture correction
mag_corr['r'], num, std = aperture_correction(
    magnitudes['r'][4], magnitudes['r'][5], 
    stellarities['r'],
    mag_min=16.0, mag_max=18.0)
print("Aperture correction for r band:")
print("Correction: {}".format(mag_corr['r']))
print("Number of source used: {}".format(num))
print("RMS: {}".format(std))
Aperture correction for r band:
Correction: -0.0465021447682048
Number of source used: 149159
RMS: 0.013977600173198289

I.d - i band¶

In [15]:
nb_plot_mag_ap_evol(magnitudes['i'], stellarities['i'], labels=apertures)
/opt/anaconda3/envs/herschelhelp_internal/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:703: RuntimeWarning: Mean of empty slice
  warnings.warn("Mean of empty slice", RuntimeWarning)
/opt/anaconda3/envs/herschelhelp_internal/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:703: RuntimeWarning: Mean of empty slice
  warnings.warn("Mean of empty slice", RuntimeWarning)
In [16]:
nb_plot_mag_vs_apcor(magnitudes['i'][4], 
                     magnitudes['i'][4], 
                     stellarities['i'])
/opt/anaconda3/envs/herschelhelp_internal/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1043: RuntimeWarning: All-NaN slice encountered
  warnings.warn("All-NaN slice encountered", RuntimeWarning)
In [17]:
# Aperture correction
mag_corr['i'] = np.nan
#mag_corr['i'], num, std = aperture_correction(
#    magnitudes['i'][4], magnitudes['i'][5], 
#    stellarities['i'],
#    mag_min=16.0, mag_max=17.5)
#print("Aperture correction for i band:")
#print("Correction: {}".format(mag_corr['i']))
#print("Number of source used: {}".format(num))
#print("RMS: {}".format(std))

I.e - z band¶

In [18]:
nb_plot_mag_ap_evol(magnitudes['z'], stellarities['z'], labels=apertures)

We will use aperture 4 as target.

In [19]:
nb_plot_mag_vs_apcor(magnitudes['z'][4], 
                     magnitudes['z'][4], 
                     stellarities['z'])

We use magnitudes between 16.0 and 17.5.

In [20]:
# Aperture correction
mag_corr['z'] = np.nan
mag_corr['z'], num, std = aperture_correction(
    magnitudes['z'][4], magnitudes['z'][4], 
    stellarities['z'],
    mag_min=15.0, mag_max=17.0)
print("Aperture correction for z band:")
print("Correction: {}".format(mag_corr['z']))
print("Number of source used: {}".format(num))
print("RMS: {}".format(std))
Aperture correction for z band:
Correction: 0.0
Number of source used: 176285
RMS: 0.0

I.f - Y band¶

In [21]:
nb_plot_mag_ap_evol(magnitudes['y'], stellarities['y'], labels=apertures)
/opt/anaconda3/envs/herschelhelp_internal/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:703: RuntimeWarning: Mean of empty slice
  warnings.warn("Mean of empty slice", RuntimeWarning)
/opt/anaconda3/envs/herschelhelp_internal/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:703: RuntimeWarning: Mean of empty slice
  warnings.warn("Mean of empty slice", RuntimeWarning)
In [22]:
nb_plot_mag_vs_apcor(magnitudes['y'][4], 
                     magnitudes['y'][4], 
                     stellarities['y'])
/opt/anaconda3/envs/herschelhelp_internal/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1043: RuntimeWarning: All-NaN slice encountered
  warnings.warn("All-NaN slice encountered", RuntimeWarning)
In [23]:
# Aperture correction
mag_corr['y'] = np.nan
#mag_corr['y'], num, std = aperture_correction(
#    magnitudes['y'][4], magnitudes['y'][5], 
#    stellarities['y'],
#    mag_min=16.0, mag_max=17.5)
#print("Aperture correction for y band:")
#print("Correction: {}".format(mag_corr['y']))
#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 [24]:
stellarities['g'][np.isclose(stellarities['g'], 1.)] = 0.9
stellarities['g'][np.isclose(stellarities['g'], 0.)] = 0.05

II - Column selection¶

In [25]:
imported_columns = OrderedDict({
        "objid": "decals_id",
        "brickid": "brickid",
        "ra": "decals_ra",
        "dec": "decals_dec",
        "decam_flux": "decam_flux_TEMP",
        "decam_flux_ivar": "decam_flux_ivar_TEMP",
        "decam_apflux": "decam_apflux_TEMP",
        "decam_apflux_ivar": "decam_apflux_ivar_TEMP",
       
    })


catalogue = Table.read("../../dmu0/dmu0_DECaLS/data/DECaLS_Herschel-Stripe-82.fits")[list(imported_columns)]

for column in imported_columns:
    catalogue[column].name = imported_columns[column]


catalogue["decals_id"] = 100000*catalogue["brickid"].astype(np.int64) + catalogue["decals_id"].astype(np.int64)
catalogue.remove_columns("brickid")
    
    
epoch = 2017

#catalogue.add_column(Column(catalogue["decam_flux_TEMP"][:,0], name="f_decam_u"))
catalogue.add_column(Column(catalogue["decam_flux_TEMP"][:,1], name="f_decam_g"))
catalogue.add_column(Column(catalogue["decam_flux_TEMP"][:,2], name="f_decam_r"))
#catalogue.add_column(Column(catalogue["decam_flux_TEMP"][:,3], name="f_decam_i"))
catalogue.add_column(Column(catalogue["decam_flux_TEMP"][:,4], name="f_decam_z"))
#catalogue.add_column(Column(catalogue["decam_flux_TEMP"][:,5], name="f_decam_y"))

#catalogue.add_column(Column(catalogue["decam_flux_ivar_TEMP"][:,0], name="ferr_decam_u"))
catalogue.add_column(Column(catalogue["decam_flux_ivar_TEMP"][:,1], name="ferr_decam_g"))
catalogue.add_column(Column(catalogue["decam_flux_ivar_TEMP"][:,2], name="ferr_decam_r"))
#catalogue.add_column(Column(catalogue["decam_flux_ivar_TEMP"][:,3], name="ferr_decam_i"))
catalogue.add_column(Column(catalogue["decam_flux_ivar_TEMP"][:,4], name="ferr_decam_z"))
#catalogue.add_column(Column(catalogue["decam_flux_ivar_TEMP"][:,5], name="ferr_decam_y"))

#For the aperture fluxes, there are 8 (0-7), we take 4 (2.0")
#DECam aperture fluxes on the co-added images in apertures of radius [0.5,0.75,1.0,1.5,2.0,3.5,5.0,7.0] arcsec in ugrizY
#catalogue.add_column(Column(catalogue["decam_apflux_TEMP"][:,0], name="f_ap_decam_u")[:,4])
catalogue.add_column(Column(catalogue["decam_apflux_TEMP"][:,1], name="f_ap_decam_g")[:,4])
catalogue.add_column(Column(catalogue["decam_apflux_TEMP"][:,2], name="f_ap_decam_r")[:,4])
#catalogue.add_column(Column(catalogue["decam_apflux_TEMP"][:,3], name="f_ap_decam_i")[:,4])
catalogue.add_column(Column(catalogue["decam_apflux_TEMP"][:,4], name="f_ap_decam_z")[:,4])
#catalogue.add_column(Column(catalogue["decam_apflux_TEMP"][:,5], name="f_ap_decam_y")[:,4])

#catalogue.add_column(Column(catalogue["decam_apflux_ivar_TEMP"][:,0], name="ferr_ap_decam_u")[:,4])
catalogue.add_column(Column(catalogue["decam_apflux_ivar_TEMP"][:,1], name="ferr_ap_decam_g")[:,4])
catalogue.add_column(Column(catalogue["decam_apflux_ivar_TEMP"][:,2], name="ferr_ap_decam_r")[:,4])
#catalogue.add_column(Column(catalogue["decam_apflux_ivar_TEMP"][:,3], name="ferr_ap_decam_i")[:,4])
catalogue.add_column(Column(catalogue["decam_apflux_ivar_TEMP"][:,4], name="ferr_ap_decam_z")[:,4])
#catalogue.add_column(Column(catalogue["decam_apflux_ivar_TEMP"][:,5], name="ferr_ap_decam_y")[:,4])

catalogue.remove_columns(["decam_flux_TEMP", 
                          "decam_flux_ivar_TEMP", 
                          "decam_apflux_TEMP", 
                          "decam_apflux_ivar_TEMP"])

# Clean table metadata
catalogue.meta = None
WARNING: UnitsWarning: '1/deg^2' did not parse as fits unit: Numeric factor not supported by FITS [astropy.units.core]
WARNING: UnitsWarning: 'nanomaggy' did not parse as fits unit: At col 0, Unit 'nanomaggy' not supported by the FITS standard.  [astropy.units.core]
WARNING: UnitsWarning: '1/nanomaggy^2' did not parse as fits unit: Numeric factor not supported by FITS [astropy.units.core]
WARNING: UnitsWarning: '1/arcsec^2' did not parse as fits unit: Numeric factor not supported by FITS [astropy.units.core]
In [26]:
flux_to_mag_vect = np.vectorize(flux_to_mag)


# Adding flux and band-flag columns
for col in catalogue.colnames:
    catalogue[col].unit = None
    if col.startswith('f_'):
        
        #Replace 0 flux with NaN and 
        catalogue[col][catalogue[col] == 0.0] = np.nan
        
        #Replace 1/sigma^2 with sigma 
        errcol = "ferr{}".format(col[1:])
        catalogue[errcol][catalogue[errcol] == 0.0] = np.nan
        catalogue[errcol] = np.sqrt(1/np.array(catalogue[errcol]))
        #catalogue[errcol][catalogue[errcol] == None] = np.nan
        
        #Replace nanomaggies with uJy
        #a nanomaggy is approximately 3.631×10-6 Jy - http://www.sdss3.org/dr8/algorithms/magnitudes.php#nmgy
        catalogue[col] = catalogue[col]  * 3.631        
        catalogue[errcol] = catalogue[errcol]  * 3.631
        
        #Compute magnitudes and errors in magnitudes. This function expects Jy so must multiply uJy by 1.e-6
        mag, error = flux_to_mag(np.array(catalogue[col])* 1.e-6, np.array(catalogue[errcol])* 1.e-6)
        
        if 'ap' in col:
            mag += mag_corr[col[-1]]
            catalogue[col],catalogue[errcol] = mag_to_flux(mag,error)
        
        catalogue.add_column(Column(mag, name="m{}".format(col[1:])))
        catalogue.add_column(Column(error, 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:])))
        
#remove units from table
for col in catalogue.colnames:
    catalogue[col].unit = None
    
catalogue.add_column(Column(data=stellarities['g'], name="decals_stellarity")) #Stellarites computed earlier
/opt/herschelhelp_internal/herschelhelp_internal/utils.py:76: RuntimeWarning: invalid value encountered in log10
  magnitudes = 2.5 * (23 - np.log10(fluxes)) - 48.6
In [27]:
catalogue[:10].show_in_notebook()
Out[27]:
<Table masked=True length=10>
idxdecals_iddecals_radecals_decf_decam_gf_decam_rf_decam_zferr_decam_gferr_decam_rferr_decam_zf_ap_decam_gf_ap_decam_rf_ap_decam_zferr_ap_decam_gferr_ap_decam_rferr_ap_decam_zm_decam_gmerr_decam_gflag_decam_gm_decam_rmerr_decam_rflag_decam_rm_decam_zmerr_decam_zflag_decam_zm_ap_decam_gmerr_ap_decam_gm_ap_decam_rmerr_ap_decam_rm_ap_decam_zmerr_ap_decam_zdecals_stellarity
0327488001310.216929045809-0.6229831551054.6585723.159956.57240.1223680.1611960.55554.72777e-062.35467e-055.4468e-051.48916e-072.06961e-077.89212e-0722.22940.0285194False20.48820.00755686False19.51850.0106611False22.21340.034198720.47020.0095429519.55960.01573170.05
1327488001320.216244192498-0.6211937518010.8953321.00071.888940.09679690.1055430.3622668.23116e-071.28561e-062.35024e-061.48917e-072.06962e-077.89211e-0724.020.117382False23.89920.114511False23.20950.208225False24.11130.1964323.62720.17478522.97220.3645910.9
2327488001330.21854126099-0.6225684117690.9398882.318942.884570.1177070.1428030.5172331.0491e-062.2166e-062.92056e-061.48916e-072.06961e-077.89217e-0723.96730.135973False22.98680.066861False22.74980.194684False23.8480.15411723.03580.10137422.73630.2933970.05
3327488001340.215834836373-0.6239668233250.5295272.545045.21480.1040190.1205880.4264987.22324e-072.39483e-065.44071e-061.48916e-072.06961e-077.89215e-0724.59030.21328False22.88580.051444False22.10690.0887981False24.25320.22383822.95180.093829322.06090.1574940.05
4327488001440.180877505286-0.6248990782720.430751.057183.184270.1167690.1407560.51644.25152e-071.19201e-063.09731e-061.48916e-072.0696e-077.89217e-0724.81440.294326False23.83960.144558False22.64250.176076False24.82860.38029623.70930.18850822.67250.2766530.05
5327488001470.132143571314-0.6248721080890.678151.337161.732840.09622730.1056720.3607691.11006e-061.7502e-062.47571e-061.48865e-072.06938e-077.89215e-0724.32170.154062False23.58450.0858026False23.30310.226045False23.78660.14560323.29230.12837422.91570.3461140.9
6327488001480.130282373055-0.6247404820530.1164330.7880290.4424470.0956410.1046820.3595124.51142e-082.51643e-071.00819e-061.48866e-072.06939e-077.89216e-0726.23480.891847False24.15860.144229False24.78530.882219False27.26423.5826625.3980.89285623.89110.8499190.9
7327488001490.131377724037-0.6247228042280.3223750.6866191.692410.09574710.1042850.3613277.21791e-079.90843e-072.59488e-061.48866e-072.06938e-077.89213e-0725.12910.322469False24.30820.164903False23.32870.231803False24.2540.22392823.910.22675722.86470.3302190.9
8327488001520.235602994924-0.6243842490450.6305091.435921.651860.1173480.1418890.5185957.72167e-071.7423e-061.84749e-061.48916e-072.06961e-077.89211e-0724.40080.202073False23.50720.107286False23.35510.340863False24.18070.2093923.29720.12897123.23350.4638050.05
9327488001530.236239923255-0.6238748009920.1067930.8640211.351220.09555330.1052330.3631545.00839e-071.19566e-061.79588e-061.48916e-072.06961e-077.89215e-0726.32860.971462False24.05870.132237False23.57320.291803False24.65080.32282523.7060.18793323.26430.4771360.9

III - Removal of duplicated sources¶

We remove duplicated objects from the input catalogues.

In [28]:
SORT_COLS = [#'merr_ap_decam_u',
             'merr_ap_decam_g',
             'merr_ap_decam_r',
             #'merr_ap_decam_i',
             'merr_ap_decam_z',
             #'merr_ap_decam_y'
]
FLAG_NAME = 'decals_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 24089956 sources.
The cleaned catalogue has 24084781 sources (5175 removed).
The cleaned catalogue has 5169 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 [29]:
gaia = Table.read("../../dmu0/dmu0_GAIA/data/GAIA_Herschel-Stripe-82.fits")
gaia_coords = SkyCoord(gaia['ra'], gaia['dec'])
In [30]:
nb_astcor_diag_plot(catalogue[RA_COL], catalogue[DEC_COL], 
                    gaia_coords.ra, gaia_coords.dec, near_ra0=True)
In [31]:
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.10549986955084023 arcsec
Dec correction: -0.0907548244944989 arcsec
In [32]:
catalogue[RA_COL] = catalogue[RA_COL] +  delta_ra.to(u.deg)
catalogue[DEC_COL] = catalogue[DEC_COL] +  delta_dec.to(u.deg)
catalogue[RA_COL].unit = u.deg
catalogue[DEC_COL].unit = u.deg
In [33]:
nb_astcor_diag_plot(catalogue[RA_COL], catalogue[DEC_COL], 
                    gaia_coords.ra, gaia_coords.dec, near_ra0=True)

IV - Flagging Gaia objects¶

In [34]:
catalogue.add_column(
    gaia_flag_column(SkyCoord(catalogue[RA_COL], catalogue[DEC_COL]), epoch, gaia)
)
In [35]:
GAIA_FLAG_NAME = "decals_flag_gaia"

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

V - Saving to disk¶

In [36]:
catalogue.write("{}/DECaLS.fits".format(OUT_DIR), overwrite=True)