Tibeuleu/FOC_Reduction
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

fix package calling and clean scripts

Browse Source
This commit is contained in:
Thibault Barnouin
2024-09-17 21:07:26 +02:00
parent 10577352d4
commit db3deac6c2
12 changed files with 345 additions and 496 deletions
Download Patch File Download Diff File Expand all files Collapse all files

240
package/src/Combine.py Executable file
View File

@@ -0,0 +1,240 @@
#!/usr/bin/python
# -*- coding:utf-8 -*-
from pathlib import Path
from sys import path as syspath
syspath.append(str(Path(__file__).parent.parent))
import numpy as np
def same_reduction(infiles):
"""
Test if infiles are pipeline productions with same parameters.
"""
from astropy.io.fits import open as fits_open
from astropy.wcs import WCS
params = {"IQU": [], "ROT": [], "SIZE": [], "TARGNAME": [], "BKG_SUB": [], "SAMPLING": [], "SMOOTH": []}
for file in infiles:
with fits_open(file) as f:
# test for presence of I, Q, U images
datatype = []
for hdu in f:
try:
datatype.append(hdu.header["datatype"])
except KeyError:
pass
test_IQU = True
for look in ["I_stokes", "Q_stokes", "U_stokes", "IQU_cov_matrix"]:
test_IQU *= look in datatype
params["IQU"].append(test_IQU)
# test for orientation and pixel size
wcs = WCS(f[0].header).celestial
if wcs.wcs.has_cd() or (wcs.wcs.cdelt[:2] == np.array([1.0, 1.0])).all():
cdelt = np.linalg.eig(wcs.wcs.cd)[0]
pc = np.dot(wcs.wcs.cd, np.diag(1.0 / cdelt))
else:
cdelt = wcs.wcs.cdelt
pc = wcs.wcs.pc
params["ROT"].append(np.round(np.arccos(pc[0, 0]), 2) if np.abs(pc[0, 0]) < 1.0 else 0.0)
params["SIZE"].append(np.round(np.max(np.abs(cdelt * 3600.0)), 2))
# look for information on reduction procedure
for key in [k for k in params.keys() if k not in ["IQU", "ROT", "SIZE"]]:
try:
params[key].append(f[0].header[key])
except KeyError:
params[key].append("null")
result = np.all(params["IQU"])
for key in [k for k in params.keys() if k != "IQU"]:
result *= np.unique(params[key]).size == 1
if np.all(params["IQU"]) and not result:
print(np.unique(params["SIZE"]))
raise ValueError("Not all observations were reduced with the same parameters, please provide the raw files.")
return result
def same_obs(infiles, data_folder):
"""
Group infiles into same observations.
"""
import astropy.units as u
from astropy.io.fits import getheader
from astropy.table import Table
from astropy.time import Time, TimeDelta
headers = [getheader("/".join([data_folder, file])) for file in infiles]
files = {}
files["PROPOSID"] = np.array([str(head["PROPOSID"]) for head in headers], dtype=str)
files["ROOTNAME"] = np.array([head["ROOTNAME"].lower() + "_c0f.fits" for head in headers], dtype=str)
files["EXPSTART"] = np.array([Time(head["EXPSTART"], format="mjd") for head in headers])
products = Table(files)
new_infiles = []
for pid in np.unique(products["PROPOSID"]):
obs = products[products["PROPOSID"] == pid].copy()
close_date = np.unique(
[[np.abs(TimeDelta(obs["EXPSTART"][i].unix - date.unix, format="sec")) < 7.0 * u.d for i in range(len(obs))] for date in obs["EXPSTART"]], axis=0
)
if len(close_date) > 1:
for date in close_date:
new_infiles.append(list(products["ROOTNAME"][np.any([products["ROOTNAME"] == dataset for dataset in obs["ROOTNAME"][date]], axis=0)]))
else:
new_infiles.append(list(products["ROOTNAME"][products["PROPOSID"] == pid]))
return new_infiles
def combine_Stokes(infiles):
"""
Combine I, Q, U from different observations of a same object.
"""
from astropy.io.fits import open as fits_open
from lib.reduction import align_data, zeropad
from scipy.ndimage import shift as sc_shift
I_array, Q_array, U_array, IQU_cov_array, data_mask, headers = [], [], [], [], [], []
shape = np.array([0, 0])
for file in infiles:
with fits_open(file) as f:
headers.append(f[0].header)
I_array.append(f["I_stokes"].data)
Q_array.append(f["Q_stokes"].data)
U_array.append(f["U_stokes"].data)
IQU_cov_array.append(f["IQU_cov_matrix"].data)
data_mask.append(f["data_mask"].data.astype(bool))
shape[0] = np.max([shape[0], f["I_stokes"].data.shape[0]])
shape[1] = np.max([shape[1], f["I_stokes"].data.shape[1]])
exposure_array = np.array([float(head["EXPTIME"]) for head in headers])
shape += np.array([5, 5])
data_mask = np.sum([zeropad(mask, shape) for mask in data_mask], axis=0).astype(bool)
I_array = np.array([zeropad(I, shape) for I in I_array])
Q_array = np.array([zeropad(Q, shape) for Q in Q_array])
U_array = np.array([zeropad(U, shape) for U in U_array])
IQU_cov_array = np.array([[[zeropad(cov[i, j], shape) for j in range(3)] for i in range(3)] for cov in IQU_cov_array])
sI_array = np.sqrt(IQU_cov_array[:, 0, 0])
sQ_array = np.sqrt(IQU_cov_array[:, 1, 1])
sU_array = np.sqrt(IQU_cov_array[:, 2, 2])
_, _, _, _, shifts, errors = align_data(I_array, headers, error_array=sI_array, data_mask=data_mask, ref_center="center", return_shifts=True)
data_mask_aligned = np.sum([sc_shift(data_mask, s, order=1, cval=0.0) for s in shifts], axis=0).astype(bool)
I_aligned, sI_aligned = (
np.array([sc_shift(I, s, order=1, cval=0.0) for I, s in zip(I_array, shifts)]),
np.array([sc_shift(sI, s, order=1, cval=0.0) for sI, s in zip(sI_array, shifts)]),
)
Q_aligned, sQ_aligned = (
np.array([sc_shift(Q, s, order=1, cval=0.0) for Q, s in zip(Q_array, shifts)]),
np.array([sc_shift(sQ, s, order=1, cval=0.0) for sQ, s in zip(sQ_array, shifts)]),
)
U_aligned, sU_aligned = (
np.array([sc_shift(U, s, order=1, cval=0.0) for U, s in zip(U_array, shifts)]),
np.array([sc_shift(sU, s, order=1, cval=0.0) for sU, s in zip(sU_array, shifts)]),
)
IQU_cov_aligned = np.array([[[sc_shift(cov[i, j], s, order=1, cval=0.0) for j in range(3)] for i in range(3)] for cov, s in zip(IQU_cov_array, shifts)])
I_combined = np.sum([exp * I for exp, I in zip(exposure_array, I_aligned)], axis=0) / exposure_array.sum()
Q_combined = np.sum([exp * Q for exp, Q in zip(exposure_array, Q_aligned)], axis=0) / exposure_array.sum()
U_combined = np.sum([exp * U for exp, U in zip(exposure_array, U_aligned)], axis=0) / exposure_array.sum()
IQU_cov_combined = np.zeros((3, 3, shape[0], shape[1]))
for i in range(3):
IQU_cov_combined[i, i] = np.sum([exp**2 * cov for exp, cov in zip(exposure_array, IQU_cov_aligned[:, i, i])], axis=0) / exposure_array.sum() ** 2
for j in [x for x in range(3) if x != i]:
IQU_cov_combined[i, j] = np.sqrt(
np.sum([exp**2 * cov**2 for exp, cov in zip(exposure_array, IQU_cov_aligned[:, i, j])], axis=0) / exposure_array.sum() ** 2
)
IQU_cov_combined[j, i] = np.sqrt(
np.sum([exp**2 * cov**2 for exp, cov in zip(exposure_array, IQU_cov_aligned[:, j, i])], axis=0) / exposure_array.sum() ** 2
)
header_combined = headers[0]
header_combined["EXPTIME"] = exposure_array.sum()
return I_combined, Q_combined, U_combined, IQU_cov_combined, data_mask_aligned, header_combined
def main(infiles, target=None, output_dir="./data/"):
""" """
from lib.fits import save_Stokes
from lib.plots import pol_map
from lib.reduction import compute_pol, rotate_Stokes
if target is None:
target = input("Target name:\n>")
prod = np.array([["/".join(filepath.split("/")[:-1]), filepath.split("/")[-1]] for filepath in infiles], dtype=str)
data_folder = prod[0][0]
files = [p[1] for p in prod]
# Reduction parameters
kwargs = {}
# Polarization map output
kwargs["SNRp_cut"] = 3.0
kwargs["SNRi_cut"] = 1.0
kwargs["flux_lim"] = 1e-19, 3e-17
kwargs["scale_vec"] = 5
kwargs["step_vec"] = 1
if not same_reduction(infiles):
from FOC_reduction import main as FOC_reduction
grouped_infiles = same_obs(files, data_folder)
new_infiles = []
for i, group in enumerate(grouped_infiles):
new_infiles.append(
FOC_reduction(target=target + "-" + str(i + 1), infiles=["/".join([data_folder, file]) for file in group], interactive=True)[0]
)
infiles = new_infiles
I_combined, Q_combined, U_combined, IQU_cov_combined, data_mask_combined, header_combined = combine_Stokes(infiles=infiles)
I_combined, Q_combined, U_combined, IQU_cov_combined, data_mask_combined, header_combined = rotate_Stokes(
I_stokes=I_combined, Q_stokes=Q_combined, U_stokes=U_combined, Stokes_cov=IQU_cov_combined, data_mask=data_mask_combined, header_stokes=header_combined
)
P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P = compute_pol(
I_stokes=I_combined, Q_stokes=Q_combined, U_stokes=U_combined, Stokes_cov=IQU_cov_combined, header_stokes=header_combined
)
filename = header_combined["FILENAME"]
figname = "_".join([target, filename[filename.find("FOC_") :], "combined"])
Stokes_combined = save_Stokes(
I_stokes=I_combined,
Q_stokes=Q_combined,
U_stokes=U_combined,
Stokes_cov=IQU_cov_combined,
P=P,
debiased_P=debiased_P,
s_P=s_P,
s_P_P=s_P_P,
PA=PA,
s_PA=s_PA,
s_PA_P=s_PA_P,
header_stokes=header_combined,
data_mask=data_mask_combined,
filename=figname,
data_folder=data_folder,
return_hdul=True,
)
pol_map(Stokes_combined, **kwargs)
return "/".join([data_folder, figname + ".fits"])
if __name__ == "__main__":
import argparse
parser = argparse.ArgumentParser(description="Combine different observations of a single object")
parser.add_argument("-t", "--target", metavar="targetname", required=False, help="the name of the target", type=str, default=None)
parser.add_argument("-f", "--files", metavar="path", required=False, nargs="*", help="the full or relative path to the data products", default=None)
parser.add_argument(
"-o", "--output_dir", metavar="directory_path", required=False, help="output directory path for the data products", type=str, default="./data"
)
args = parser.parse_args()
exitcode = main(target=args.target, infiles=args.files, output_dir=args.output_dir)
print("Written to: ", exitcode)

