In [1]:
import numpy as np
from astropy.table import Table, join
from astropy import units as u
from astropy.coordinates import SkyCoord, search_around_sky
from IPython.display import clear_output
import pickle
import os
In [2]:
from mltier1 import (get_center, Field, MultiMLEstimatorU, parallel_process, get_sigma_all)
In [3]:
%load_ext autoreload
In [4]:
%autoreload
In [5]:
from IPython.display import clear_output
In [6]:
%pylab inline
In [7]:
save_intermediate = True
plot_intermediate = True
In [8]:
idp = "idata/gaus_pdf_v0.8"
In [9]:
if not os.path.isdir(idp):
os.makedirs(idp)
In [10]:
# Busy week Edinburgh 2017
ra_down = 172.09
ra_up = 187.5833
dec_down = 46.106
dec_up = 56.1611
In [11]:
# Busy week Hatfield 2017
ra_down = 170.
ra_up = 190.
dec_down = 46.8
dec_up = 55.9
In [12]:
# Full field July 2017
ra_down = 160.
ra_up = 232.
dec_down = 42.
dec_up = 62.
In [13]:
field = Field(170.0, 190.0, 46.8, 55.9)
In [14]:
field_full = Field(160.0, 232.0, 42.0, 62.0)
In [15]:
combined_all = Table.read("pw.fits")
In [16]:
#lofar_all = Table.read("data/LOFAR_HBA_T1_DR1_catalog_v0.9.gaus.fits")
lofar_all = Table.read("data/LOFAR_HBA_T1_DR1_catalog_v0.99.gaus.fits")
In [17]:
np.array(combined_all.colnames)
Out[17]:
In [18]:
np.array(lofar_all.colnames)
Out[18]:
The following line has been corrected in the latest versions to use all the sources, including the extended. Hence the running of the "-extended" version of this notebook is no longer necessary.
In [19]:
lofar = field_full.filter_catalogue(lofar_all, colnames=("RA", "DEC"))
In [20]:
combined = field_full.filter_catalogue(combined_all,
colnames=("ra", "dec"))
In [21]:
combined["colour"] = combined["i"] - combined["W1mag"]
In [22]:
combined_aux_index = np.arange(len(combined))
In [23]:
coords_combined = SkyCoord(combined['ra'],
combined['dec'],
unit=(u.deg, u.deg),
frame='icrs')
In [24]:
coords_lofar = SkyCoord(lofar['RA'],
lofar['DEC'],
unit=(u.deg, u.deg),
frame='icrs')
In [25]:
combined_matched = (~np.isnan(combined["i"]) & ~np.isnan(combined["W1mag"])) # Matched i-W1 sources
In [26]:
combined_panstarrs = (~np.isnan(combined["i"]) & np.isnan(combined["W1mag"])) # Sources with only i-band
In [27]:
combined_wise =(np.isnan(combined["i"]) & ~np.isnan(combined["W1mag"])) # Sources with only W1-band
In [28]:
combined_i = combined_matched | combined_panstarrs
combined_w1 = combined_matched | combined_wise
#combined_only_i = combined_panstarrs & ~combined_matched
#combined_only_w1 = combined_wise & ~combined_matched
In [29]:
print("Total - ", len(combined))
print("i and W1 - ", np.sum(combined_matched))
print("Only i - ", np.sum(combined_panstarrs))
print("With i - ", np.sum(combined_i))
print("Only W1 - ", np.sum(combined_wise))
print("With W1 - ", np.sum(combined_w1))
In [30]:
colour_limits = [0.0, 0.5, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, 3.0, 3.5, 4.0]
In [31]:
# Start with the W1-only, i-only and "less than lower colour" bins
colour_bin_def = [{"name":"only W1", "condition": combined_wise},
{"name":"only i", "condition": combined_panstarrs},
{"name":"-inf to {}".format(colour_limits[0]),
"condition": (combined["colour"] < colour_limits[0])}]
# Get the colour bins
for i in range(len(colour_limits)-1):
name = "{} to {}".format(colour_limits[i], colour_limits[i+1])
condition = ((combined["colour"] >= colour_limits[i]) &
(combined["colour"] < colour_limits[i+1]))
colour_bin_def.append({"name":name, "condition":condition})
# Add the "more than higher colour" bin
colour_bin_def.append({"name":"{} to inf".format(colour_limits[-1]),
"condition": (combined["colour"] >= colour_limits[-1])})
In [32]:
combined["category"] = np.nan
for i in range(len(colour_bin_def)):
combined["category"][colour_bin_def[i]["condition"]] = i
In [33]:
np.sum(np.isnan(combined["category"]))
Out[33]:
We get the number of sources of the combined catalogue in each colour category. It will be used at a later stage to compute the $Q_0$ values
In [34]:
numbers_combined_bins = np.array([np.sum(a["condition"]) for a in colour_bin_def])
In [35]:
numbers_combined_bins
Out[35]:
In [36]:
bin_list, centers, Q_0_colour, n_m, q_m = pickle.load(open("lofar_params.pckl", "rb"))
In [37]:
likelihood_ratio_function = MultiMLEstimatorU(Q_0_colour, n_m, q_m, centers)
In [38]:
radius = 15
In [39]:
selection = ~np.isnan(combined["category"]) # Avoid the dreaded sources with no actual data
catalogue = combined[selection]
In [40]:
def apply_ml(i, likelihood_ratio_function):
idx_0 = idx_i[idx_lofar == i]
d2d_0 = d2d[idx_lofar == i]
category = catalogue["category"][idx_0].astype(int)
mag = catalogue["i"][idx_0]
mag[category == 0] = catalogue["W1mag"][idx_0][category == 0]
lofar_ra = lofar[i]["RA"]
lofar_dec = lofar[i]["DEC"]
lofar_pa = lofar[i]["PA"]
lofar_maj_err = lofar[i]["E_Maj"]
lofar_min_err = lofar[i]["E_Min"]
c_ra = catalogue["ra"][idx_0]
c_dec = catalogue["dec"][idx_0]
c_ra_err = catalogue["raErr"][idx_0]
c_dec_err = catalogue["decErr"][idx_0]
sigma, sigma_maj, sigma_min = get_sigma_all(lofar_maj_err, lofar_min_err, lofar_pa,
lofar_ra, lofar_dec,
c_ra, c_dec, c_ra_err, c_dec_err)
lr_0 = likelihood_ratio_function(mag, d2d_0.arcsec, sigma, sigma_maj, sigma_min, category)
chosen_index = np.argmax(lr_0)
result = [combined_aux_index[selection][idx_0[chosen_index]], # Index
(d2d_0.arcsec)[chosen_index], # distance
lr_0[chosen_index]] # LR
return result
In [41]:
idx_lofar, idx_i, d2d, d3d = search_around_sky(
coords_lofar, coords_combined[selection], radius*u.arcsec)
In [42]:
idx_lofar_unique = np.unique(idx_lofar)
In [43]:
def ml(i):
return apply_ml(i, likelihood_ratio_function)
In [44]:
res = parallel_process(idx_lofar_unique, ml, n_jobs=8)
In [45]:
lofar["lr"] = np.nan # Likelihood ratio
lofar["lr_dist"] = np.nan # Distance to the selected source
lofar["lr_index"] = np.nan # Index of the PanSTARRS source in combined
In [46]:
(lofar["lr_index"][idx_lofar_unique],
lofar["lr_dist"][idx_lofar_unique],
lofar["lr"][idx_lofar_unique]) = list(map(list, zip(*res)))
In [47]:
total_sources = len(idx_lofar_unique)
combined_aux_index = np.arange(len(combined))
In [48]:
lofar["lrt"] = lofar["lr"]
lofar["lrt"][np.isnan(lofar["lr"])] = 0
In [49]:
threshold_sel = 0.358
In [50]:
lofar["lr_index_sel"] = lofar["lr_index"]
lofar["lr_index_sel"][lofar["lrt"] < threshold_sel] = np.nan
In [51]:
combined["lr_index_sel"] = combined_aux_index.astype(float)
In [52]:
pwl = join(lofar, combined, join_type='left', keys='lr_index_sel')
In [53]:
pwl_columns = pwl.colnames
In [54]:
for col in pwl_columns:
fv = pwl[col].fill_value
if (isinstance(fv, np.float64) and (fv != 1e+20)):
print(col, fv)
pwl[col].fill_value = 1e+20
In [55]:
columns_save = ['Source_Name', 'RA', 'E_RA', 'E_RA_tot', 'DEC', 'E_DEC',
'E_DEC_tot', 'Peak_flux', 'E_Peak_flux', 'E_Peak_flux_tot',
'Total_flux', 'E_Total_flux', 'E_Total_flux_tot', 'Maj', 'E_Maj',
'Min', 'E_Min', 'DC_Maj', 'E_DC_Maj', 'DC_Min', 'E_DC_Min', 'PA',
'E_PA', 'DC_PA', 'E_DC_PA', 'Isl_rms', 'S_Code', 'Mosaic_ID',
'Isl_id', 'Gaus_id',
'AllWISE', 'objID', 'ra', 'dec', 'raErr', 'decErr',
'W1mag', 'W1magErr', 'i', 'iErr', 'colour', 'category',
'lr', 'lr_dist']
In [56]:
pwl[columns_save].filled().write('lofar_gaus_pw_pdf.fits', format="fits")
In [ ]: