From 85cacfdb51693cb4b151c0defd95ae4aefb7caa8 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Mon, 18 Sep 2023 16:54:39 +0200 Subject: [PATCH] add the possibility to fix flux limits in output maps --- src/FOC_reduction.py | 21 +++++++++++---------- src/analysis.py | 13 ++++++++----- src/lib/plots.py | 45 +++++++++++++++++++++++++++++++------------- 3 files changed, 51 insertions(+), 28 deletions(-) diff --git a/src/FOC_reduction.py b/src/FOC_reduction.py index f27f145..30e8874 100755 --- a/src/FOC_reduction.py +++ b/src/FOC_reduction.py @@ -63,6 +63,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data"): # Polarization map output SNRp_cut = 3. #P measurments with SNR>3 SNRi_cut = 30. #I measurments with SNR>30, which implies an uncertainty in P of 4.7%. + flux_lim = [5e-19,5e-14] #lowest and highest flux displayed on plot, defaults to bkg and maximum in cut if None vec_scale = 2.0 step_vec = 1 #plot all vectors in the array. if step_vec = 2, then every other vector will be plotted # if step_vec = 0 then all vectors are displayed at full length @@ -182,19 +183,19 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data"): print("PA_bkg = {0:.1f} ± {1:.1f} °".format(PA_bkg[0,0],np.ceil(s_PA_bkg[0,0]*10.)/10.)) # Plot polarization map (Background is either total Flux, Polarization degree or Polarization degree error). if px_scale.lower() not in ['full','integrate'] and final_display: - proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype]), plots_folder=plots_folder) - proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"I"]), plots_folder=plots_folder, display='Intensity') - proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"P_flux"]), plots_folder=plots_folder, display='Pol_Flux') - proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"P"]), plots_folder=plots_folder, display='Pol_deg') - proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"PA"]), plots_folder=plots_folder, display='Pol_ang') - proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"I_err"]), plots_folder=plots_folder, display='I_err') - proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"P_err"]), plots_folder=plots_folder, display='Pol_deg_err') - proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"SNRi"]), plots_folder=plots_folder, display='SNRi') - proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"SNRp"]), plots_folder=plots_folder, display='SNRp') + proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype]), plots_folder=plots_folder) + proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"I"]), plots_folder=plots_folder, display='Intensity') + proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"P_flux"]), plots_folder=plots_folder, display='Pol_Flux') + proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"P"]), plots_folder=plots_folder, display='Pol_deg') + proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"PA"]), plots_folder=plots_folder, display='Pol_ang') + proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"I_err"]), plots_folder=plots_folder, display='I_err') + proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"P_err"]), plots_folder=plots_folder, display='Pol_deg_err') + proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"SNRi"]), plots_folder=plots_folder, display='SNRi') + proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"SNRp"]), plots_folder=plots_folder, display='SNRp') elif final_display: proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, savename="_".join([figname,figtype]), plots_folder=plots_folder, display='integrate') elif px_scale.lower() not in ['full', 'integrate']: - pol_map = proj_plots.pol_map(Stokes_test, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut) + pol_map = proj_plots.pol_map(Stokes_test, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim) return 0 diff --git a/src/analysis.py b/src/analysis.py index 949d83d..24775b7 100755 --- a/src/analysis.py +++ b/src/analysis.py @@ -3,11 +3,12 @@ from getopt import getopt, error as get_error from sys import argv arglist = argv[1:] -options = "hf:p:i:" -long_options = ["help","fits=","snrp=","snri="] +options = "hf:p:i:l:" +long_options = ["help","fits=","snrp=","snri=","lim="] fits_path = None SNRp_cut, SNRi_cut = 3, 30 +flux_lim = None out_txt = None try: @@ -15,13 +16,15 @@ try: for curr_arg, curr_val in arg: if curr_arg in ("-h", "--help"): - print("python3 analysis.py -f -p -i ") + 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)) @@ -30,7 +33,7 @@ if not fits_path is None: 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) + 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 ") + print("python3 analysis.py -f -p -i -l ") diff --git a/src/lib/plots.py b/src/lib/plots.py index a1d0bee..5aa607e 100755 --- a/src/lib/plots.py +++ b/src/lib/plots.py @@ -219,7 +219,7 @@ def plot_Stokes(Stokes, savename=None, plots_folder=""): def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_cut=30., - step_vec=1, vec_scale=2., savename=None, plots_folder="", display="default"): + flux_lim=None, step_vec=1, vec_scale=2., savename=None, plots_folder="", display="default"): """ Plots polarization map from Stokes HDUList. ---------- @@ -239,6 +239,10 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c Cut that should be applied to the signal-to-noise ratio on I. Any SNR < SNRi_cut won't be displayed. Defaults to 30. This value implies an uncertainty in P of 4.7% + flux_lim : float list, optional + Limits that should be applied to the flux colorbar. + Defaults to None, limits are computed on the background value and the + maximum value in the cut. step_vec : int, optional Number of steps between each displayed polarization vector. If step_vec = 2, every other vector will be displayed. @@ -320,10 +324,13 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c if display.lower() in ['intensity']: # If no display selected, show intensity map display='i' - if mask.sum() > 0.: - vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov.data[0,0][mask])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) + if flux_lim is None: + if mask.sum() > 0.: + vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov.data[0,0][mask])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) + else: + vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov.data[0,0][stkI.data > 0.])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) else: - vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov.data[0,0][stkI.data > 0.])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) + vmin, vmax = flux_lim im = ax.imshow(stkI.data*convert_flux, norm=LogNorm(vmin,vmax), aspect='equal', cmap='inferno', alpha=1.) cbar = plt.colorbar(im, cax=cbar_ax, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") levelsI = np.linspace(vmax*0.01, vmax*0.99, 10) @@ -334,10 +341,13 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c # Display polarisation flux display='pf' pf_mask = (stkI.data > 0.) * (pol.data > 0.) - if mask.sum() > 0.: - vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov.data[0,0][mask])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) + if flux_lim is None: + if mask.sum() > 0.: + vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov.data[0,0][mask])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) + else: + vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov.data[0,0][stkI.data > 0.])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) else: - vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov.data[0,0][stkI.data > 0.])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) + vmin, vmax = flux_lim im = ax.imshow(stkI.data*convert_flux*pol.data, norm=LogNorm(vmin,vmax), aspect='equal', cmap='inferno', alpha=1.) cbar = plt.colorbar(im, cax=cbar_ax, label=r"$F_{\lambda} \cdot P$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") levelsPf = np.linspace(vmax*0.01, vmax*0.99, 10) @@ -1353,13 +1363,14 @@ class pol_map(object): """ Class to interactively study polarization maps. """ - def __init__(self, Stokes, SNRp_cut=3., SNRi_cut=30., selection=None): + def __init__(self, Stokes, SNRp_cut=3., SNRi_cut=30., flux_lim=None, selection=None): if type(Stokes) == str: Stokes = fits.open(Stokes) self.Stokes = deepcopy(Stokes) self.SNRp_cut = SNRp_cut self.SNRi_cut = SNRi_cut + self.flux_lim = flux_lim self.SNRi = deepcopy(self.SNRi_cut) self.SNRp = deepcopy(self.SNRp_cut) self.region = None @@ -1782,18 +1793,26 @@ class pol_map(object): self.north_dir = AnchoredDirectionArrows(ax.transAxes, "E", "N", length=-0.08, fontsize=0.025, loc=1, aspect_ratio=-1, sep_y=0.01, sep_x=0.01, back_length=0., head_length=10., head_width=10., angle=-self.Stokes[0].header['orientat'], color='white', text_props={'ec': None, 'fc': 'w', 'alpha': 1, 'lw': 0.4}, arrow_props={'ec': None,'fc':'w','alpha': 1,'lw': 1}) ax.add_artist(self.north_dir) - def display(self, fig=None, ax=None): + def display(self, fig=None, ax=None, flux_lim=None): norm = None if self.display_selection is None: self.display_selection = "total_flux" + if flux_lim is None: + flux_lim = self.flux_lim if self.display_selection.lower() in ['total_flux']: self.data = self.I*self.convert_flux - vmin, vmax = 1./2.*np.median(self.data[self.data > 0.]), np.max(self.data[self.data > 0.]) + if flux_lim is None: + vmin, vmax = 1./2.*np.median(self.data[self.data > 0.]), np.max(self.data[self.data > 0.]) + else: + vmin, vmax = flux_lim norm = LogNorm(vmin, vmax) label = r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]" elif self.display_selection.lower() in ['pol_flux']: self.data = self.I*self.convert_flux*self.P - vmin, vmax = 1./2.*np.median(self.I[self.I > 0.]*self.convert_flux), np.max(self.I[self.I > 0.]*self.convert_flux) + if flux_lim is None: + vmin, vmax = 1./2.*np.median(self.I[self.I > 0.]*self.convert_flux), np.max(self.I[self.I > 0.]*self.convert_flux) + else: + vmin, vmax = flux_lim norm = LogNorm(vmin, vmax) label = r"$F_{\lambda} \cdot P$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]" elif self.display_selection.lower() in ['pol_deg']: @@ -1853,11 +1872,11 @@ class pol_map(object): ax = self.ax if hasattr(self, 'quiver'): self.quiver.remove() - self.quiver = ax.quiver(X, Y, XY_U, XY_V, units='xy', scale=1./self.vec_scale, scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.1, linewidth=0.5, color='white',edgecolor='black') + self.quiver = ax.quiver(X, Y, XY_U, XY_V, units='xy', scale=1./self.vec_scale, scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.15, linewidth=0.5, color='white',edgecolor='black') fig.canvas.draw_idle() return self.quiver else: - ax.quiver(X, Y, XY_U, XY_V, units='xy', scale=1./self.vec_scale, scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.1, linewidth=0.5, color='white',edgecolor='black') + ax.quiver(X, Y, XY_U, XY_V, units='xy', scale=1./self.vec_scale, scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.15, linewidth=0.5, color='white',edgecolor='black') fig.canvas.draw_idle() def pol_int(self, fig=None, ax=None):