0
package/src/__init__.py
View File

40
package/src/analysis.py
View File

@@ -1,40 +0,0 @@
#!/usr/bin/python
from getopt import error as get_error
from getopt import getopt
from sys import argv
arglist = argv[1:]
options = "hf:p:i:l:"
long_options = ["help", "fits=", "snrp=", "snri=", "lim="]
fits_path = None
SNRp_cut, SNRi_cut = 3, 3
flux_lim = None
out_txt = None
try:
arg, val = getopt(arglist, options, long_options)
for curr_arg, curr_val in arg:
if curr_arg in ("-h", "--help"):
print("python3 analysis.py -f <path_to_reduced_fits> -p <SNRp_cut> -i <SNRi_cut> -l <flux_lim>")
elif curr_arg in ("-f", "--fits"):
fits_path = str(curr_val)
elif curr_arg in ("-p", "--snrp"):
SNRp_cut = int(curr_val)
elif curr_arg in ("-i", "--snri"):
SNRi_cut = int(curr_val)
elif curr_arg in ("-l", "--lim"):
flux_lim = list("".join(curr_val).split(","))
except get_error as err:
print(str(err))
if fits_path is not None:
from astropy.io import fits
from lib.plots import pol_map
Stokes_UV = fits.open(fits_path)
p = pol_map(Stokes_UV, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim)
else:
print("python3 analysis.py -f <path_to_reduced_fits> -p <SNRp_cut> -i <SNRi_cut> -l <flux_lim>")

214
package/src/comparison_Kishimoto.py
View File

@@ -1,214 +0,0 @@
#!/usr/bin/python
from src.lib.background import gauss, bin_centers
from src.lib.deconvolve import zeropad
from src.lib.reduction import align_data
from src.lib.plots import princ_angle
from matplotlib.colors import LogNorm
from os.path import join as path_join
from astropy.io import fits
from astropy.wcs import WCS
from scipy.ndimage import shift
from scipy.optimize import curve_fit
import numpy as np
import matplotlib.pyplot as plt
root_dir = path_join('/home/t.barnouin/Documents/Thesis/HST')
root_dir_K = path_join(root_dir, 'Kishimoto', 'output')
root_dir_S = path_join(root_dir, 'FOC_Reduction', 'output')
root_dir_data_S = path_join(root_dir, 'FOC_Reduction', 'data', 'NGC1068', '5144')
root_dir_plot_S = path_join(root_dir, 'FOC_Reduction', 'plots', 'NGC1068', '5144', 'notaligned')
filename_S = "NGC1068_FOC_b10.00pixel_not_aligned.fits"
plt.rcParams.update({'font.size': 15})
SNRi_cut = 30.
SNRp_cut = 3.
data_K = {}
data_S = {}
for d, i in zip(['I', 'Q', 'U', 'P', 'PA', 'sI', 'sQ', 'sU', 'sP', 'sPA'], [0, 1, 2, 5, 8, (3, 0, 0), (3, 1, 1), (3, 2, 2), 6, 9]):
data_K[d] = np.loadtxt(path_join(root_dir_K, d+'.txt'))
with fits.open(path_join(root_dir_data_S, filename_S)) as f:
if not type(i) is int:
data_S[d] = np.sqrt(f[i[0]].data[i[1], i[2]])
else:
data_S[d] = f[i].data
if i == 0:
header = f[i].header
wcs = WCS(header)
convert_flux = header['photflam']
bkg_S = np.median(data_S['I'])/3
bkg_K = np.median(data_K['I'])/3
# zeropad data to get same size of array
shape = data_S['I'].shape
for d in data_K:
data_K[d] = zeropad(data_K[d], shape)
# shift array to get same information in same pixel
data_arr, error_ar, heads, data_msk, shifts, shifts_err = align_data(np.array([data_S['I'], data_K['I']]), [header, header], error_array=np.array(
[data_S['sI'], data_K['sI']]), background=np.array([bkg_S, bkg_K]), upsample_factor=10., ref_center='center', return_shifts=True)
for d in data_K:
data_K[d] = shift(data_K[d], shifts[1], order=1, cval=0.)
# compute pol components from shifted array
for d in [data_S, data_K]:
for i in d:
d[i][np.isnan(d[i])] = 0.
d['P'] = np.where(np.logical_and(np.isfinite(d['I']), d['I'] > 0.), np.sqrt(d['Q']**2+d['U']**2)/d['I'], 0.)
d['sP'] = np.where(np.logical_and(np.isfinite(d['I']), d['I'] > 0.), np.sqrt((d['Q']**2*d['sQ']**2+d['U']**2*d['sU']**2) /
(d['Q']**2+d['U']**2)+((d['Q']/d['I'])**2+(d['U']/d['I'])**2)*d['sI']**2)/d['I'], 0.)
d['d_P'] = np.where(np.logical_and(np.isfinite(d['P']), np.isfinite(d['sP'])), np.sqrt(d['P']**2-d['sP']**2), 0.)
d['PA'] = 0.5*np.arctan2(d['U'], d['Q'])+np.pi
d['SNRp'] = np.zeros(d['d_P'].shape)
d['SNRp'][d['sP'] > 0.] = d['d_P'][d['sP'] > 0.]/d['sP'][d['sP'] > 0.]
d['SNRi'] = np.zeros(d['I'].shape)
d['SNRi'][d['sI'] > 0.] = d['I'][d['sI'] > 0.]/d['sI'][d['sI'] > 0.]
d['mask'] = np.logical_and(d['SNRi'] > SNRi_cut, d['SNRp'] > SNRp_cut)
data_S['mask'], data_K['mask'] = np.logical_and(data_S['mask'], data_K['mask']), np.logical_and(data_S['mask'], data_K['mask'])
#
# Compute histogram of measured polarization in cut
#
bins = int(data_S['mask'].sum()/5)
bin_size = 1./bins
mod_p = np.linspace(0., 1., 300)
for d in [data_S, data_K]:
d['hist'], d['bin_edges'] = np.histogram(d['d_P'][d['mask']], bins=bins, range=(0., 1.))
d['binning'] = bin_centers(d['bin_edges'])
peak, bins_fwhm = d['binning'][np.argmax(d['hist'])], d['binning'][d['hist'] > d['hist'].max()/2.]
fwhm = bins_fwhm[1]-bins_fwhm[0]
p0 = [d['hist'].max(), peak, fwhm]
try:
popt, pcov = curve_fit(gauss, d['binning'], d['hist'], p0=p0)
except RuntimeError:
popt = p0
d['hist_chi2'] = np.sum((d['hist']-gauss(d['binning'], *popt))**2)/d['hist'].size
d['hist_popt'] = popt
fig_p, ax_p = plt.subplots(num="Polarization degree histogram", figsize=(10, 6), constrained_layout=True)
ax_p.errorbar(data_S['binning'], data_S['hist'], xerr=bin_size/2., fmt='b.', ecolor='b', label='P through this pipeline')
ax_p.plot(mod_p, gauss(mod_p, *data_S['hist_popt']), 'b--', label='mean = {1:.2f}, stdev = {2:.2f}'.format(*data_S['hist_popt']))
ax_p.errorbar(data_K['binning'], data_K['hist'], xerr=bin_size/2., fmt='r.', ecolor='r', label="P through Kishimoto's pipeline")
ax_p.plot(mod_p, gauss(mod_p, *data_K['hist_popt']), 'r--', label='mean = {1:.2f}, stdev = {2:.2f}'.format(*data_K['hist_popt']))
ax_p.set(xlabel="Polarization degree", ylabel="Counts", title="Histogram of polarization degree computed in the cut for both pipelines.")
ax_p.legend()
fig_p.savefig(path_join(root_dir_plot_S, "NGC1068_K_pol_deg.png"), bbox_inches="tight", dpi=300)
#
# Compute angular difference between the maps in cut
#
dtheta = np.where(data_S['mask'], 0.5*np.arctan((np.sin(2*data_S['PA'])*np.cos(2*data_K['PA'])-np.cos(2*data_S['PA']) *
np.cos(2*data_K['PA']))/(np.cos(2*data_S['PA'])*np.cos(2*data_K['PA'])+np.cos(2*data_S['PA'])*np.sin(2*data_K['PA']))), np.nan)
fig_pa = plt.figure(num="Polarization degree alignement")
ax_pa = fig_pa.add_subplot(111, projection=wcs)
cbar_ax_pa = fig_pa.add_axes([0.88, 0.12, 0.01, 0.75])
ax_pa.set_title(r"Degree of alignement $\zeta$ of the polarization angles from the 2 pipelines in the cut")
im_pa = ax_pa.imshow(np.cos(2*dtheta), vmin=-1., vmax=1., origin='lower', cmap='bwr', label=r"$\zeta$ between this pipeline and Kishimoto's")
cbar_pa = plt.colorbar(im_pa, cax=cbar_ax_pa, label=r"$\zeta = \cos\left( 2 \cdot \delta\theta_P \right)$")
ax_pa.coords[0].set_axislabel('Right Ascension (J2000)')
ax_pa.coords[1].set_axislabel('Declination (J2000)')
fig_pa.savefig(path_join(root_dir_plot_S, "NGC1068_K_pol_ang.png"), bbox_inches="tight", dpi=300)
#
# Compute power uncertainty difference between the maps in cut
#
eta = np.where(data_S['mask'], np.abs(data_K['d_P']-data_S['d_P'])/np.sqrt(data_S['sP']**2+data_K['sP']**2)/2., np.nan)
fig_dif_p = plt.figure(num="Polarization power difference ratio")
ax_dif_p = fig_dif_p.add_subplot(111, projection=wcs)
cbar_ax_dif_p = fig_dif_p.add_axes([0.88, 0.12, 0.01, 0.75])
ax_dif_p.set_title(r"Degree of difference $\eta$ of the polarization from the 2 pipelines in the cut")
im_dif_p = ax_dif_p.imshow(eta, vmin=0., vmax=2., origin='lower', cmap='bwr_r', label=r"$\eta$ between this pipeline and Kishimoto's")
cbar_dif_p = plt.colorbar(im_dif_p, cax=cbar_ax_dif_p, label=r"$\eta = \frac{2 \left|P^K-P^S\right|}{\sqrt{{\sigma^K_P}^2+{\sigma^S_P}^2}}$")
ax_dif_p.coords[0].set_axislabel('Right Ascension (J2000)')
ax_dif_p.coords[1].set_axislabel('Declination (J2000)')
fig_dif_p.savefig(path_join(root_dir_plot_S, "NGC1068_K_pol_diff.png"), bbox_inches="tight", dpi=300)
#
# Compute angle uncertainty difference between the maps in cut
#
eta = np.where(data_S['mask'], np.abs(data_K['PA']-data_S['PA'])/np.sqrt(data_S['sPA']**2+data_K['sPA']**2)/2., np.nan)
fig_dif_pa = plt.figure(num="Polarization angle difference ratio")
ax_dif_pa = fig_dif_pa.add_subplot(111, projection=wcs)
cbar_ax_dif_pa = fig_dif_pa.add_axes([0.88, 0.12, 0.01, 0.75])
ax_dif_pa.set_title(r"Degree of difference $\eta$ of the polarization from the 2 pipelines in the cut")
im_dif_pa = ax_dif_pa.imshow(eta, vmin=0., vmax=2., origin='lower', cmap='bwr_r', label=r"$\eta$ between this pipeline and Kishimoto's")
cbar_dif_pa = plt.colorbar(im_dif_pa, cax=cbar_ax_dif_pa,
label=r"$\eta = \frac{2 \left|\theta_P^K-\theta_P^S\right|}{\sqrt{{\sigma^K_{\theta_P}}^2+{\sigma^S_{\theta_P}}^2}}$")
ax_dif_pa.coords[0].set_axislabel('Right Ascension (J2000)')
ax_dif_pa.coords[1].set_axislabel('Declination (J2000)')
fig_dif_pa.savefig(path_join(root_dir_plot_S, "NGC1068_K_polang_diff.png"), bbox_inches="tight", dpi=300)
# display both polarization maps to check consistency
# plt.rcParams.update({'font.size': 15})
fig = plt.figure(num="Polarization maps comparison", figsize=(10, 10))
ax = fig.add_subplot(111, projection=wcs)
fig.subplots_adjust(right=0.85)
cbar_ax = fig.add_axes([0.88, 0.12, 0.01, 0.75])
for d in [data_S, data_K]:
d['X'], d['Y'] = np.meshgrid(np.arange(d['I'].shape[1]), np.arange(d['I'].shape[0]))
d['xy_U'], d['xy_V'] = np.where(d['mask'], d['d_P']*np.cos(np.pi/2.+d['PA']), np.nan), np.where(d['mask'], d['d_P']*np.sin(np.pi/2.+d['PA']), np.nan)
im0 = ax.imshow(data_S['I']*convert_flux, norm=LogNorm(data_S['I'][data_S['I'] > 0].min()*convert_flux, data_S['I']
[data_S['I'] > 0].max()*convert_flux), origin='lower', cmap='gray', label=r"$I_{STOKES}$ through this pipeline")
quiv0 = ax.quiver(data_S['X'], data_S['Y'], data_S['xy_U'], data_S['xy_V'], units='xy', angles='uv', scale=0.5, scale_units='xy',
pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.2, color='b', alpha=0.75, label="PA through this pipeline")
quiv1 = ax.quiver(data_K['X'], data_K['Y'], data_K['xy_U'], data_K['xy_V'], units='xy', angles='uv', scale=0.5, scale_units='xy',
pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.1, color='r', alpha=0.75, label="PA through Kishimoto's pipeline")
ax.set_title(r"$SNR_P \geq$ "+str(SNRi_cut)+r"$\; & \; SNR_I \geq $"+str(SNRp_cut))
# ax.coords.grid(True, color='white', ls='dotted', alpha=0.5)
ax.coords[0].set_axislabel('Right Ascension (J2000)')
ax.coords[0].set_axislabel_position('b')
ax.coords[0].set_ticklabel_position('b')
ax.coords[1].set_axislabel('Declination (J2000)')
ax.coords[1].set_axislabel_position('l')
ax.coords[1].set_ticklabel_position('l')
# ax.axis('equal')
cbar = plt.colorbar(im0, cax=cbar_ax, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]")
ax.legend(loc='upper right')
fig.savefig(path_join(root_dir_plot_S, "NGC1068_K_comparison.png"), bbox_inches="tight", dpi=300)
# compute integrated polarization parameters on a specific cut
for d in [data_S, data_K]:
d['I_dil'] = np.sum(d['I'][d['mask']])
d['sI_dil'] = np.sqrt(np.sum(d['sI'][d['mask']]**2))
d['Q_dil'] = np.sum(d['Q'][d['mask']])
d['sQ_dil'] = np.sqrt(np.sum(d['sQ'][d['mask']]**2))
d['U_dil'] = np.sum(d['U'][d['mask']])
d['sU_dil'] = np.sqrt(np.sum(d['sU'][d['mask']]**2))
d['P_dil'] = np.sqrt(d['Q_dil']**2+d['U_dil']**2)/d['I_dil']
d['sP_dil'] = np.sqrt((d['Q_dil']**2*d['sQ_dil']**2+d['U_dil']**2*d['sU_dil']**2)/(d['Q_dil']**2+d['U_dil']**2) +
((d['Q_dil']/d['I_dil'])**2+(d['U_dil']/d['I_dil'])**2)*d['sI_dil']**2)/d['I_dil']
d['d_P_dil'] = np.sqrt(d['P_dil']**2-d['sP_dil']**2)
d['PA_dil'] = princ_angle((90./np.pi)*np.arctan2(d['U_dil'], d['Q_dil']))
d['sPA_dil'] = princ_angle((90./(np.pi*(d['Q_dil']**2+d['U_dil']**2)))*np.sqrt(d['Q_dil']**2*d['sU_dil']**2+d['U_dil']**2*d['sU_dil']**2))
print('From this pipeline :\n', "P = {0:.2f} ± {1:.2f} %\n".format(
data_S['d_P_dil']*100., data_S['sP_dil']*100.), "PA = {0:.2f} ± {1:.2f} °".format(data_S['PA_dil'], data_S['sPA_dil']))
print("From Kishimoto's pipeline :\n", "P = {0:.2f} ± {1:.2f} %\n".format(
data_K['d_P_dil']*100., data_K['sP_dil']*100.), "PA = {0:.2f} ± {1:.2f} °".format(data_K['PA_dil'], data_K['sPA_dil']))
# compare different types of error
print("This pipeline : average sI/I={0:.2f} ; sQ/Q={1:.2f} ; sU/U={2:.2f} ; sP/P={3:.2f}".format(np.mean(data_S['sI'][data_S['mask']]/data_S['I'][data_S['mask']]), np.mean(
data_S['sQ'][data_S['mask']]/data_S['Q'][data_S['mask']]), np.mean(data_S['sU'][data_S['mask']]/data_S['U'][data_S['mask']]), np.mean(data_S['sP'][data_S['mask']]/data_S['P'][data_S['mask']])))
print("Kishimoto's pipeline : average sI/I={0:.2f} ; sQ/Q={1:.2f} ; sU/U={2:.2f} ; sP/P={3:.2f}".format(np.mean(data_K['sI'][data_S['mask']]/data_K['I'][data_S['mask']]), np.mean(
data_K['sQ'][data_S['mask']]/data_K['Q'][data_S['mask']]), np.mean(data_K['sU'][data_S['mask']]/data_K['U'][data_S['mask']]), np.mean(data_K['sP'][data_S['mask']]/data_K['P'][data_S['mask']])))
for d, i in zip(['I', 'Q', 'U', 'P', 'PA', 'sI', 'sQ', 'sU', 'sP', 'sPA'], [0, 1, 2, 5, 8, (3, 0, 0), (3, 1, 1), (3, 2, 2), 6, 9]):
data_K[d] = np.loadtxt(path_join(root_dir_K, d+'.txt'))
with fits.open(path_join(root_dir_data_S, filename_S)) as f:
if not type(i) is int:
data_S[d] = np.sqrt(f[i[0]].data[i[1], i[2]])
else:
data_S[d] = f[i].data
if i == 0:
header = f[i].header
# from Kishimoto's pipeline : IQU_dir, IQU_shift, IQU_stat, IQU_trans
# from my pipeline : raw_bg, raw_flat, raw_psf, raw_shift, raw_wav, IQU_dir
# but errors from my pipeline are propagated all along, how to compare then ?
plt.show()

