From db3deac6c268ef2974e2ae8ad45531f380cee328 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Tue, 17 Sep 2024 21:07:26 +0200 Subject: [PATCH] fix package calling and clean scripts --- package/FOC_reduction.py | 40 +++- package/__init__.py | 3 +- package/lib/plots.py | 314 +++++++++++++++++--------- package/{ => src}/Combine.py | 5 +- package/src/__init__.py | 0 package/src/analysis.py | 40 ---- package/src/comparison_Kishimoto.py | 214 ------------------ package/src/emission_center.py | 77 +++++++ package/src/get_cdelt.py | 5 + package/{ => src}/overplot_IC5063.py | 22 +- package/{ => src}/overplot_MRK463E.py | 12 +- package/test_center.py | 109 --------- 12 files changed, 345 insertions(+), 496 deletions(-) rename package/{ => src}/Combine.py (98%) delete mode 100644 package/src/__init__.py delete mode 100755 package/src/analysis.py delete mode 100755 package/src/comparison_Kishimoto.py create mode 100755 package/src/emission_center.py rename package/{ => src}/overplot_IC5063.py (66%) rename package/{ => src}/overplot_MRK463E.py (66%) delete mode 100644 package/test_center.py diff --git a/package/FOC_reduction.py b/package/FOC_reduction.py index 04345f3..58d73c2 100755 --- a/package/FOC_reduction.py +++ b/package/FOC_reduction.py @@ -3,6 +3,10 @@ """ Main script where are progressively added the steps for the FOC pipeline reduction. """ +from pathlib import Path +from sys import path as syspath + +syspath.append(str(Path(__file__).parent.parent)) # Project libraries from copy import deepcopy @@ -61,7 +65,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= rotate_North = True # Polarization map output - SNRp_cut = 3.0 # P measurments with SNR>3 + P_cut = 0.99 # if >=1.0 cut on the signal-to-noise else cut on the confidence level in Q, U SNRi_cut = 1.0 # I measurments with SNR>30, which implies an uncertainty in P of 4.7%. flux_lim = None # lowest and highest flux displayed on plot, defaults to bkg and maximum in cut if None scale_vec = 3 @@ -292,7 +296,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= proj_plots.polarization_map( deepcopy(Stokes_hdul), data_mask, - SNRp_cut=SNRp_cut, + P_cut=P_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, @@ -303,7 +307,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= proj_plots.polarization_map( deepcopy(Stokes_hdul), data_mask, - SNRp_cut=SNRp_cut, + P_cut=P_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, @@ -315,7 +319,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= proj_plots.polarization_map( deepcopy(Stokes_hdul), data_mask, - SNRp_cut=SNRp_cut, + P_cut=P_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, @@ -327,7 +331,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= proj_plots.polarization_map( deepcopy(Stokes_hdul), data_mask, - SNRp_cut=SNRp_cut, + P_cut=P_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, @@ -339,7 +343,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= proj_plots.polarization_map( deepcopy(Stokes_hdul), data_mask, - SNRp_cut=SNRp_cut, + P_cut=P_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, @@ -351,7 +355,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= proj_plots.polarization_map( deepcopy(Stokes_hdul), data_mask, - SNRp_cut=SNRp_cut, + P_cut=P_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, @@ -363,7 +367,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= proj_plots.polarization_map( deepcopy(Stokes_hdul), data_mask, - SNRp_cut=SNRp_cut, + P_cut=P_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, @@ -375,7 +379,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= proj_plots.polarization_map( deepcopy(Stokes_hdul), data_mask, - SNRp_cut=SNRp_cut, + P_cut=P_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, @@ -387,7 +391,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= proj_plots.polarization_map( deepcopy(Stokes_hdul), data_mask, - SNRp_cut=SNRp_cut, + P_cut=P_cut if P_cut >= 1. else 3., SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, @@ -396,12 +400,24 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= plots_folder=plots_folder, display="SNRp", ) + proj_plots.polarization_map( + deepcopy(Stokes_hdul), + data_mask, + P_cut=P_cut if P_cut < 1. else 0.99, + SNRi_cut=SNRi_cut, + flux_lim=flux_lim, + step_vec=step_vec, + scale_vec=scale_vec, + savename="_".join([figname, "confP"]), + plots_folder=plots_folder, + display="confp", + ) elif not interactive: proj_plots.polarization_map( - deepcopy(Stokes_hdul), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, savename=figname, plots_folder=plots_folder, display="integrate" + deepcopy(Stokes_hdul), data_mask, P_cut=P_cut, SNRi_cut=SNRi_cut, savename=figname, plots_folder=plots_folder, display="integrate" ) elif pxscale.lower() not in ["full", "integrate"]: - proj_plots.pol_map(Stokes_hdul, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, scale_vec=scale_vec, flux_lim=flux_lim) + proj_plots.pol_map(Stokes_hdul, P_cut=P_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, scale_vec=scale_vec, flux_lim=flux_lim) return outfiles diff --git a/package/__init__.py b/package/__init__.py index 094aa13..a8ca1bd 100644 --- a/package/__init__.py +++ b/package/__init__.py @@ -1,3 +1,2 @@ -from . import lib -from . import src +from .lib import * from . import FOC_reduction diff --git a/package/lib/plots.py b/package/lib/plots.py index 3a694c6..df6354f 100755 --- a/package/lib/plots.py +++ b/package/lib/plots.py @@ -9,7 +9,7 @@ prototypes : - plot_Stokes(Stokes, savename, plots_folder) Plot the I/Q/U maps from the Stokes HDUList. - - polarization_map(Stokes, data_mask, rectangle, SNRp_cut, SNRi_cut, step_vec, savename, plots_folder, display) -> fig, ax + - polarization_map(Stokes, data_mask, rectangle, P_cut, SNRi_cut, step_vec, savename, plots_folder, display) -> fig, ax Plots polarization map of polarimetric parameters saved in an HDUList. class align_maps(map, other_map, **kwargs) @@ -36,7 +36,7 @@ prototypes : class aperture(img, cdelt, radius, fig, ax) Class to interactively simulate aperture integration. - class pol_map(Stokes, SNRp_cut, SNRi_cut, selection) + class pol_map(Stokes, P_cut, SNRi_cut, selection) Class to interactively study polarization maps making use of the cropping and selecting tools. """ @@ -60,10 +60,7 @@ from mpl_toolkits.axes_grid1.anchored_artists import ( ) from scipy.ndimage import zoom as sc_zoom -try: - from .utils import PCconf, princ_angle, rot2D, sci_not -except ImportError: - from utils import PCconf, princ_angle, rot2D, sci_not +from .utils import PCconf, princ_angle, rot2D, sci_not def plot_obs(data_array, headers, rectangle=None, savename=None, plots_folder="", **kwargs): @@ -216,11 +213,11 @@ def polarization_map( Stokes, data_mask=None, rectangle=None, - SNRp_cut=3.0, - SNRi_cut=3.0, + P_cut=0.99, + SNRi_cut=1.0, flux_lim=None, step_vec=1, - scale_vec=2.0, + scale_vec=3.0, savename=None, plots_folder="", display="default", @@ -236,9 +233,9 @@ def polarization_map( Array of parameters for matplotlib.patches.Rectangle objects that will be displayed on each output image. If None, no rectangle displayed. Defaults to None. - SNRp_cut : float, optional + P_cut : float, optional Cut that should be applied to the signal-to-noise ratio on P. - Any SNR < SNRp_cut won't be displayed. + Any SNR < P_cut won't be displayed. Defaults to 3. SNRi_cut : float, optional Cut that should be applied to the signal-to-noise ratio on I. @@ -276,15 +273,24 @@ def polarization_map( """ # Get data stkI = Stokes["I_stokes"].data.copy() + stkQ = Stokes["Q_stokes"].data.copy() + stkU = Stokes["U_stokes"].data.copy() stk_cov = Stokes["IQU_cov_matrix"].data.copy() pol = Stokes["Pol_deg_debiased"].data.copy() pol_err = Stokes["Pol_deg_err"].data.copy() pang = Stokes["Pol_ang"].data.copy() - try: - if data_mask is None: + if data_mask is None: + try: data_mask = Stokes["Data_mask"].data.astype(bool).copy() - except KeyError: - data_mask = np.ones(stkI.shape).astype(bool) + except KeyError: + data_mask = np.zeros(stkI.shape).astype(bool) + data_mask[stkI > 0.0] = True + + # Compute confidence level map + QN, UN, QN_ERR, UN_ERR = np.full((4, stkI.shape[0], stkI.shape[1]), np.nan) + for nflux, sflux in zip([QN, UN, QN_ERR, UN_ERR], [stkQ, stkU, np.sqrt(stk_cov[1, 1]), np.sqrt(stk_cov[2, 2])]): + nflux[stkI > 0.0] = sflux[stkI > 0.0] / stkI[stkI > 0.0] + confP = PCconf(QN, UN, QN_ERR, UN_ERR) for dataset in [stkI, pol, pol_err, pang]: dataset[np.logical_not(data_mask)] = np.nan @@ -302,15 +308,23 @@ def polarization_map( # Compute SNR and apply cuts poldata, pangdata = pol.copy(), pang.copy() - maskP = pol_err > 0 - SNRp = np.ones(pol.shape) * np.nan - SNRp[maskP] = pol[maskP] / pol_err[maskP] + SNRi = np.full(stkI.shape, np.nan) + SNRi[stk_cov[0, 0] > 0.0] = stkI[stk_cov[0, 0] > 0.0] / np.sqrt(stk_cov[0, 0][stk_cov[0, 0] > 0.0]) + maskI = np.zeros(stkI.shape, dtype=bool) + maskI[stk_cov[0, 0] > 0.0] = SNRi[stk_cov[0, 0] > 0.0] > SNRi_cut - maskI = stk_cov[0, 0] > 0 - SNRi = np.ones(stkI.shape) * np.nan - SNRi[maskI] = stkI[maskI] / np.sqrt(stk_cov[0, 0][maskI]) + SNRp = np.full(pol.shape, np.nan) + SNRp[pol_err > 0.0] = pol[pol_err > 0.0] / pol_err[pol_err > 0.0] + maskP = np.zeros(pol.shape, dtype=bool) - mask = (SNRp > SNRp_cut) * (SNRi > SNRi_cut) + if P_cut >= 1.0: + # MaskP on the signal-to-noise value + maskP[pol_err > 0.0] = SNRp[pol_err > 0.0] > P_cut + else: + # MaskP on the confidence value + maskP = confP > P_cut + + mask = np.logical_and(maskI, maskP) poldata[np.logical_not(mask)] = np.nan pangdata[np.logical_not(mask)] = np.nan @@ -381,8 +395,8 @@ def polarization_map( elif display.lower() in ["s_p", "pol_err", "pol_deg_err"]: # Display polarization degree error map display = "s_p" - if (SNRp > SNRp_cut).any(): - vmin, vmax = 0.0, np.max([pol_err[SNRp > SNRp_cut].max(), 1.0]) * 100.0 + if (SNRp > P_cut).any(): + vmin, vmax = 0.0, np.max([pol_err[SNRp > P_cut].max(), 1.0]) * 100.0 im = ax.imshow(pol_err * 100.0, vmin=vmin, vmax=vmax, aspect="equal", cmap="inferno_r", alpha=1.0) else: vmin, vmax = 0.0, 100.0 @@ -400,30 +414,39 @@ def polarization_map( else: im = ax.imshow(np.sqrt(stk_cov[0, 0]) * convert_flux, aspect="equal", cmap="inferno", alpha=1.0) fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=r"$\sigma_I$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") - elif display.lower() in ["snr", "snri"]: + elif display.lower() in ["snri"]: # Display I_stokes signal-to-noise map display = "snri" vmin, vmax = 0.0, np.max(SNRi[np.isfinite(SNRi)]) if vmax * 0.99 > SNRi_cut: im = ax.imshow(SNRi, vmin=vmin, vmax=vmax, aspect="equal", cmap="inferno", alpha=1.0) - levelsSNRi = np.linspace(SNRi_cut, vmax * 0.99, 5) + levelsSNRi = np.linspace(SNRi_cut, vmax * 0.99, 5).astype(int) print("SNRi contour levels : ", levelsSNRi) ax.contour(SNRi, levels=levelsSNRi, colors="grey", linewidths=0.5) else: im = ax.imshow(SNRi, aspect="equal", cmap="inferno", alpha=1.0) fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=r"$I_{Stokes}/\sigma_{I}$") - elif display.lower() in ["snrp"]: + elif display.lower() in ["snr", "snrp"]: # Display polarization degree signal-to-noise map display = "snrp" vmin, vmax = 0.0, np.max(SNRp[np.isfinite(SNRp)]) - if vmax * 0.99 > SNRp_cut: + if vmax * 0.99 > P_cut: im = ax.imshow(SNRp, vmin=vmin, vmax=vmax, aspect="equal", cmap="inferno", alpha=1.0) - levelsSNRp = np.linspace(SNRp_cut, vmax * 0.99, 5) + levelsSNRp = np.linspace(P_cut, vmax * 0.99, 5).astype(int) print("SNRp contour levels : ", levelsSNRp) ax.contour(SNRp, levels=levelsSNRp, colors="grey", linewidths=0.5) else: im = ax.imshow(SNRp, aspect="equal", cmap="inferno", alpha=1.0) fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=r"$P/\sigma_{P}$") + elif display.lower() in ["conf", "confp"]: + # Display polarization degree signal-to-noise map + display = "confp" + vmin, vmax = 0.0, 1.0 + im = ax.imshow(confP, vmin=vmin, vmax=vmax, aspect="equal", cmap="inferno", alpha=1.0) + levelsconfp = np.array([0.500, 0.900, 0.990, 0.999]) + print("confp contour levels : ", levelsconfp) + ax.contour(confP, levels=levelsconfp, colors="grey", linewidths=0.5) + fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=r"$Conf_{P}$") else: # Defaults to intensity map if mask.sum() > 0.0: @@ -463,7 +486,7 @@ def polarization_map( arrow_props={"ec": "k", "fc": "w", "alpha": 1, "lw": 1}, ) - if display.lower() in ["i", "s_i", "snri", "pf", "p", "pa", "s_p", "snrp"]: + if display.lower() in ["i", "s_i", "snri", "pf", "p", "pa", "s_p", "snrp", "confp"]: if step_vec == 0: poldata[np.isfinite(poldata)] = 1.0 / 2.0 step_vec = 1 @@ -536,7 +559,6 @@ def polarization_map( savename += ".pdf" fig.savefig(path_join(plots_folder, savename), bbox_inches="tight", dpi=150, facecolor="None") - plt.show() return fig, ax @@ -765,14 +787,14 @@ class align_maps(object): x = event.xdata y = event.ydata - self.cr_map.set(data=[x, y]) + self.cr_map.set(data=[[x], [y]]) self.fig_align.canvas.draw_idle() if (event.inaxes is not None) and (event.inaxes == self.other_ax): x = event.xdata y = event.ydata - self.cr_other.set(data=[x, y]) + self.cr_other.set(data=[[x], [y]]) self.fig_align.canvas.draw_idle() def reset_align(self, event): @@ -843,16 +865,23 @@ class overplot_radio(align_maps): Inherit from class align_maps in order to get the same WCS on both maps. """ - def overplot(self, levels=None, SNRp_cut=3.0, SNRi_cut=3.0, scale_vec=2, savename=None, **kwargs): + def overplot(self, levels=None, P_cut=0.99, SNRi_cut=1.0, scale_vec=2, savename=None, **kwargs): self.Stokes_UV = self.map self.wcs_UV = self.map_wcs # Get Data obj = self.Stokes_UV[0].header["targname"] stkI = self.Stokes_UV["I_STOKES"].data + stkQ = self.Stokes_UV["Q_STOKES"].data + stkU = self.Stokes_UV["U_STOKES"].data stk_cov = self.Stokes_UV["IQU_COV_MATRIX"].data pol = deepcopy(self.Stokes_UV["POL_DEG_DEBIASED"].data) pol_err = self.Stokes_UV["POL_DEG_ERR"].data pang = self.Stokes_UV["POL_ANG"].data + # Compute confidence level map + QN, UN, QN_ERR, UN_ERR = np.full((4, stkI.shape[0], stkI.shape[1]), np.nan) + for nflux, sflux in zip([QN, UN, QN_ERR, UN_ERR], [stkQ, stkU, np.sqrt(stk_cov[1, 1]), np.sqrt(stk_cov[2, 2])]): + nflux[stkI > 0.0] = sflux[stkI > 0.0] / stkI[stkI > 0.0] + confP = PCconf(QN, UN, QN_ERR, UN_ERR) other_data = self.other_data self.other_convert = 1.0 @@ -864,12 +893,17 @@ class overplot_radio(align_maps): self.map_convert = self.Stokes_UV[0].header["photflam"] # Compute SNR and apply cuts - pol[pol == 0.0] = np.nan - SNRp = pol / pol_err - SNRp[np.isnan(SNRp)] = 0.0 - pol[SNRp < SNRp_cut] = np.nan - SNRi = stkI / np.sqrt(stk_cov[0, 0]) - SNRi[np.isnan(SNRi)] = 0.0 + maskP = np.zeros(pol.shape, dtype=bool) + if P_cut >= 1.0: + SNRp = np.zeros(pol.shape) + SNRp[pol_err > 0.0] = pol[pol_err > 0.0] / pol_err[pol_err > 0.0] + maskP = SNRp > P_cut + else: + maskP = confP > P_cut + pol[np.logical_not(maskP)] = np.nan + + SNRi = np.zeros(stkI.shape) + SNRi[stk_cov[0, 0] > 0.0] = stkI[stk_cov[0, 0] > 0.0] / np.sqrt(stk_cov[0, 0][stk_cov[0, 0] > 0.0]) pol[SNRi < SNRi_cut] = np.nan plt.rcParams.update({"font.size": 16}) @@ -1025,7 +1059,7 @@ class overplot_radio(align_maps): (0, 0), (0, 1), arrowstyle="-", fc="w", ec="k", lw=2 ) labels.append("{0:s} contour".format(self.other_observer)) - handles.append(Rectangle((0, 0), 1, 1, fill=False, lw=2, ec=other_cont.collections[0].get_edgecolor()[0])) + handles.append(Rectangle((0, 0), 1, 1, fill=False, lw=2, ec=other_cont.get_edgecolor()[0])) self.legend = self.ax_overplot.legend( handles=handles, labels=labels, bbox_to_anchor=(0.0, 1.02, 1.0, 0.102), loc="lower left", mode="expand", borderaxespad=0.0 ) @@ -1037,10 +1071,10 @@ class overplot_radio(align_maps): self.fig_overplot.canvas.draw() - def plot(self, levels=None, SNRp_cut=3.0, SNRi_cut=3.0, savename=None, **kwargs) -> None: + def plot(self, levels=None, P_cut=0.99, SNRi_cut=1.0, savename=None, **kwargs) -> None: while not self.aligned: self.align() - self.overplot(levels=levels, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, savename=savename, **kwargs) + self.overplot(levels=levels, P_cut=P_cut, SNRi_cut=SNRi_cut, savename=savename, **kwargs) plt.show(block=True) @@ -1050,16 +1084,23 @@ class overplot_chandra(align_maps): Inherit from class align_maps in order to get the same WCS on both maps. """ - def overplot(self, levels=None, SNRp_cut=3.0, SNRi_cut=3.0, scale_vec=2, zoom=1, savename=None, **kwargs): + def overplot(self, levels=None, P_cut=0.99, SNRi_cut=1.0, scale_vec=2, zoom=1, savename=None, **kwargs): self.Stokes_UV = self.map self.wcs_UV = self.map_wcs # Get Data obj = self.Stokes_UV[0].header["targname"] stkI = self.Stokes_UV["I_STOKES"].data + stkQ = self.Stokes_UV["Q_STOKES"].data + stkU = self.Stokes_UV["U_STOKES"].data stk_cov = self.Stokes_UV["IQU_COV_MATRIX"].data pol = deepcopy(self.Stokes_UV["POL_DEG_DEBIASED"].data) pol_err = self.Stokes_UV["POL_DEG_ERR"].data pang = self.Stokes_UV["POL_ANG"].data + # Compute confidence level map + QN, UN, QN_ERR, UN_ERR = np.full((4, stkI.shape[0], stkI.shape[1]), np.nan) + for nflux, sflux in zip([QN, UN, QN_ERR, UN_ERR], [stkQ, stkU, np.sqrt(stk_cov[1, 1]), np.sqrt(stk_cov[2, 2])]): + nflux[stkI > 0.0] = sflux[stkI > 0.0] / stkI[stkI > 0.0] + confP = PCconf(QN, UN, QN_ERR, UN_ERR) other_data = deepcopy(self.other_data) other_wcs = self.other_wcs.deepcopy() @@ -1070,12 +1111,17 @@ class overplot_chandra(align_maps): self.other_unit = "counts" # Compute SNR and apply cuts - pol[pol == 0.0] = np.nan - SNRp = pol / pol_err - SNRp[np.isnan(SNRp)] = 0.0 - pol[SNRp < SNRp_cut] = np.nan - SNRi = stkI / np.sqrt(stk_cov[0, 0]) - SNRi[np.isnan(SNRi)] = 0.0 + maskP = np.zeros(pol.shape, dtype=bool) + if P_cut >= 1.0: + SNRp = np.zeros(pol.shape) + SNRp[pol_err > 0.0] = pol[pol_err > 0.0] / pol_err[pol_err > 0.0] + maskP = SNRp > P_cut + else: + maskP = confP > P_cut + pol[np.logical_not(maskP)] = np.nan + + SNRi = np.zeros(stkI.shape) + SNRi[stk_cov[0, 0] > 0.0] = stkI[stk_cov[0, 0] > 0.0] / np.sqrt(stk_cov[0, 0][stk_cov[0, 0] > 0.0]) pol[SNRi < SNRi_cut] = np.nan plt.rcParams.update({"font.size": 16}) @@ -1224,7 +1270,7 @@ class overplot_chandra(align_maps): (0, 0), (0, 1), arrowstyle="-", fc="w", ec="k", lw=2 ) labels.append("{0:s} contour in counts".format(self.other_observer)) - handles.append(Rectangle((0, 0), 1, 1, fill=False, lw=2, ec=other_cont.collections[0].get_edgecolor()[0])) + handles.append(Rectangle((0, 0), 1, 1, fill=False, lw=2, ec=other_cont.get_edgecolor()[0])) self.legend = self.ax_overplot.legend( handles=handles, labels=labels, bbox_to_anchor=(0.0, 1.02, 1.0, 0.102), loc="lower left", mode="expand", borderaxespad=0.0 ) @@ -1236,10 +1282,10 @@ class overplot_chandra(align_maps): self.fig_overplot.canvas.draw() - def plot(self, levels=None, SNRp_cut=3.0, SNRi_cut=3.0, zoom=1, savename=None, **kwargs) -> None: + def plot(self, levels=None, P_cut=0.99, SNRi_cut=1.0, zoom=1, savename=None, **kwargs) -> None: while not self.aligned: self.align() - self.overplot(levels=levels, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, zoom=zoom, savename=savename, **kwargs) + self.overplot(levels=levels, P_cut=P_cut, SNRi_cut=SNRi_cut, zoom=zoom, savename=savename, **kwargs) plt.show(block=True) @@ -1249,26 +1295,38 @@ class overplot_pol(align_maps): Inherit from class align_maps in order to get the same WCS on both maps. """ - def overplot(self, levels=None, SNRp_cut=3.0, SNRi_cut=3.0, scale_vec=2.0, savename=None, **kwargs): + def overplot(self, levels=None, P_cut=0.99, SNRi_cut=1.0, scale_vec=2.0, savename=None, **kwargs): self.Stokes_UV = self.map self.wcs_UV = self.map_wcs # Get Data obj = self.Stokes_UV[0].header["targname"] stkI = self.Stokes_UV["I_STOKES"].data + stkQ = self.Stokes_UV["Q_STOKES"].data + stkU = self.Stokes_UV["U_STOKES"].data stk_cov = self.Stokes_UV["IQU_COV_MATRIX"].data pol = deepcopy(self.Stokes_UV["POL_DEG_DEBIASED"].data) pol_err = self.Stokes_UV["POL_DEG_ERR"].data pang = self.Stokes_UV["POL_ANG"].data + # Compute confidence level map + QN, UN, QN_ERR, UN_ERR = np.full((4, stkI.shape[0], stkI.shape[1]), np.nan) + for nflux, sflux in zip([QN, UN, QN_ERR, UN_ERR], [stkQ, stkU, np.sqrt(stk_cov[1, 1]), np.sqrt(stk_cov[2, 2])]): + nflux[stkI > 0.0] = sflux[stkI > 0.0] / stkI[stkI > 0.0] + confP = PCconf(QN, UN, QN_ERR, UN_ERR) other_data = self.other_data # Compute SNR and apply cuts - pol[pol == 0.0] = np.nan - SNRp = pol / pol_err - SNRp[np.isnan(SNRp)] = 0.0 - pol[SNRp < SNRp_cut] = np.nan - SNRi = stkI / np.sqrt(stk_cov[0, 0]) - SNRi[np.isnan(SNRi)] = 0.0 + maskP = np.zeros(pol.shape, dtype=bool) + if P_cut >= 1.0: + SNRp = np.zeros(pol.shape) + SNRp[pol_err > 0.0] = pol[pol_err > 0.0] / pol_err[pol_err > 0.0] + maskP = SNRp > P_cut + else: + maskP = confP > P_cut + pol[np.logical_not(maskP)] = np.nan + + SNRi = np.zeros(stkI.shape) + SNRi[stk_cov[0, 0] > 0.0] = stkI[stk_cov[0, 0] > 0.0] / np.sqrt(stk_cov[0, 0][stk_cov[0, 0] > 0.0]) pol[SNRi < SNRi_cut] = np.nan plt.rcParams.update({"font.size": 16}) @@ -1437,7 +1495,7 @@ class overplot_pol(align_maps): (0, 0), (0, 1), arrowstyle="-", fc="w", ec="k", lw=2 ) labels.append("{0:s} Stokes I contour".format(self.map_observer)) - handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=cont_stkI.collections[0].get_edgecolor()[0])) + handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=cont_stkI.get_edgecolor()[0])) self.legend = self.ax_overplot.legend( handles=handles, labels=labels, bbox_to_anchor=(0.0, 1.02, 1.0, 0.102), loc="lower left", mode="expand", borderaxespad=0.0 ) @@ -1449,10 +1507,10 @@ class overplot_pol(align_maps): self.fig_overplot.canvas.draw() - def plot(self, levels=None, SNRp_cut=3.0, SNRi_cut=3.0, scale_vec=2.0, savename=None, **kwargs) -> None: + def plot(self, levels=None, P_cut=0.99, SNRi_cut=1.0, scale_vec=2.0, savename=None, **kwargs) -> None: while not self.aligned: self.align() - self.overplot(levels=levels, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, scale_vec=scale_vec, savename=savename, **kwargs) + self.overplot(levels=levels, P_cut=P_cut, SNRi_cut=SNRi_cut, scale_vec=scale_vec, savename=savename, **kwargs) plt.show(block=True) def add_vector(self, position="center", pol_deg=1.0, pol_ang=0.0, **kwargs): @@ -1462,7 +1520,12 @@ class overplot_pol(align_maps): position = self.other_wcs.world_to_pixel(position) u, v = pol_deg * np.cos(np.radians(pol_ang) + np.pi / 2.0), pol_deg * np.sin(np.radians(pol_ang) + np.pi / 2.0) - for key, value in [["scale", [["scale", self.scale_vec]]], ["width", [["width", 0.1]]], ["color", [["color", "k"]]], ["edgecolor", [["edgecolor", "w"]]]]: + for key, value in [ + ["scale", [["scale", self.scale_vec]]], + ["width", [["width", 0.1]]], + ["color", [["color", "k"]]], + ["edgecolor", [["edgecolor", "w"]]], + ]: try: _ = kwargs[key] except KeyError: @@ -1495,7 +1558,7 @@ class align_pol(object): self.kwargs = kwargs - def single_plot(self, curr_map, wcs, v_lim=None, ax_lim=None, SNRp_cut=3.0, SNRi_cut=3.0, savename=None, **kwargs): + def single_plot(self, curr_map, wcs, v_lim=None, ax_lim=None, P_cut=3.0, SNRi_cut=3.0, savename=None, **kwargs): # Get data stkI = curr_map["I_STOKES"].data stk_cov = curr_map["IQU_COV_MATRIX"].data @@ -1518,7 +1581,7 @@ class align_pol(object): SNRi = np.zeros(stkI.shape) SNRi[maskI] = stkI[maskI] / np.sqrt(stk_cov[0, 0][maskI]) - mask = (SNRp > SNRp_cut) * (SNRi > SNRi_cut) * (pol >= 0.0) + mask = (SNRp > P_cut) * (SNRi > SNRi_cut) * (pol >= 0.0) pol[mask] = np.nan # Plot the map @@ -1627,7 +1690,7 @@ class align_pol(object): self.wcs, self.wcs_other[i] = curr_align.align() self.aligned[i] = curr_align.aligned - def plot(self, SNRp_cut=3.0, SNRi_cut=3.0, savename=None, **kwargs): + def plot(self, P_cut=3.0, SNRi_cut=3.0, savename=None, **kwargs): while not self.aligned.all(): self.align() eps = 1e-35 @@ -1673,13 +1736,13 @@ class align_pol(object): ) v_lim = np.array([vmin, vmax]) - fig, ax = self.single_plot(self.ref_map, self.wcs, v_lim=v_lim, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, savename=savename + "_0", **kwargs) + fig, ax = self.single_plot(self.ref_map, self.wcs, v_lim=v_lim, P_cut=P_cut, SNRi_cut=SNRi_cut, savename=savename + "_0", **kwargs) x_lim, y_lim = ax.get_xlim(), ax.get_ylim() ax_lim = np.array([self.wcs.pixel_to_world(x_lim[i], y_lim[i]) for i in range(len(x_lim))]) for i, curr_map in enumerate(self.other_maps): self.single_plot( - curr_map, self.wcs_other[i], v_lim=v_lim, ax_lim=ax_lim, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, savename=savename + "_" + str(i + 1), **kwargs + curr_map, self.wcs_other[i], v_lim=v_lim, ax_lim=ax_lim, P_cut=P_cut, SNRi_cut=SNRi_cut, savename=savename + "_" + str(i + 1), **kwargs ) @@ -1726,6 +1789,7 @@ class crop_map(object): self.mask_alpha = 0.75 self.rect_selector = RectangleSelector(self.ax, self.onselect_crop, button=[1]) self.embedded = True + self.ax.set(xlabel="Right Ascension (J2000)", ylabel="Declination (J2000)") self.display(self.data, self.wcs, self.map_convert, **self.kwargs) self.extent = np.array([0.0, self.data.shape[0], 0.0, self.data.shape[1]]) @@ -2024,7 +2088,7 @@ class image_lasso_selector(object): self.mask = np.zeros(self.img.shape[:2], dtype=bool) self.mask[self.indices] = True if hasattr(self, "cont"): - for coll in self.cont.collections: + for coll in self.cont: coll.remove() self.cont = self.ax.contour(self.mask.astype(float), levels=[0.5], colors="white", linewidths=1) if not self.embedded: @@ -2137,7 +2201,7 @@ class slit(object): for p in self.pix: self.mask[tuple(p)] = (np.abs(np.dot(rot2D(-self.angle), p - self.rect.get_center()[::-1])) < (self.height / 2.0, self.width / 2.0)).all() if hasattr(self, "cont"): - for coll in self.cont.collections: + for coll in self.cont: try: coll.remove() except AttributeError: @@ -2240,7 +2304,7 @@ class aperture(object): x0, y0 = self.circ.center self.mask = np.sqrt((xx - x0) ** 2 + (yy - y0) ** 2) < self.radius if hasattr(self, "cont"): - for coll in self.cont.collections: + for coll in self.cont: try: coll.remove() except AttributeError: @@ -2258,21 +2322,22 @@ class pol_map(object): Class to interactively study polarization maps. """ - def __init__(self, Stokes, SNRp_cut=3.0, SNRi_cut=3.0, step_vec=1, scale_vec=3.0, flux_lim=None, selection=None, pa_err=False): + def __init__(self, Stokes, P_cut=0.99, SNRi_cut=1.0, step_vec=1, scale_vec=3.0, flux_lim=None, selection=None, pa_err=False): if isinstance(Stokes, str): Stokes = fits.open(Stokes) self.Stokes = deepcopy(Stokes) - self.SNRp_cut = SNRp_cut + self.P_cut = P_cut self.SNRi_cut = SNRi_cut self.flux_lim = flux_lim self.SNRi = deepcopy(self.SNRi_cut) - self.SNRp = deepcopy(self.SNRp_cut) + self.SNRp = deepcopy(self.P_cut) self.region = None self.data = None self.display_selection = selection self.step_vec = step_vec self.scale_vec = scale_vec self.pa_err = pa_err + self.conf = PCconf(self.Q / self.I, self.U / self.I, np.sqrt(self.IQU_cov[1, 1]) / self.I, np.sqrt(self.IQU_cov[2, 2]) / self.I) # Get data self.targ = self.Stokes[0].header["targname"] @@ -2292,18 +2357,22 @@ class pol_map(object): # Display integrated values in ROI self.pol_int() - # Set axes for sliders (SNRp_cut, SNRi_cut) - ax_I_cut = self.fig.add_axes([0.120, 0.080, 0.230, 0.01]) - ax_P_cut = self.fig.add_axes([0.120, 0.055, 0.230, 0.01]) - ax_vec_sc = self.fig.add_axes([0.240, 0.030, 0.110, 0.01]) - ax_snr_reset = self.fig.add_axes([0.080, 0.020, 0.05, 0.02]) + # Set axes for sliders (P_cut, SNRi_cut) + ax_I_cut = self.fig.add_axes([0.120, 0.080, 0.220, 0.01]) + self.ax_P_cut = self.fig.add_axes([0.120, 0.055, 0.220, 0.01]) + ax_vec_sc = self.fig.add_axes([0.260, 0.030, 0.090, 0.01]) + ax_snr_reset = self.fig.add_axes([0.060, 0.020, 0.05, 0.02]) + ax_snr_conf = self.fig.add_axes([0.115, 0.020, 0.05, 0.02]) SNRi_max = np.max(self.I[self.IQU_cov[0, 0] > 0.0] / np.sqrt(self.IQU_cov[0, 0][self.IQU_cov[0, 0] > 0.0])) SNRp_max = np.max(self.P[self.s_P > 0.0] / self.s_P[self.s_P > 0.0]) s_I_cut = Slider(ax_I_cut, r"$SNR^{I}_{cut}$", 1.0, int(SNRi_max * 0.95), valstep=1, valinit=self.SNRi_cut) - s_P_cut = Slider(ax_P_cut, r"$SNR^{P}_{cut}$", 1.0, int(SNRp_max * 0.95), valstep=1, valinit=self.SNRp_cut) - s_vec_sc = Slider(ax_vec_sc, r"Vectors scale", 0.0, 10.0, valstep=1, valinit=self.scale_vec) + self.s_P_cut = Slider(self.ax_P_cut, r"$Conf^{P}_{cut}$", 0.50, 1.0, valstep=[0.500, 0.900, 0.990, 0.999], valinit=self.P_cut if P_cut <= 1.0 else 0.99) + s_vec_sc = Slider(ax_vec_sc, r"Vec scale", 0.0, 10.0, valstep=1, valinit=self.scale_vec) b_snr_reset = Button(ax_snr_reset, "Reset") b_snr_reset.label.set_fontsize(8) + self.snr_conf = 1 + b_snr_conf = Button(ax_snr_conf, "Conf") + b_snr_conf.label.set_fontsize(8) def update_snri(val): self.SNRi = val @@ -2325,13 +2394,34 @@ class pol_map(object): def reset_snr(event): s_I_cut.reset() - s_P_cut.reset() + self.s_P_cut.reset() s_vec_sc.reset() + def toggle_snr_conf(event=None): + self.ax_P_cut.remove() + self.ax_P_cut = self.fig.add_axes([0.120, 0.055, 0.220, 0.01]) + if self.snr_conf: + self.snr_conf = 0 + b_snr_conf.label.set_text("Conf") + self.s_P_cut = Slider(self.ax_P_cut, r"$SNR^{P}_{cut}$", 1.0, int(SNRp_max * 0.95), valstep=1, valinit=self.P_cut if P_cut >= 1.0 else 3) + else: + self.snr_conf = 1 + b_snr_conf.label.set_text("SNR") + self.s_P_cut = Slider( + self.ax_P_cut, r"$Conf^{P}_{cut}$", 0.50, 1.0, valstep=[0.500, 0.900, 0.990, 0.999], valinit=self.P_cut if P_cut <= 1.0 else 0.99 + ) + self.s_P_cut.on_changed(update_snrp) + update_snrp(self.s_P_cut.val) + self.fig.canvas.draw_idle() + s_I_cut.on_changed(update_snri) - s_P_cut.on_changed(update_snrp) + self.s_P_cut.on_changed(update_snrp) s_vec_sc.on_changed(update_vecsc) b_snr_reset.on_clicked(reset_snr) + b_snr_conf.on_clicked(toggle_snr_conf) + + if self.P_cut >= 1.0: + toggle_snr_conf() # Set axe for ROI selection ax_select = self.fig.add_axes([0.375, 0.070, 0.05, 0.02]) @@ -2349,7 +2439,7 @@ class pol_map(object): self.selected = False self.region = deepcopy(self.select_instance.mask.astype(bool)) self.select_instance.displayed.remove() - for coll in self.select_instance.cont.collections: + for coll in self.select_instance.cont: coll.remove() self.select_instance.lasso.set_active(False) self.set_data_mask(deepcopy(self.region)) @@ -2393,7 +2483,7 @@ class pol_map(object): self.select_instance.update_mask() self.region = deepcopy(self.select_instance.mask.astype(bool)) self.select_instance.displayed.remove() - for coll in self.select_instance.cont.collections[:]: + for coll in self.select_instance.cont: coll.remove() self.select_instance.circ.set_visible(False) self.set_data_mask(deepcopy(self.region)) @@ -2451,7 +2541,7 @@ class pol_map(object): self.select_instance.update_mask() self.region = deepcopy(self.select_instance.mask.astype(bool)) self.select_instance.displayed.remove() - for coll in self.select_instance.cont.collections[:]: + for coll in self.select_instance.cont: coll.remove() self.select_instance.rect.set_visible(False) self.set_data_mask(deepcopy(self.region)) @@ -2768,8 +2858,11 @@ class pol_map(object): def cut(self): s_I = np.sqrt(self.IQU_cov[0, 0]) SNRp_mask, SNRi_mask = np.zeros(self.P.shape).astype(bool), np.zeros(self.I.shape).astype(bool) - SNRp_mask[self.s_P > 0.0] = self.P[self.s_P > 0.0] / self.s_P[self.s_P > 0.0] > self.SNRp SNRi_mask[s_I > 0.0] = self.I[s_I > 0.0] / s_I[s_I > 0.0] > self.SNRi + if self.SNRp >= 1.0: + SNRp_mask[self.s_P > 0.0] = self.P[self.s_P > 0.0] / self.s_P[self.s_P > 0.0] > self.SNRp + else: + SNRp_mask = self.conf > self.SNRp return np.logical_and(SNRi_mask, SNRp_mask) def ax_cosmetics(self, ax=None): @@ -3171,7 +3264,7 @@ class pol_map(object): ) if hasattr(self, "cont"): - for coll in self.cont.collections: + for coll in self.cont: try: coll.remove() except AttributeError: @@ -3242,15 +3335,13 @@ if __name__ == "__main__": parser = argparse.ArgumentParser(description="Interactively plot the pipeline products") 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( - "-p", "--snrp", metavar="snrp_cut", required=False, help="The cut in signal-to-noise for the polarization degree", type=float, default=3.0 - ) + parser.add_argument("-p", "--pcut", metavar="p_cut", required=False, help="The cut in signal-to-noise for the polarization degree", type=float, default=3.0) parser.add_argument("-i", "--snri", metavar="snri_cut", required=False, help="The cut in signal-to-noise for the intensity", type=float, default=3.0) parser.add_argument( "-st", "--step-vec", metavar="step_vec", required=False, help="Quantity of vectors to be shown, 1 is all, 2 is every other, etc.", type=int, default=1 ) parser.add_argument( - "-sc", "--scale-vec", metavar="scale_vec", required=False, help="Size of the 100% polarization vector in pixel units", type=float, default=3.0 + "-sc", "--scale-vec", metavar="scale_vec", required=False, help="Size of the 100%% polarization vector in pixel units", type=float, default=3.0 ) parser.add_argument("-pa", "--pang-err", action="store_true", required=False, help="Whether the polarization angle uncertainties should be displayed") parser.add_argument("-l", "--lim", metavar="flux_lim", nargs=2, required=False, help="Limits for the intensity map", type=float, default=None) @@ -3265,7 +3356,7 @@ if __name__ == "__main__": polarization_map( Stokes_UV, Stokes_UV["DATA_MASK"].data.astype(bool), - SNRp_cut=args.snrp, + P_cut=args.p_cut, SNRi_cut=args.snri, flux_lim=args.lim, step_vec=args.step_vec, @@ -3276,7 +3367,7 @@ if __name__ == "__main__": polarization_map( Stokes_UV, Stokes_UV["DATA_MASK"].data.astype(bool), - SNRp_cut=args.snrp, + P_cut=args.p_cut, SNRi_cut=args.snri, flux_lim=args.lim, step_vec=args.step_vec, @@ -3288,7 +3379,7 @@ if __name__ == "__main__": polarization_map( Stokes_UV, Stokes_UV["DATA_MASK"].data.astype(bool), - SNRp_cut=args.snrp, + P_cut=args.p_cut, SNRi_cut=args.snri, flux_lim=args.lim, step_vec=args.step_vec, @@ -3300,7 +3391,7 @@ if __name__ == "__main__": polarization_map( Stokes_UV, Stokes_UV["DATA_MASK"].data.astype(bool), - SNRp_cut=args.snrp, + P_cut=args.p_cut, SNRi_cut=args.snri, flux_lim=args.lim, step_vec=args.step_vec, @@ -3312,7 +3403,7 @@ if __name__ == "__main__": polarization_map( Stokes_UV, Stokes_UV["DATA_MASK"].data.astype(bool), - SNRp_cut=args.snrp, + P_cut=args.p_cut, SNRi_cut=args.snri, flux_lim=args.lim, step_vec=args.step_vec, @@ -3324,7 +3415,7 @@ if __name__ == "__main__": polarization_map( Stokes_UV, Stokes_UV["DATA_MASK"].data.astype(bool), - SNRp_cut=args.snrp, + P_cut=args.p_cut, SNRi_cut=args.snri, flux_lim=args.lim, step_vec=args.step_vec, @@ -3336,7 +3427,7 @@ if __name__ == "__main__": polarization_map( Stokes_UV, Stokes_UV["DATA_MASK"].data.astype(bool), - SNRp_cut=args.snrp, + P_cut=args.p_cut, SNRi_cut=args.snri, flux_lim=args.lim, step_vec=args.step_vec, @@ -3348,7 +3439,7 @@ if __name__ == "__main__": polarization_map( Stokes_UV, Stokes_UV["DATA_MASK"].data.astype(bool), - SNRp_cut=args.snrp, + P_cut=args.p_cut, SNRi_cut=args.snri, flux_lim=args.lim, step_vec=args.step_vec, @@ -3360,7 +3451,7 @@ if __name__ == "__main__": polarization_map( Stokes_UV, Stokes_UV["DATA_MASK"].data.astype(bool), - SNRp_cut=args.snrp, + P_cut=args.p_cut if args.p_cut >= 1.0 else 3.0, SNRi_cut=args.snri, flux_lim=args.lim, step_vec=args.step_vec, @@ -3369,12 +3460,21 @@ if __name__ == "__main__": plots_folder=args.static_pdf, display="SNRp", ) - else: - pol_map( - Stokes_UV, SNRp_cut=args.snrp, SNRi_cut=args.snri, step_vec=args.step_vec, scale_vec=args.scale_vec, flux_lim=args.lim, pa_err=args.pang_err + polarization_map( + Stokes_UV, + Stokes_UV["DATA_MASK"].data.astype(bool), + P_cut=args.p_cut if args.p_cut < 1.0 else 0.99, + SNRi_cut=args.snri, + flux_lim=args.lim, + step_vec=args.step_vec, + scale_vec=args.scale_vec, + savename="_".join([Stokes_UV[0].header["FILENAME"], "confP"]), + plots_folder=args.static_pdf, + display="confp", ) - + else: + pol_map(Stokes_UV, P_cut=args.p_cut, SNRi_cut=args.snri, step_vec=args.step_vec, scale_vec=args.scale_vec, flux_lim=args.lim, pa_err=args.pang_err) else: print( - "python3 plots.py -f -p -i -st -sc -l -pa --pdf " + "python3 plots.py -f -p -i -st -sc -l -pa --pdf " ) diff --git a/package/Combine.py b/package/src/Combine.py similarity index 98% rename from package/Combine.py rename to package/src/Combine.py index b3871f1..87b9632 100755 --- a/package/Combine.py +++ b/package/src/Combine.py @@ -1,6 +1,9 @@ #!/usr/bin/python # -*- coding:utf-8 -*- -# Project libraries +from pathlib import Path +from sys import path as syspath + +syspath.append(str(Path(__file__).parent.parent)) import numpy as np diff --git a/package/src/__init__.py b/package/src/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/package/src/analysis.py b/package/src/analysis.py deleted file mode 100755 index 815eaa3..0000000 --- a/package/src/analysis.py +++ /dev/null @@ -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 -p -i -l ") - 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 -p -i -l ") diff --git a/package/src/comparison_Kishimoto.py b/package/src/comparison_Kishimoto.py deleted file mode 100755 index 75b7073..0000000 --- a/package/src/comparison_Kishimoto.py +++ /dev/null @@ -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() diff --git a/package/src/emission_center.py b/package/src/emission_center.py new file mode 100755 index 0000000..829fb54 --- /dev/null +++ b/package/src/emission_center.py @@ -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) diff --git a/package/src/get_cdelt.py b/package/src/get_cdelt.py index b7054c6..d6ea688 100755 --- a/package/src/get_cdelt.py +++ b/package/src/get_cdelt.py @@ -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): diff --git a/package/overplot_IC5063.py b/package/src/overplot_IC5063.py similarity index 66% rename from package/overplot_IC5063.py rename to package/src/overplot_IC5063.py index 6a4fb1e..f9ac6c3 100755 --- a/package/overplot_IC5063.py +++ b/package/src/overplot_IC5063.py @@ -1,4 +1,10 @@ #!/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 @@ -18,33 +24,33 @@ 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, SNRp_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/18GHz_overplot.pdf", vec_scale=None) +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, SNRp_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/24GHz_overplot.pdf", vec_scale=None) +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, SNRp_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/103GHz_overplot.pdf", vec_scale=None) +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, SNRp_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/229GHz_overplot.pdf", vec_scale=None) +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, SNRp_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/357GHz_overplot.pdf", vec_scale=None) +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(SNRp_cut=3.0, SNRi_cut=80.0, savename='./plots/IC5063/S2_overplot.pdf', norm=LogNorm(vmin=5e-20,vmax=5e-18)) +# 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( - SNRp_cut=2.0, + P_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/IR_overplot.pdf", - vec_scale=None, + 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", ) diff --git a/package/overplot_MRK463E.py b/package/src/overplot_MRK463E.py similarity index 66% rename from package/overplot_MRK463E.py rename to package/src/overplot_MRK463E.py index 5c3411d..a571a7a 100755 --- a/package/overplot_MRK463E.py +++ b/package/src/overplot_MRK463E.py @@ -1,20 +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.10arcsec.fits") +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, SNRp_cut=3.0, SNRi_cut=3.0, vec_scale=5, zoom=1, savename="./plots/MRK463E/Chandra_overplot.pdf") +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, SNRp_cut=3.0, SNRi_cut=3.0, vec_scale=5, norm=LogNorm(8.5e-18, 2.5e-15), savename="./plots/MRK463E/IR_overplot.pdf") +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") diff --git a/package/test_center.py b/package/test_center.py deleted file mode 100644 index bedfda8..0000000 --- a/package/test_center.py +++ /dev/null @@ -1,109 +0,0 @@ -import matplotlib.pyplot as plt -import numpy as np -from astropy.io.fits import open as fits_open -from astropy.wcs import WCS -from lib.utils import CenterConf, PCconf -from matplotlib.colors import LogNorm -from matplotlib.patches import Rectangle - -levelssnr = np.array([3.0, 4.0]) -levelsconf = np.array([0.99]) - -NGC1068 = fits_open("./data/NGC1068/5144/NGC1068_FOC_b0.05arcsec_c0.07arcsec.fits") -NGC1068conf = PCconf( - NGC1068["Q_STOKES"].data / NGC1068["I_STOKES"].data, - NGC1068["U_STOKES"].data / NGC1068["I_STOKES"].data, - np.sqrt(NGC1068["IQU_COV_MATRIX"].data[1, 1]) / NGC1068["I_STOKES"].data, - np.sqrt(NGC1068["IQU_COV_MATRIX"].data[2, 2]) / NGC1068["I_STOKES"].data, -) -NGC1068mask = NGC1068["DATA_MASK"].data.astype(bool) -NGC1068snr = np.full(NGC1068mask.shape, np.nan) -NGC1068snr[NGC1068["POL_DEG_ERR"].data > 0.0] = ( - NGC1068["POL_DEG_DEBIASED"].data[NGC1068["POL_DEG_ERR"].data > 0.0] / NGC1068["POL_DEG_ERR"].data[NGC1068["POL_DEG_ERR"].data > 0.0] -) - -NGC1068centconf, NGC1068center = CenterConf(NGC1068conf > 0.99, NGC1068["POL_ANG"].data, NGC1068["POL_ANG_ERR"].data) -NGC1068pos = WCS(NGC1068[0].header).pixel_to_world(*NGC1068center) - -figngc, axngc = plt.subplots(1, 2, layout="tight", figsize=(18, 9), subplot_kw=dict(projection=WCS(NGC1068[0].header)), sharex=True, sharey=True) - -axngc[0].set(xlabel="RA", ylabel="DEC", title="NGC1069 intensity map with SNR and confidence contours") -vmin, vmax = ( - 0.5 * np.median(NGC1068["I_STOKES"].data[NGC1068mask]) * NGC1068[0].header["PHOTFLAM"], - np.max(NGC1068["I_STOKES"].data[NGC1068mask]) * NGC1068[0].header["PHOTFLAM"], -) -imngc = axngc[0].imshow(NGC1068["I_STOKES"].data * NGC1068["I_STOKES"].header["PHOTFLAM"], norm=LogNorm(vmin, vmax), cmap="inferno") -ngcsnrcont = axngc[0].contour(NGC1068snr, levelssnr, colors="b") -ngcconfcont = axngc[0].contour(NGC1068conf, levelsconf, colors="r") -ngcconfcenter = axngc[0].plot(*NGC1068center, marker="+",color="gray", label="Best confidence for center: {0}".format(NGC1068pos.to_string('hmsdms'))) -ngcconfcentcont = axngc[0].contour(NGC1068centconf, [0.01], colors="gray") -handles, labels = axngc[0].get_legend_handles_labels() -labels.append("SNR contours") -handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=ngcsnrcont.collections[0].get_edgecolor()[0])) -labels.append("CONF99 contour") -handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=ngcconfcont.collections[0].get_edgecolor()[0])) -labels.append("Center CONF99 contour") -handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=ngcconfcentcont.collections[0].get_edgecolor()[0])) -axngc[0].legend(handles=handles, labels=labels) - -axngc[1].set(xlabel="RA", ylabel="DEC", title="Location of the nucleus confidence map") -ngccent = axngc[1].imshow(NGC1068centconf, vmin=0.0, cmap="inferno") -ngccentcont = axngc[1].contour(NGC1068centconf, [0.01], colors="gray") -ngccentcenter = axngc[1].plot(*NGC1068center, marker="+",color="gray", label="Best confidence for center: {0}".format(NGC1068pos.to_string('hmsdms'))) -handles, labels = axngc[1].get_legend_handles_labels() -labels.append("CONF99 contour") -handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=ngccentcont.collections[0].get_edgecolor()[0])) -axngc[1].legend(handles=handles, labels=labels) - -figngc.savefig("NGC1068_center.pdf", dpi=150, facecolor="None") - -################################################################################################### - -MRK463E = fits_open("./data/MRK463E/5960/MRK463E_FOC_b0.05arcsec_c0.07arcsec.fits") -MRK463Econf = PCconf( - MRK463E["Q_STOKES"].data / MRK463E["I_STOKES"].data, - MRK463E["U_STOKES"].data / MRK463E["I_STOKES"].data, - np.sqrt(MRK463E["IQU_COV_MATRIX"].data[1, 1]) / MRK463E["I_STOKES"].data, - np.sqrt(MRK463E["IQU_COV_MATRIX"].data[2, 2]) / MRK463E["I_STOKES"].data, -) -MRK463Emask = MRK463E["DATA_MASK"].data.astype(bool) -MRK463Esnr = np.full(MRK463Emask.shape, np.nan) -MRK463Esnr[MRK463E["POL_DEG_ERR"].data > 0.0] = ( - MRK463E["POL_DEG_DEBIASED"].data[MRK463E["POL_DEG_ERR"].data > 0.0] / MRK463E["POL_DEG_ERR"].data[MRK463E["POL_DEG_ERR"].data > 0.0] -) - -MRK463Ecentconf, MRK463Ecenter = CenterConf(MRK463Econf > 0.99, MRK463E["POL_ANG"].data, MRK463E["POL_ANG_ERR"].data) -MRK463Epos = WCS(MRK463E[0].header).pixel_to_world(*MRK463Ecenter) - -figmrk, axmrk = plt.subplots(1, 2, layout="tight", figsize=(18, 9), subplot_kw=dict(projection=WCS(MRK463E[0].header)), sharex=True, sharey=True) - -axmrk[0].set(xlabel="RA", ylabel="DEC", title="NGC1069 intensity map with SNR and confidence contours") -vmin, vmax = ( - 0.5 * np.median(MRK463E["I_STOKES"].data[MRK463Emask]) * MRK463E[0].header["PHOTFLAM"], - np.max(MRK463E["I_STOKES"].data[MRK463Emask]) * MRK463E[0].header["PHOTFLAM"], -) -immrk = axmrk[0].imshow(MRK463E["I_STOKES"].data * MRK463E["I_STOKES"].header["PHOTFLAM"], norm=LogNorm(vmin, vmax), cmap="inferno") -mrksnrcont = axmrk[0].contour(MRK463Esnr, levelssnr, colors="b") -mrkconfcont = axmrk[0].contour(MRK463Econf, levelsconf, colors="r") -mrkconfcenter = axmrk[0].plot(*MRK463Ecenter, marker="+",color="gray", label="Best confidence for center: {0}".format(MRK463Epos.to_string('hmsdms'))) -mrkconfcentcont = axmrk[0].contour(MRK463Ecentconf, [0.01], colors="gray") -handles, labels = axmrk[0].get_legend_handles_labels() -labels.append("SNR contours") -handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=mrksnrcont.collections[0].get_edgecolor()[0])) -labels.append("CONF99 contour") -handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=mrkconfcont.collections[0].get_edgecolor()[0])) -labels.append("Center CONF99 contour") -handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=mrkconfcentcont.collections[0].get_edgecolor()[0])) -axmrk[0].legend(handles=handles, labels=labels) - -axmrk[1].set(xlabel="RA", ylabel="DEC", title="Location of the nucleus confidence map") -mrkcent = axmrk[1].imshow(MRK463Ecentconf, vmin=0.0, cmap="inferno") -mrkcentcont = axmrk[1].contour(MRK463Ecentconf, [0.01], colors="gray") -mrkcentcenter = axmrk[1].plot(*MRK463Ecenter, marker="+",color="gray", label="Best confidence for center: {0}".format(MRK463Epos.to_string('hmsdms'))) -handles, labels = axmrk[1].get_legend_handles_labels() -labels.append("CONF99 contour") -handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=mrkcentcont.collections[0].get_edgecolor()[0])) -axmrk[1].legend(handles=handles, labels=labels) - -figmrk.savefig("MRK463E_center.pdf", dpi=150, facecolor="None") -plt.show()