77
package/src/emission_center.py Executable file
View File

@@ -0,0 +1,77 @@
#!/usr/bin/python
# -*- coding:utf-8 -*-
from pathlib import Path
from sys import path as syspath
syspath.append(str(Path(__file__).parent.parent))
def main(infile, target=None, output_dir=None):
from os.path import join as pathjoin
import numpy as np
from astropy.io.fits import open as fits_open
from astropy.wcs import WCS
from lib.plots import polarization_map
from lib.utils import CenterConf, PCconf
from matplotlib.patches import Rectangle
from matplotlib.pyplot import show
output = []
levelssnr = np.array([3.0, 4.0])
levelsconf = np.array([0.99])
Stokes = fits_open(infile)
stkI = Stokes["I_STOKES"].data
QN, UN, QN_ERR, UN_ERR = np.full((4, stkI.shape[0], stkI.shape[1]), np.nan)
for sflux, nflux in zip(
[Stokes["Q_STOKES"].data, Stokes["U_STOKES"].data, np.sqrt(Stokes["IQU_COV_MATRIX"].data[1, 1]), np.sqrt(Stokes["IQU_COV_MATRIX"].data[2, 2])],
[QN, UN, QN_ERR, UN_ERR],
):
nflux[stkI > 0.0] = sflux[stkI > 0.0] / stkI[stkI > 0.0]
Stokesconf = PCconf(QN, UN, QN_ERR, UN_ERR)
Stokesmask = Stokes["DATA_MASK"].data.astype(bool)
Stokessnr = np.zeros(Stokesmask.shape)
Stokessnr[Stokes["POL_DEG_ERR"].data > 0.0] = (
Stokes["POL_DEG_DEBIASED"].data[Stokes["POL_DEG_ERR"].data > 0.0] / Stokes["POL_DEG_ERR"].data[Stokes["POL_DEG_ERR"].data > 0.0]
)
Stokescentconf, Stokescenter = CenterConf(Stokesconf > 0.99, Stokes["POL_ANG"].data, Stokes["POL_ANG_ERR"].data)
Stokespos = WCS(Stokes[0].header).pixel_to_world(*Stokescenter)
if target is None:
target = Stokes[0].header["TARGNAME"]
fig, ax = polarization_map(Stokes, P_cut=0.99, step_vec=2, scale_vec=5, display="i")
snrcont = ax.contour(Stokessnr, levelssnr, colors="b")
confcont = ax.contour(Stokesconf, levelsconf, colors="r")
confcenter = ax.plot(*Stokescenter, marker="+", color="gray", label="Best confidence for center: {0}".format(Stokespos.to_string("hmsdms")))
confcentcont = ax.contour(Stokescentconf, [0.01], colors="gray")
handles, labels = ax.get_legend_handles_labels()
labels.append(r"$SNR_P \geq$ 3 and 4 contours")
handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=snrcont.get_edgecolor()[0]))
labels.append(r"Polarization $Conf_{99\%}$ contour")
handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=confcont.get_edgecolor()[0]))
labels.append(r"Center $Conf_{99\%}$ contour")
handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=confcentcont.get_edgecolor()[0]))
ax.legend(handles=handles, labels=labels, bbox_to_anchor=(0.0, -0.12, 1.0, 0.102), loc="lower left", mode="expand", borderaxespad=0.0)
show()
if output_dir is not None:
filename = pathjoin(output_dir, "%s_center.pdf" % target)
fig.savefig(filename, dpi=150, facecolor="None")
output.append(filename)
return output
if __name__ == "__main__":
import argparse
parser = argparse.ArgumentParser(description="Look for the center of emission for a given reduced observation")
parser.add_argument("-t", "--target", metavar="targetname", required=False, help="the name of the target", type=str, default=None)
parser.add_argument("-f", "--file", metavar="path", required=False, help="The full or relative path to the data product", type=str, default=None)
parser.add_argument("-o", "--output_dir", metavar="directory_path", required=False, help="output directory path for the plots", type=str, default="./data")
args = parser.parse_args()
exitcode = main(infile=args.file, target=args.target, output_dir=args.output_dir)
print("Written to: ", exitcode)

5
package/src/get_cdelt.py
View File

@@ -1,4 +1,9 @@
#!/usr/bin/python
# -*- coding:utf-8 -*-
from pathlib import Path
from sys import path as syspath
syspath.append(str(Path(__file__).parent.parent))
def main(infiles=None):

56
package/src/overplot_IC5063.py Executable file
View File

@@ -0,0 +1,56 @@
#!/usr/bin/python3
# -*- coding:utf-8 -*-
from pathlib import Path
from sys import path as syspath
syspath.append(str(Path(__file__).parent.parent))
import numpy as np
from astropy.io import fits
from lib.plots import overplot_pol, overplot_radio
from matplotlib.colors import LogNorm
Stokes_UV = fits.open("./data/IC5063/5918/IC5063_FOC_b0.10arcsec_c0.20arcsec.fits")
Stokes_18GHz = fits.open("./data/IC5063/radio/IC5063_18GHz.fits")
Stokes_24GHz = fits.open("./data/IC5063/radio/IC5063_24GHz.fits")
Stokes_103GHz = fits.open("./data/IC5063/radio/IC5063_103GHz.fits")
Stokes_229GHz = fits.open("./data/IC5063/radio/IC5063_229GHz.fits")
Stokes_357GHz = fits.open("./data/IC5063/radio/IC5063_357GHz.fits")
# Stokes_S2 = fits.open("./data/IC5063/POLARIZATION_COMPARISON/S2_rot_crop.fits")
Stokes_IR = fits.open("./data/IC5063/IR/u2e65g01t_c0f_rot.fits")
# levelsMorganti = np.array([1.,2.,3.,8.,16.,32.,64.,128.])
levelsMorganti = np.logspace(-0.1249, 1.97, 7) / 100.0
levels18GHz = levelsMorganti * Stokes_18GHz[0].data.max()
A = overplot_radio(Stokes_UV, Stokes_18GHz)
A.plot(levels=levels18GHz, P_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/18GHz_overplot.pdf", scale_vec=None)
levels24GHz = levelsMorganti * Stokes_24GHz[0].data.max()
B = overplot_radio(Stokes_UV, Stokes_24GHz)
B.plot(levels=levels24GHz, P_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/24GHz_overplot.pdf", scale_vec=None)
levels103GHz = levelsMorganti * Stokes_103GHz[0].data.max()
C = overplot_radio(Stokes_UV, Stokes_103GHz)
C.plot(levels=levels103GHz, P_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/103GHz_overplot.pdf", scale_vec=None)
levels229GHz = levelsMorganti * Stokes_229GHz[0].data.max()
D = overplot_radio(Stokes_UV, Stokes_229GHz)
D.plot(levels=levels229GHz, P_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/229GHz_overplot.pdf", scale_vec=None)
levels357GHz = levelsMorganti * Stokes_357GHz[0].data.max()
E = overplot_radio(Stokes_UV, Stokes_357GHz)
E.plot(levels=levels357GHz, P_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/357GHz_overplot.pdf", scale_vec=None)
# F = overplot_pol(Stokes_UV, Stokes_S2)
# F.plot(P_cut=3.0, SNRi_cut=80.0, savename='./plots/IC5063/S2_overplot.pdf', norm=LogNorm(vmin=5e-20,vmax=5e-18))
G = overplot_pol(Stokes_UV, Stokes_IR, cmap="inferno")
G.plot(
P_cut=2.0,
SNRi_cut=10.0,
savename="./plots/IC5063/IR_overplot.pdf",
scale_vec=None,
norm=LogNorm(Stokes_IR[0].data.max() * Stokes_IR[0].header["photflam"] / 1e3, Stokes_IR[0].data.max() * Stokes_IR[0].header["photflam"]),
cmap="inferno_r",
)

26
package/src/overplot_MRK463E.py Executable file
View File

@@ -0,0 +1,26 @@
#!/usr/bin/python3
# -*- coding:utf-8 -*-
from pathlib import Path
from sys import path as syspath
syspath.append(str(Path(__file__).parent.parent))
import numpy as np
from astropy.io import fits
from lib.plots import overplot_chandra, overplot_pol
from matplotlib.colors import LogNorm
Stokes_UV = fits.open("./data/MRK463E/5960/MRK463E_FOC_b0.05arcsec_c0.07arcsec.fits")
Stokes_IR = fits.open("./data/MRK463E/WFPC2/IR_rot_crop.fits")
Stokes_Xr = fits.open("./data/MRK463E/Chandra/X_ray_crop.fits")
levels = np.geomspace(1.0, 99.0, 7)
A = overplot_chandra(Stokes_UV, Stokes_Xr, norm=LogNorm())
A.plot(levels=levels, P_cut=0.99, SNRi_cut=1.0, scale_vec=5, zoom=1, savename="./plots/MRK463E/Chandra_overplot.pdf")
A.write_to(path1="./data/MRK463E/FOC_data_Chandra.fits", path2="./data/MRK463E/Chandra_data.fits", suffix="aligned")
levels = np.array([0.8, 2, 5, 10, 20, 50]) / 100.0 * Stokes_UV[0].header["photflam"]
B = overplot_pol(Stokes_UV, Stokes_IR, norm=LogNorm())
B.plot(levels=levels, P_cut=0.99, SNRi_cut=1.0, scale_vec=5, norm=LogNorm(8.5e-18, 2.5e-15), savename="./plots/MRK463E/IR_overplot.pdf")
B.write_to(path1="./data/MRK463E/FOC_data_WFPC.fits", path2="./data/MRK463E/WFPC_data.fits", suffix="aligned")