From 483fd93f42eabd8a16afc3f5ba88884a9a47c69b Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Fri, 22 Mar 2024 15:28:14 +0100 Subject: [PATCH 1/2] modify interaction, fix use of global variables --- src/FOC_reduction.py | 60 +++++++++++++------------- src/lib/plots.py | 97 ++++++++++++++++++++++-------------------- src/lib/reduction.py | 84 +++++++++++++++++------------------- src/lib/utils.py | 3 +- src/overplot_IC5063.py | 88 +++++++++++++------------------------- 5 files changed, 151 insertions(+), 181 deletions(-) mode change 100644 => 100755 src/lib/utils.py diff --git a/src/FOC_reduction.py b/src/FOC_reduction.py index 2683df8..16eeaa9 100755 --- a/src/FOC_reduction.py +++ b/src/FOC_reduction.py @@ -11,10 +11,11 @@ import lib.fits as proj_fits # Functions to handle fits files import lib.reduction as proj_red # Functions used in reduction pipeline import lib.plots as proj_plots # Functions for plotting data from lib.query import retrieve_products, path_exists, system +from lib.utils import sci_not, princ_angle from matplotlib.colors import LogNorm -def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=0, interactive=0): +def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=False, interactive=False): # Reduction parameters # Deconvolution deconvolve = False @@ -34,7 +35,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= # Background estimation error_sub_type = 'freedman-diaconis' # sqrt, sturges, rice, scott, freedman-diaconis (default) or shape (example (51, 51)) subtract_error = 0.50 - display_bkg = False + display_bkg = True # Data binning rebin = True @@ -56,13 +57,9 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= rotate_data = False # rotation to North convention can give erroneous results rotate_stokes = True - # Final crop - crop = True # Crop to desired ROI - interactive = True # Whether to output to intercative analysis tool - # 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%. + SNRi_cut = 3. # 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 vec_scale = 3 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 @@ -93,6 +90,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= data_array, headers = proj_fits.get_obs_data(infiles, data_folder=data_folder, compute_flux=True) figname = "_".join([target, "FOC"]) + figtype = "" if rebin: if px_scale not in ['full']: figtype = "".join(["b", "{0:.2f}".format(pxsize), px_scale]) # additionnal informations @@ -124,8 +122,8 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= data_array, headers, error_array=error_array, background=background, upsample_factor=10, ref_center=align_center, return_shifts=False) if display_align: - proj_plots.plot_obs(data_array, headers, vmin=data_array[data_array > 0.].min( - )*headers[0]['photflam'], vmax=data_array[data_array > 0.].max()*headers[0]['photflam'], savename="_".join([figname, str(align_center)]), plots_folder=plots_folder) + proj_plots.plot_obs(data_array, headers, savename="_".join([figname, str(align_center)]), plots_folder=plots_folder, norm=LogNorm( + vmin=data_array[data_array > 0.].min()*headers[0]['photflam'], vmax=data_array[data_array > 0.].max()*headers[0]['photflam'])) # Rebin data to desired pixel size. if rebin: @@ -140,8 +138,8 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= # Plot array for checking output if display_data and px_scale.lower() not in ['full', 'integrate']: - proj_plots.plot_obs(data_array, headers, vmin=data_array[data_array > 0.].min( - )*headers[0]['photflam'], vmax=data_array[data_array > 0.].max()*headers[0]['photflam'], savename="_".join([figname, "rebin"]), plots_folder=plots_folder) + proj_plots.plot_obs(data_array, headers, savename="_".join([figname, "rebin"]), plots_folder=plots_folder, norm=LogNorm( + vmin=data_array[data_array > 0.].min()*headers[0]['photflam'], vmax=data_array[data_array > 0.].max()*headers[0]['photflam'])) background = np.array([np.array(bkg).reshape(1, 1) for bkg in background]) background_error = np.array([np.array(np.sqrt((bkg-background[np.array([h['filtnam1'] == head['filtnam1'] for h in headers], dtype=bool)].mean()) @@ -171,51 +169,52 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= # Step 4: # Save image to FITS. + figname = "_".join([figname, figtype]) if figtype != "" else figname Stokes_test = proj_fits.save_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P, - headers, data_mask, "_".join([figname, figtype]), data_folder=data_folder, return_hdul=True) + headers, data_mask, figname, data_folder=data_folder, return_hdul=True) data_mask = Stokes_test[-1].data.astype(bool) # Step 5: # crop to desired region of interest (roi) if crop: - figtype += "_crop" + figname += "_crop" stokescrop = proj_plots.crop_Stokes(deepcopy(Stokes_test), norm=LogNorm()) stokescrop.crop() - stokescrop.write_to("/".join([data_folder, "_".join([figname, figtype+".fits"])])) + stokescrop.write_to("/".join([data_folder, figname+".fits"])) Stokes_test, data_mask, headers = stokescrop.hdul_crop, stokescrop.data_mask, [dataset.header for dataset in stokescrop.hdul_crop] - print("F_int({0:.0f} Angs) = ({1} ± {2})e{3} ergs.cm^-2.s^-1.Angs^-1".format(headers[0]['photplam'], *proj_plots.sci_not( + print("F_int({0:.0f} Angs) = ({1} ± {2})e{3} ergs.cm^-2.s^-1.Angs^-1".format(headers[0]['photplam'], *sci_not( Stokes_test[0].data[data_mask].sum()*headers[0]['photflam'], np.sqrt(Stokes_test[3].data[0, 0][data_mask].sum())*headers[0]['photflam'], 2, out=int))) print("P_int = {0:.1f} ± {1:.1f} %".format(headers[0]['p_int']*100., np.ceil(headers[0]['p_int_err']*1000.)/10.)) - print("PA_int = {0:.1f} ± {1:.1f} °".format(headers[0]['pa_int'], np.ceil(headers[0]['pa_int_err']*10.)/10.)) + print("PA_int = {0:.1f} ± {1:.1f} °".format(princ_angle(headers[0]['pa_int']), princ_angle(np.ceil(headers[0]['pa_int_err']*10.)/10.))) # Background values - print("F_bkg({0:.0f} Angs) = ({1} ± {2})e{3} ergs.cm^-2.s^-1.Angs^-1".format(headers[0]['photplam'], *proj_plots.sci_not( + print("F_bkg({0:.0f} Angs) = ({1} ± {2})e{3} ergs.cm^-2.s^-1.Angs^-1".format(headers[0]['photplam'], *sci_not( I_bkg[0, 0]*headers[0]['photflam'], np.sqrt(S_cov_bkg[0, 0][0, 0])*headers[0]['photflam'], 2, out=int))) print("P_bkg = {0:.1f} ± {1:.1f} %".format(debiased_P_bkg[0, 0]*100., np.ceil(s_P_bkg[0, 0]*1000.)/10.)) - print("PA_bkg = {0:.1f} ± {1:.1f} °".format(PA_bkg[0, 0], np.ceil(s_PA_bkg[0, 0]*10.)/10.)) + print("PA_bkg = {0:.1f} ± {1:.1f} °".format(princ_angle(PA_bkg[0, 0]), princ_angle(np.ceil(s_PA_bkg[0, 0]*10.)/10.))) # Plot polarisation map (Background is either total Flux, Polarization degree or Polarization degree error). if px_scale.lower() not in ['full', 'integrate'] and not interactive: proj_plots.polarisation_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) + step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname]), plots_folder=plots_folder) proj_plots.polarisation_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') + vec_scale=vec_scale, savename="_".join([figname, "I"]), plots_folder=plots_folder, display='Intensity') proj_plots.polarisation_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') + vec_scale=vec_scale, savename="_".join([figname, "P_flux"]), plots_folder=plots_folder, display='Pol_Flux') proj_plots.polarisation_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') + vec_scale=vec_scale, savename="_".join([figname, "P"]), plots_folder=plots_folder, display='Pol_deg') proj_plots.polarisation_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') + vec_scale=vec_scale, savename="_".join([figname, "PA"]), plots_folder=plots_folder, display='Pol_ang') proj_plots.polarisation_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') + vec_scale=vec_scale, savename="_".join([figname, "I_err"]), plots_folder=plots_folder, display='I_err') proj_plots.polarisation_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') + vec_scale=vec_scale, savename="_".join([figname, "P_err"]), plots_folder=plots_folder, display='Pol_deg_err') proj_plots.polarisation_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') + vec_scale=vec_scale, savename="_".join([figname, "SNRi"]), plots_folder=plots_folder, display='SNRi') proj_plots.polarisation_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') + vec_scale=vec_scale, savename="_".join([figname, "SNRp"]), plots_folder=plots_folder, display='SNRp') elif not interactive: proj_plots.polarisation_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, - savename="_".join([figname, figtype]), plots_folder=plots_folder, display='integrate') + savename=figname, plots_folder=plots_folder, display='integrate') elif px_scale.lower() not in ['full', 'integrate']: proj_plots.pol_map(Stokes_test, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim) @@ -231,9 +230,8 @@ if __name__ == "__main__": 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") - parser.add_argument('-c', '--crop', metavar='crop_boolean', required=False, help='whether to crop the analysis region', type=int, default=0) - parser.add_argument('-i', '--interactive', metavar='interactive_boolean', required=False, - help='whether to output to the interactive analysis tool', type=int, default=0) + parser.add_argument('-c', '--crop', action='store_true', required=False, help='whether to crop the analysis region') + parser.add_argument('-i', '--interactive', action='store_true', required=False, help='whether to output to the interactive analysis tool') args = parser.parse_args() exitcode = main(target=args.target, proposal_id=args.proposal_id, infiles=args.files, output_dir=args.output_dir, crop=args.crop, interactive=args.interactive) diff --git a/src/lib/plots.py b/src/lib/plots.py index 5a141ec..604f674 100755 --- a/src/lib/plots.py +++ b/src/lib/plots.py @@ -58,8 +58,7 @@ from lib.fits import save_Stokes from lib.utils import rot2D, princ_angle, sci_not -def plot_obs(data_array, headers, shape=None, vmin=None, vmax=None, rectangle=None, - savename=None, plots_folder=""): +def plot_obs(data_array, headers, rectangle=None, savename=None, plots_folder="", **kwargs): """ Plots raw observation imagery with some information on the instrument and filters. @@ -70,10 +69,6 @@ def plot_obs(data_array, headers, shape=None, vmin=None, vmax=None, rectangle=No single observation with multiple polarizers of an instrument headers : header list List of headers corresponding to the images in data_array - shape : array-like of length 2, optional - Shape of the display, with shape = [#row, #columns]. If None, defaults - to the optimal square. - Defaults to None. vmin : float, optional Min pixel value that should be displayed. Defaults to 0. @@ -94,41 +89,49 @@ def plot_obs(data_array, headers, shape=None, vmin=None, vmax=None, rectangle=No Defaults to current folder. """ plt.rcParams.update({'font.size': 10}) - if shape is None: - shape = np.array([np.ceil(np.sqrt(data_array.shape[0])).astype(int), ]*2) + nb_obs = np.max([np.sum([head['filtnam1'] == curr_pol for head in headers]) for curr_pol in ['POL0', 'POL60', 'POL120']]) + shape = np.array((3, nb_obs)) fig, ax = plt.subplots(shape[0], shape[1], figsize=(10, 10), dpi=200, sharex=True, sharey=True) - - for i, (axe, data, head) in enumerate(zip(ax.flatten(), data_array, headers)): + r_pol = dict(pol0=0, pol60=1, pol120=2) + c_pol = dict(pol0=0, pol60=0, pol120=0) + for i, (data, head) in enumerate(zip(data_array, headers)): instr = head['instrume'] rootname = head['rootname'] exptime = head['exptime'] filt = head['filtnam1'] convert = head['photflam'] + r_ax, c_ax = r_pol[filt.lower()], c_pol[filt.lower()] + c_pol[filt.lower()] += 1 # plots - if vmin is None or vmax is None: - vmin, vmax = convert*data[data > 0.].min()/10., convert*data[data > 0.].max() - # im = axe.imshow(convert*data, vmin=vmin, vmax=vmax, origin='lower', cmap='gray') + vmin, vmax = convert*data[data > 0.].min()/10., convert*data[data > 0.].max() + for key, value in [["cmap", [["cmap", "gray"]]], ["norm", [["norm", LogNorm(vmin, vmax)]]]]: + try: + _ = kwargs[key] + except KeyError: + for key_i, val_i in value: + kwargs[key_i] = val_i + # im = ax[r_ax][c_ax].imshow(convert*data, origin='lower', **kwargs) data[data*convert < vmin*10.] = vmin*10./convert - im = axe.imshow(convert*data, norm=LogNorm(vmin, vmax), origin='lower', cmap='gray') + im = ax[r_ax][c_ax].imshow(convert*data, origin='lower', **kwargs) if rectangle is not None: x, y, width, height, angle, color = rectangle[i] - axe.add_patch(Rectangle((x, y), width, height, angle=angle, - edgecolor=color, fill=False)) + ax[r_ax][c_ax].add_patch(Rectangle((x, y), width, height, angle=angle, + edgecolor=color, fill=False)) # position of centroid - axe.plot([data.shape[1]/2, data.shape[1]/2], [0, data.shape[0]-1], '--', lw=1, - color='grey', alpha=0.5) - axe.plot([0, data.shape[1]-1], [data.shape[1]/2, data.shape[1]/2], '--', lw=1, - color='grey', alpha=0.5) - axe.annotate(instr+":"+rootname, color='white', fontsize=5, xy=(0.01, 1.00), - xycoords='axes fraction', verticalalignment='top', - horizontalalignment='left') - axe.annotate(filt, color='white', fontsize=10, xy=(0.01, 0.01), - xycoords='axes fraction', verticalalignment='bottom', - horizontalalignment='left') - axe.annotate(exptime, color='white', fontsize=5, xy=(1.00, 0.01), - xycoords='axes fraction', verticalalignment='bottom', - horizontalalignment='right') + ax[r_ax][c_ax].plot([data.shape[1]/2, data.shape[1]/2], [0, data.shape[0]-1], '--', lw=1, + color='grey', alpha=0.5) + ax[r_ax][c_ax].plot([0, data.shape[1]-1], [data.shape[1]/2, data.shape[1]/2], '--', lw=1, + color='grey', alpha=0.5) + ax[r_ax][c_ax].annotate(instr+":"+rootname, color='white', fontsize=5, xy=(0.01, 1.00), + xycoords='axes fraction', verticalalignment='top', + horizontalalignment='left') + ax[r_ax][c_ax].annotate(filt, color='white', fontsize=10, xy=(0.01, 0.01), + xycoords='axes fraction', verticalalignment='bottom', + horizontalalignment='left') + ax[r_ax][c_ax].annotate(exptime, color='white', fontsize=5, xy=(1.00, 0.01), + xycoords='axes fraction', verticalalignment='bottom', + horizontalalignment='right') fig.subplots_adjust(hspace=0.01, wspace=0.01, right=1.02) fig.colorbar(im, ax=ax[:, :], location='right', shrink=0.75, aspect=50, pad=0.025, label=r"Flux [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") @@ -306,7 +309,7 @@ def polarisation_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c vmin, vmax = flux_lim im = ax.imshow(stkI.data*convert_flux, norm=LogNorm(vmin, vmax), aspect='equal', cmap='inferno', alpha=1.) fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") - levelsI = np.array([2., 5., 10., 20., 50., 90.])/100.*vmax + levelsI = np.array([0.8, 2., 5., 10., 20., 50.])/100.*vmax print("Total flux contour levels : ", levelsI) ax.contour(stkI.data*convert_flux, levels=levelsI, colors='grey', linewidths=0.5) elif display.lower() in ['pol_flux']: @@ -391,9 +394,8 @@ def polarisation_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]") # Get integrated values from header - n_pix = stkI.data[data_mask].size I_diluted = stkI.data[data_mask].sum() - I_diluted_err = np.sqrt(n_pix)*np.sqrt(np.sum(stk_cov.data[0, 0][data_mask])) + I_diluted_err = np.sqrt(np.sum(stk_cov.data[0, 0][data_mask])) P_diluted = Stokes[0].header['P_int'] P_diluted_err = Stokes[0].header['P_int_err'] @@ -730,12 +732,10 @@ class overplot_radio(align_maps): # Display other map as contours if levels is None: - levels = np.logspace(np.log(3)/np.log(10), 2., 5)/100.*other_data[other_data > 0.].max() + levels = np.logspace(0., 1.9, 5)/100.*other_data[other_data > 0.].max() other_cont = self.ax_overplot.contour( other_data*self.other_convert, transform=self.ax_overplot.get_transform(self.other_wcs.celestial), levels=levels*self.other_convert, colors='grey') self.ax_overplot.clabel(other_cont, inline=True, fontsize=5) - other_proxy = Rectangle((0, 0), 1, 1, fc='w', ec=other_cont.collections[0].get_edgecolor()[0], label=r"{0:s} contour".format(self.other_observer)) - self.ax_overplot.add_patch(other_proxy) self.ax_overplot.set_xlabel(label="Right Ascension (J2000)") self.ax_overplot.set_ylabel(label="Declination (J2000)", labelpad=-1) @@ -761,6 +761,8 @@ class overplot_radio(align_maps): handles, labels = self.ax_overplot.get_legend_handles_labels() handles[np.argmax([li == "{0:s} polarisation map".format(self.map_observer) for li in labels]) ] = FancyArrowPatch((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])) self.legend = self.ax_overplot.legend(handles=handles, labels=labels, bbox_to_anchor=( 0., 1.02, 1., .102), loc='lower left', mode="expand", borderaxespad=0.) @@ -854,9 +856,6 @@ class overplot_chandra(align_maps): levels *= other_data.max()/self.other_data.max() other_cont = self.ax_overplot.contour(other_data*self.other_convert, transform=self.ax_overplot.get_transform(other_wcs), levels=levels, colors='grey') self.ax_overplot.clabel(other_cont, inline=True, fontsize=8) - other_proxy = Rectangle((0, 0), 1., 1., fc='w', ec=other_cont.collections[0].get_edgecolor()[ - 0], lw=2, label=r"{0:s} contour in counts".format(self.other_observer)) - self.ax_overplot.add_patch(other_proxy) self.ax_overplot.set_xlabel(label="Right Ascension (J2000)") self.ax_overplot.set_ylabel(label="Declination (J2000)", labelpad=-1) @@ -881,6 +880,8 @@ class overplot_chandra(align_maps): handles, labels = self.ax_overplot.get_legend_handles_labels() handles[np.argmax([li == "{0:s} polarisation map".format(self.map_observer) for li in labels]) ] = FancyArrowPatch((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])) self.legend = self.ax_overplot.legend(handles=handles, labels=labels, bbox_to_anchor=( 0., 1.02, 1., .102), loc='lower left', mode="expand", borderaxespad=0.) @@ -972,8 +973,6 @@ class overplot_pol(align_maps): cont_stkI = self.ax_overplot.contour(stkI*self.map_convert, levels=levels, colors='grey', alpha=0.75, transform=self.ax_overplot.get_transform(self.wcs_UV)) # self.ax_overplot.clabel(cont_stkI, inline=True, fontsize=5) - cont_proxy = Rectangle((0, 0), 1, 1, fc='w', ec=cont_stkI.collections[0].get_edgecolor()[0], label="{0:s} Stokes I contour".format(self.map_observer)) - self.ax_overplot.add_patch(cont_proxy) # Display pixel scale and North direction fontprops = fm.FontProperties(size=16) @@ -1001,6 +1000,8 @@ class overplot_pol(align_maps): handles, labels = self.ax_overplot.get_legend_handles_labels() handles[np.argmax([li == "{0:s} polarisation map".format(self.map_observer) for li in labels]) ] = FancyArrowPatch((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])) self.legend = self.ax_overplot.legend(handles=handles, labels=labels, bbox_to_anchor=( 0., 1.02, 1., .102), loc='lower left', mode="expand", borderaxespad=0.) @@ -2351,10 +2352,9 @@ class pol_map(object): def pol_int(self, fig=None, ax=None): if self.region is None: - n_pix = self.I.size s_I = np.sqrt(self.IQU_cov[0, 0]) I_reg = self.I.sum() - I_reg_err = np.sqrt(n_pix)*np.sqrt(np.sum(s_I**2)) + I_reg_err = np.sqrt(np.sum(s_I**2)) P_reg = self.Stokes[0].header['P_int'] P_reg_err = self.Stokes[0].header['P_int_err'] PA_reg = self.Stokes[0].header['PA_int'] @@ -2386,7 +2386,6 @@ class pol_map(object): Q_cut_err**2 + Q_cut**2*U_cut_err**2 - 2.*Q_cut*U_cut*QU_cut_err))) else: - n_pix = self.I[self.region].size s_I = np.sqrt(self.IQU_cov[0, 0]) s_Q = np.sqrt(self.IQU_cov[1, 1]) s_U = np.sqrt(self.IQU_cov[2, 2]) @@ -2442,15 +2441,19 @@ class pol_map(object): ax = self.ax if hasattr(self, 'an_int'): self.an_int.remove() - self.an_int = ax.annotate(r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(self.pivot_wav, sci_not(I_reg*self.map_convert, I_reg_err*self.map_convert, 2))+"\n"+r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_reg*100., np.ceil( - P_reg_err*1000.)/10.)+"\n"+r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_reg, np.ceil(PA_reg_err*10.)/10.), color='white', fontsize=12, xy=(0.01, 1.00), xycoords='axes fraction', path_effects=[pe.withStroke(linewidth=0.5, foreground='k')], verticalalignment='top', horizontalalignment='left') + self.str_int = r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(self.pivot_wav, sci_not(I_reg*self.map_convert, I_reg_err*self.map_convert, 2))+"\n"+r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_reg*100., np.ceil(P_reg_err*1000.)/10.)+"\n"+r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_reg, np.ceil(PA_reg_err*10.)/10.) + self.str_cut = "" + # self.str_cut = "\n"+r"$F_{{\lambda}}^{{cut}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(self.pivot_wav, sci_not(I_cut*self.map_convert, I_cut_err*self.map_convert, 2))+"\n"+r"$P^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_cut*100., np.ceil(P_cut_err*1000.)/10.)+"\n"+r"$\theta_{{P}}^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_cut, np.ceil(PA_cut_err*10.)/10.) + self.an_int = ax.annotate(self.str_int+self.str_cut, color='white', fontsize=12, xy=(0.01, 1.00), xycoords='axes fraction', path_effects=[pe.withStroke(linewidth=0.5, foreground='k')], verticalalignment='top', horizontalalignment='left') if self.region is not None: self.cont = ax.contour(self.region.astype(float), levels=[0.5], colors='white', linewidths=0.8) fig.canvas.draw_idle() return self.an_int else: - ax.annotate(r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(self.pivot_wav, sci_not(I_reg*self.map_convert, I_reg_err*self.map_convert, 2))+"\n"+r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_reg*100., np.ceil(P_reg_err*1000.)/10.) + - "\n"+r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_reg, np.ceil(PA_reg_err*10.)/10.), color='white', fontsize=12, xy=(0.01, 1.00), xycoords='axes fraction', path_effects=[pe.withStroke(linewidth=0.5, foreground='k')], verticalalignment='top', horizontalalignment='left') + str_int = r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(self.pivot_wav, sci_not(I_reg*self.map_convert, I_reg_err*self.map_convert, 2))+"\n"+r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_reg*100., np.ceil(P_reg_err*1000.)/10.)+"\n"+r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_reg, np.ceil(PA_reg_err*10.)/10.) + str_cut = "" + # str_cut = "\n"+r"$F_{{\lambda}}^{{cut}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(self.pivot_wav, sci_not(I_cut*self.map_convert, I_cut_err*self.map_convert, 2))+"\n"+r"$P^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_cut*100., np.ceil(P_cut_err*1000.)/10.)+"\n"+r"$\theta_{{P}}^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_cut, np.ceil(PA_cut_err*10.)/10.) + ax.annotate(str_int+str_cut, color='white', fontsize=12, xy=(0.01, 1.00), xycoords='axes fraction', path_effects=[pe.withStroke(linewidth=0.5, foreground='k')], verticalalignment='top', horizontalalignment='left') if self.region is not None: ax.contour(self.region.astype(float), levels=[0.5], colors='white', linewidths=0.8) fig.canvas.draw_idle() diff --git a/src/lib/reduction.py b/src/lib/reduction.py index 6f504bf..6019b98 100755 --- a/src/lib/reduction.py +++ b/src/lib/reduction.py @@ -12,7 +12,7 @@ prototypes : Homogeneously deconvolve a data array using a chosen deconvolution algorithm. - get_error(data_array, headers, error_array, data_mask, sub_type, display, savename, plots_folder, return_background) -> data_array, error_array, headers (, background) - Compute the error (noise) on each image of the input array. + Compute the uncertainties from the different error sources on each image of the input array. - rebin_array(data_array, error_array, headers, pxsize, scale, operation, data_mask) -> rebinned_data, rebinned_error, rebinned_headers, Dxy (, data_mask) Homegeneously rebin a data array given a target pixel size in scale units. @@ -729,8 +729,8 @@ def align_data(data_array, headers, error_array=None, background=None, upsample_ if do_shift: shift, error, _ = phase_cross_correlation(ref_data/ref_data.max(), image/image.max(), upsample_factor=upsample_factor) else: - shift = pol_shift[headers[i]['filtnam1'].lower()] - error = sigma_shift[headers[i]['filtnam1'].lower()] + shift = globals["pol_shift"][headers[i]['filtnam1'].lower()] + error = globals["sigma_shift"][headers[i]['filtnam1'].lower()] # Rescale image to requested output rescaled_image[i, res_shift[0]:res_shift[0]+shape[1], res_shift[1]:res_shift[1]+shape[2]] = deepcopy(image) rescaled_error[i, res_shift[0]:res_shift[0]+shape[1], res_shift[1]:res_shift[1]+shape[2]] = deepcopy(error_array[i]) @@ -1101,16 +1101,16 @@ def compute_Stokes(data_array, error_array, data_mask, headers, FWHM=None, scale same_filt3 = np.array([filt3 == header['filtnam3'] for header in headers]).all() same_filt4 = np.array([filt4 == header['filtnam4'] for header in headers]).all() if (same_filt2 and same_filt3 and same_filt4): - transmit2, transmit3, transmit4 = trans2[filt2.lower()], trans3[filt3.lower()], trans4[filt4.lower()] + transmit2, transmit3, transmit4 = globals["trans2"][filt2.lower()], globals["trans3"][filt3.lower()], globals["trans4"][filt4.lower()] else: print("WARNING : All images in data_array are not from the same \ band filter, the limiting transmittance will be taken.") - transmit2 = np.min([trans2[header['filtnam2'].lower()] for header in headers]) - transmit3 = np.min([trans3[header['filtnam3'].lower()] for header in headers]) - transmit4 = np.min([trans4[header['filtnam4'].lower()] for header in headers]) + transmit2 = np.min([globals["trans2"][header['filtnam2'].lower()] for header in headers]) + transmit3 = np.min([globals["trans3"][header['filtnam3'].lower()] for header in headers]) + transmit4 = np.min([globals["trans4"][header['filtnam4'].lower()] for header in headers]) if transmitcorr: transmit *= transmit2*transmit3*transmit4 - pol_eff = np.array([pol_efficiency['pol0'], pol_efficiency['pol60'], pol_efficiency['pol120']]) + pol_eff = np.array([globals["pol_efficiency"]['pol0'], globals["pol_efficiency"]['pol60'], globals["pol_efficiency"]['pol120']]) # Calculating correction factor corr = np.array([1.0*h['photflam']/h['exptime'] for h in pol_headers])*pol_headers[0]['exptime']/pol_headers[0]['photflam'] @@ -1122,9 +1122,11 @@ def compute_Stokes(data_array, error_array, data_mask, headers, FWHM=None, scale coeff_stokes = np.zeros((3, 3)) # Coefficients linking each polarizer flux to each Stokes parameter for i in range(3): - coeff_stokes[0, i] = pol_eff[(i+1) % 3]*pol_eff[(i+2) % 3]*np.sin(-2.*theta[(i+1) % 3]+2.*theta[(i+2) % 3])*2./transmit[i] - coeff_stokes[1, i] = (-pol_eff[(i+1) % 3]*np.sin(2.*theta[(i+1) % 3]) + pol_eff[(i+2) % 3]*np.sin(2.*theta[(i+2) % 3]))*2./transmit[i] - coeff_stokes[2, i] = (pol_eff[(i+1) % 3]*np.cos(2.*theta[(i+1) % 3]) - pol_eff[(i+2) % 3]*np.cos(2.*theta[(i+2) % 3]))*2./transmit[i] + coeff_stokes[0, i] = pol_eff[(i+1) % 3]*pol_eff[(i+2) % 3]*np.sin(-2.*globals["theta"][(i+1) % 3]+2.*globals["theta"][(i+2) % 3])*2./transmit[i] + coeff_stokes[1, i] = (-pol_eff[(i+1) % 3]*np.sin(2.*globals["theta"][(i+1) % 3]) + + pol_eff[(i+2) % 3]*np.sin(2.*globals["theta"][(i+2) % 3]))*2./transmit[i] + coeff_stokes[2, i] = (pol_eff[(i+1) % 3]*np.cos(2.*globals["theta"][(i+1) % 3]) - + pol_eff[(i+2) % 3]*np.cos(2.*globals["theta"][(i+2) % 3]))*2./transmit[i] # Normalization parameter for Stokes parameters computation A = (coeff_stokes[0, :]*transmit/2.).sum() @@ -1173,34 +1175,34 @@ def compute_Stokes(data_array, error_array, data_mask, headers, FWHM=None, scale s_U2_stat = np.sum([coeff_stokes[2, i]**2*sigma_flux[i]**2 for i in range(len(sigma_flux))], axis=0) # Compute the derivative of each Stokes parameter with respect to the polarizer orientation - dI_dtheta1 = 2.*pol_eff[0]/A*(pol_eff[2]*np.cos(-2.*theta[2]+2.*theta[0])*(pol_flux[1]-I_stokes) - - pol_eff[1]*np.cos(-2.*theta[0]+2.*theta[1])*(pol_flux[2]-I_stokes)) - dI_dtheta2 = 2.*pol_eff[1]/A*(pol_eff[0]*np.cos(-2.*theta[0]+2.*theta[1])*(pol_flux[2]-I_stokes) - - pol_eff[2]*np.cos(-2.*theta[1]+2.*theta[2])*(pol_flux[0]-I_stokes)) - dI_dtheta3 = 2.*pol_eff[2]/A*(pol_eff[1]*np.cos(-2.*theta[1]+2.*theta[2])*(pol_flux[0]-I_stokes) - - pol_eff[0]*np.cos(-2.*theta[2]+2.*theta[0])*(pol_flux[1]-I_stokes)) + dI_dtheta1 = 2.*pol_eff[0]/A*(pol_eff[2]*np.cos(-2.*globals()["theta"][2]+2.*globals()["theta"][0])*(pol_flux[1]-I_stokes) - + pol_eff[1]*np.cos(-2.*globals()["theta"][0]+2.*globals()["theta"][1])*(pol_flux[2]-I_stokes)) + dI_dtheta2 = 2.*pol_eff[1]/A*(pol_eff[0]*np.cos(-2.*globals()["theta"][0]+2.*globals()["theta"][1])*(pol_flux[2]-I_stokes) - + pol_eff[2]*np.cos(-2.*globals()["theta"][1]+2.*globals()["theta"][2])*(pol_flux[0]-I_stokes)) + dI_dtheta3 = 2.*pol_eff[2]/A*(pol_eff[1]*np.cos(-2.*globals()["theta"][1]+2.*globals()["theta"][2])*(pol_flux[0]-I_stokes) - + pol_eff[0]*np.cos(-2.*globals()["theta"][2]+2.*globals()["theta"][0])*(pol_flux[1]-I_stokes)) dI_dtheta = np.array([dI_dtheta1, dI_dtheta2, dI_dtheta3]) - dQ_dtheta1 = 2.*pol_eff[0]/A*(np.cos(2.*theta[0])*(pol_flux[1]-pol_flux[2]) - (pol_eff[2]*np.cos(-2. * - theta[2]+2.*theta[0]) - pol_eff[1]*np.cos(-2.*theta[0]+2.*theta[1]))*Q_stokes) - dQ_dtheta2 = 2.*pol_eff[1]/A*(np.cos(2.*theta[1])*(pol_flux[2]-pol_flux[0]) - (pol_eff[0]*np.cos(-2. * - theta[0]+2.*theta[1]) - pol_eff[2]*np.cos(-2.*theta[1]+2.*theta[2]))*Q_stokes) - dQ_dtheta3 = 2.*pol_eff[2]/A*(np.cos(2.*theta[2])*(pol_flux[0]-pol_flux[1]) - (pol_eff[1]*np.cos(-2. * - theta[1]+2.*theta[2]) - pol_eff[0]*np.cos(-2.*theta[2]+2.*theta[0]))*Q_stokes) + dQ_dtheta1 = 2.*pol_eff[0]/A*(np.cos(2.*globals()["theta"][0])*(pol_flux[1]-pol_flux[2]) - (pol_eff[2]*np.cos(-2.*globals() + ["theta"][2]+2.*globals()["theta"][0]) - pol_eff[1]*np.cos(-2.*globals()["theta"][0]+2.*globals()["theta"][1]))*Q_stokes) + dQ_dtheta2 = 2.*pol_eff[1]/A*(np.cos(2.*globals()["theta"][1])*(pol_flux[2]-pol_flux[0]) - (pol_eff[0]*np.cos(-2.*globals() + ["theta"][0]+2.*globals()["theta"][1]) - pol_eff[2]*np.cos(-2.*globals()["theta"][1]+2.*globals()["theta"][2]))*Q_stokes) + dQ_dtheta3 = 2.*pol_eff[2]/A*(np.cos(2.*globals()["theta"][2])*(pol_flux[0]-pol_flux[1]) - (pol_eff[1]*np.cos(-2.*globals() + ["theta"][1]+2.*globals()["theta"][2]) - pol_eff[0]*np.cos(-2.*globals()["theta"][2]+2.*globals()["theta"][0]))*Q_stokes) dQ_dtheta = np.array([dQ_dtheta1, dQ_dtheta2, dQ_dtheta3]) - dU_dtheta1 = 2.*pol_eff[0]/A*(np.sin(2.*theta[0])*(pol_flux[1]-pol_flux[2]) - (pol_eff[2]*np.cos(-2. * - theta[2]+2.*theta[0]) - pol_eff[1]*np.cos(-2.*theta[0]+2.*theta[1]))*U_stokes) - dU_dtheta2 = 2.*pol_eff[1]/A*(np.sin(2.*theta[1])*(pol_flux[2]-pol_flux[0]) - (pol_eff[0]*np.cos(-2. * - theta[0]+2.*theta[1]) - pol_eff[2]*np.cos(-2.*theta[1]+2.*theta[2]))*U_stokes) - dU_dtheta3 = 2.*pol_eff[2]/A*(np.sin(2.*theta[2])*(pol_flux[0]-pol_flux[1]) - (pol_eff[1]*np.cos(-2. * - theta[1]+2.*theta[2]) - pol_eff[0]*np.cos(-2.*theta[2]+2.*theta[0]))*U_stokes) + dU_dtheta1 = 2.*pol_eff[0]/A*(np.sin(2.*globals()["theta"][0])*(pol_flux[1]-pol_flux[2]) - (pol_eff[2]*np.cos(-2.*globals() + ["theta"][2]+2.*globals()["theta"][0]) - pol_eff[1]*np.cos(-2.*globals()["theta"][0]+2.*globals()["theta"][1]))*U_stokes) + dU_dtheta2 = 2.*pol_eff[1]/A*(np.sin(2.*globals()["theta"][1])*(pol_flux[2]-pol_flux[0]) - (pol_eff[0]*np.cos(-2.*globals() + ["theta"][0]+2.*globals()["theta"][1]) - pol_eff[2]*np.cos(-2.*globals()["theta"][1]+2.*globals()["theta"][2]))*U_stokes) + dU_dtheta3 = 2.*pol_eff[2]/A*(np.sin(2.*globals()["theta"][2])*(pol_flux[0]-pol_flux[1]) - (pol_eff[1]*np.cos(-2.*globals() + ["theta"][1]+2.*globals()["theta"][2]) - pol_eff[0]*np.cos(-2.*globals()["theta"][2]+2.*globals()["theta"][0]))*U_stokes) dU_dtheta = np.array([dU_dtheta1, dU_dtheta2, dU_dtheta3]) # Compute the uncertainty associated with the polarizers' orientation (see Kishimoto 1999) - s_I2_axis = np.sum([dI_dtheta[i]**2 * sigma_theta[i]**2 for i in range(len(sigma_theta))], axis=0) - s_Q2_axis = np.sum([dQ_dtheta[i]**2 * sigma_theta[i]**2 for i in range(len(sigma_theta))], axis=0) - s_U2_axis = np.sum([dU_dtheta[i]**2 * sigma_theta[i]**2 for i in range(len(sigma_theta))], axis=0) + s_I2_axis = np.sum([dI_dtheta[i]**2 * globals()["sigma_theta"][i]**2 for i in range(len(globals()["sigma_theta"]))], axis=0) + s_Q2_axis = np.sum([dQ_dtheta[i]**2 * globals()["sigma_theta"][i]**2 for i in range(len(globals()["sigma_theta"]))], axis=0) + s_U2_axis = np.sum([dU_dtheta[i]**2 * globals()["sigma_theta"][i]**2 for i in range(len(globals()["sigma_theta"]))], axis=0) # np.savetxt("output/sI_dir.txt", np.sqrt(s_I2_axis)) # np.savetxt("output/sQ_dir.txt", np.sqrt(s_Q2_axis)) # np.savetxt("output/sU_dir.txt", np.sqrt(s_U2_axis)) @@ -1212,7 +1214,6 @@ def compute_Stokes(data_array, error_array, data_mask, headers, FWHM=None, scale # Compute integrated values for P, PA before any rotation mask = np.logical_and(data_mask.astype(bool), (I_stokes > 0.)) - n_pix = I_stokes[mask].size I_diluted = I_stokes[mask].sum() Q_diluted = Q_stokes[mask].sum() U_diluted = U_stokes[mask].sum() @@ -1224,12 +1225,10 @@ def compute_Stokes(data_array, error_array, data_mask, headers, FWHM=None, scale QU_diluted_err = np.sqrt(np.sum(Stokes_cov[1, 2][mask]**2)) P_diluted = np.sqrt(Q_diluted**2+U_diluted**2)/I_diluted - P_diluted_err = (1./I_diluted)*np.sqrt((Q_diluted**2*Q_diluted_err**2 + U_diluted**2*U_diluted_err**2 + 2.*Q_diluted*U_diluted*QU_diluted_err)/(Q_diluted**2 + U_diluted ** - 2) + ((Q_diluted/I_diluted)**2 + (U_diluted/I_diluted)**2)*I_diluted_err**2 - 2.*(Q_diluted/I_diluted)*IQ_diluted_err - 2.*(U_diluted/I_diluted)*IU_diluted_err) + P_diluted_err = (1./I_diluted)*np.sqrt((Q_diluted**2*Q_diluted_err**2 + U_diluted**2*U_diluted_err**2 + 2.*Q_diluted*U_diluted*QU_diluted_err)/(Q_diluted**2 + U_diluted**2) + ((Q_diluted/I_diluted)**2 + (U_diluted/I_diluted)**2)*I_diluted_err**2 - 2.*(Q_diluted/I_diluted)*IQ_diluted_err - 2.*(U_diluted/I_diluted)*IU_diluted_err) PA_diluted = princ_angle((90./np.pi)*np.arctan2(U_diluted, Q_diluted)) - PA_diluted_err = (90./(np.pi*(Q_diluted**2 + U_diluted**2)))*np.sqrt(U_diluted**2*Q_diluted_err ** - 2 + Q_diluted**2*U_diluted_err**2 - 2.*Q_diluted*U_diluted*QU_diluted_err) + PA_diluted_err = (90./(np.pi*(Q_diluted**2 + U_diluted**2)))*np.sqrt(U_diluted**2*Q_diluted_err**2 + Q_diluted**2*U_diluted_err**2 - 2.*Q_diluted*U_diluted*QU_diluted_err) for header in headers: header['P_int'] = (P_diluted, 'Integrated polarisation degree') @@ -1307,7 +1306,7 @@ def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers): s_PA[np.isnan(s_PA)] = fmax # Catch expected "OverflowWarning" as wrong pixel have an overflowing error - with warnings.catch_warnings(record=True) as w: + with warnings.catch_warnings(record=True) as _: mask2 = P**2 >= s_P**2 debiased_P = np.zeros(I_stokes.shape) debiased_P[mask2] = np.sqrt(P[mask2]**2 - s_P[mask2]**2) @@ -1474,7 +1473,6 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers, # Compute updated integrated values for P, PA mask = deepcopy(new_data_mask).astype(bool) - n_pix = new_I_stokes[mask].size I_diluted = new_I_stokes[mask].sum() Q_diluted = new_Q_stokes[mask].sum() U_diluted = new_U_stokes[mask].sum() @@ -1486,12 +1484,10 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers, QU_diluted_err = np.sqrt(np.sum(new_Stokes_cov[1, 2][mask]**2)) P_diluted = np.sqrt(Q_diluted**2+U_diluted**2)/I_diluted - P_diluted_err = (1./I_diluted)*np.sqrt((Q_diluted**2*Q_diluted_err**2 + U_diluted**2*U_diluted_err**2 + 2.*Q_diluted*U_diluted*QU_diluted_err)/(Q_diluted**2 + U_diluted ** - 2) + ((Q_diluted/I_diluted)**2 + (U_diluted/I_diluted)**2)*I_diluted_err**2 - 2.*(Q_diluted/I_diluted)*IQ_diluted_err - 2.*(U_diluted/I_diluted)*IU_diluted_err) + P_diluted_err = (1./I_diluted)*np.sqrt((Q_diluted**2*Q_diluted_err**2 + U_diluted**2*U_diluted_err**2 + 2.*Q_diluted*U_diluted*QU_diluted_err)/(Q_diluted**2 + U_diluted**2) + ((Q_diluted/I_diluted)**2 + (U_diluted/I_diluted)**2)*I_diluted_err**2 - 2.*(Q_diluted/I_diluted)*IQ_diluted_err - 2.*(U_diluted/I_diluted)*IU_diluted_err) PA_diluted = princ_angle((90./np.pi)*np.arctan2(U_diluted, Q_diluted)) - PA_diluted_err = (90./(np.pi*(Q_diluted**2 + U_diluted**2)))*np.sqrt(U_diluted**2*Q_diluted_err ** - 2 + Q_diluted**2*U_diluted_err**2 - 2.*Q_diluted*U_diluted*QU_diluted_err) + PA_diluted_err = (90./(np.pi*(Q_diluted**2 + U_diluted**2)))*np.sqrt(U_diluted**2*Q_diluted_err**2 + Q_diluted**2*U_diluted_err**2 - 2.*Q_diluted*U_diluted*QU_diluted_err) for header in new_headers: header['P_int'] = (P_diluted, 'Integrated polarisation degree') @@ -1573,6 +1569,6 @@ def rotate_data(data_array, error_array, data_mask, headers, ang): new_header[key] = val new_headers.append(new_header) - globals()['theta'] = theta - alpha + globals()['theta'] = globals()["theta"] - alpha return new_data_array, new_error_array, new_data_mask, new_headers diff --git a/src/lib/utils.py b/src/lib/utils.py old mode 100644 new mode 100755 index c668142..51a4568 --- a/src/lib/utils.py +++ b/src/lib/utils.py @@ -10,7 +10,8 @@ def rot2D(ang): def princ_angle(ang): """ - Return the principal angle in the 0° to 360° quadrant. + Return the principal angle in the 0° to 180° quadrant as PA is always + defined at p/m 180°. """ if not isinstance(ang, np.ndarray): A = np.array([ang]) diff --git a/src/overplot_IC5063.py b/src/overplot_IC5063.py index 61d7962..3cb55c3 100755 --- a/src/overplot_IC5063.py +++ b/src/overplot_IC5063.py @@ -1,72 +1,44 @@ #!/usr/bin/python3 -from os import system as command from astropy.io import fits import numpy as np -from copy import deepcopy -from lib.plots import overplot_radio, overplot_pol, align_pol +from lib.plots import overplot_radio, overplot_pol 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_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.,1.97,5)/100. -# -# 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_forced.pdf',vec_scale=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_forced.pdf',vec_scale=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_forced.pdf',vec_scale=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_forced.pdf',vec_scale=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_forced.pdf',vec_scale=None) -## +levelsMorganti = np.logspace(-0.1249, 1.97, 7)/100. + +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) + +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) + +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) + +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) + +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) + # F = overplot_pol(Stokes_UV, Stokes_S2) -# F.plot(SNRp_cut=3.0, SNRi_cut=80.0, savename='./plots/IC5063/S2_overplot_forced.pdf', norm=LogNorm(vmin=5e-20,vmax=5e-18)) +# F.plot(SNRp_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, SNRi_cut=10.0, savename='./plots/IC5063/IR_overplot_forced.pdf', vec_scale=None, +G.plot(SNRp_cut=2.0, SNRi_cut=10.0, savename='./plots/IC5063/IR_overplot.pdf', vec_scale=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') - -# data_folder1 = "./data/M87/POS1/" -# plots_folder1 = "./plots/M87/POS1/" -# basename1 = "M87_020_log" -# M87_1_95 = fits.open(data_folder1+"M87_POS1_1995_FOC_combine_FWHM020.fits") -# M87_1_96 = fits.open(data_folder1+"M87_POS1_1996_FOC_combine_FWHM020.fits") -# M87_1_97 = fits.open(data_folder1+"M87_POS1_1997_FOC_combine_FWHM020.fits") -# M87_1_98 = fits.open(data_folder1+"M87_POS1_1998_FOC_combine_FWHM020.fits") -# M87_1_99 = fits.open(data_folder1+"M87_POS1_1999_FOC_combine_FWHM020.fits") - -# H = align_pol(np.array([M87_1_95,M87_1_96,M87_1_97,M87_1_98,M87_1_99]), norm=LogNorm()) -# H.plot(SNRp_cut=5.0, SNRi_cut=50.0, savename=plots_folder1+'animated_loop/'+basename1, norm=LogNorm()) -# command("convert -delay 50 -loop 0 {0:s}animated_loop/{1:s}*.pdf {0:s}animated_loop/{1:s}.gif".format(plots_folder1, basename1)) - -# data_folder3 = "./data/M87/POS3/" -# plots_folder3 = "./plots/M87/POS3/" -# basename3 = "M87_020_log" -# M87_3_95 = fits.open(data_folder3+"M87_POS3_1995_FOC_combine_FWHM020.fits") -# M87_3_96 = fits.open(data_folder3+"M87_POS3_1996_FOC_combine_FWHM020.fits") -# M87_3_97 = fits.open(data_folder3+"M87_POS3_1997_FOC_combine_FWHM020.fits") -# M87_3_98 = fits.open(data_folder3+"M87_POS3_1998_FOC_combine_FWHM020.fits") -# M87_3_99 = fits.open(data_folder3+"M87_POS3_1999_FOC_combine_FWHM020.fits") - -# I = align_pol(np.array([M87_3_95,M87_3_96,M87_3_97,M87_3_98,M87_3_99]), norm=LogNorm()) -# I.plot(SNRp_cut=5.0, SNRi_cut=50.0, savename=plots_folder3+'animated_loop/'+basename3, norm=LogNorm()) -# command("convert -delay 20 -loop 0 {0:s}animated_loop/{1:s}*.pdf {0:s}animated_loop/{1:s}.gif".format(plots_folder3, basename3)) From 26fb83dd182c632fa499339803d39b6d54fe77fa Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Tue, 26 Mar 2024 18:03:57 +0100 Subject: [PATCH 2/2] fix global variables --- src/FOC_reduction.py | 2 +- src/lib/plots.py | 5 ++--- src/lib/reduction.py | 20 ++++++++++---------- 3 files changed, 13 insertions(+), 14 deletions(-) diff --git a/src/FOC_reduction.py b/src/FOC_reduction.py index 16eeaa9..3700749 100755 --- a/src/FOC_reduction.py +++ b/src/FOC_reduction.py @@ -35,7 +35,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= # Background estimation error_sub_type = 'freedman-diaconis' # sqrt, sturges, rice, scott, freedman-diaconis (default) or shape (example (51, 51)) subtract_error = 0.50 - display_bkg = True + display_bkg = False # Data binning rebin = True diff --git a/src/lib/plots.py b/src/lib/plots.py index 604f674..6675845 100755 --- a/src/lib/plots.py +++ b/src/lib/plots.py @@ -54,7 +54,6 @@ from astropy.wcs import WCS from astropy.io import fits from astropy.coordinates import SkyCoord from scipy.ndimage import zoom as sc_zoom -from lib.fits import save_Stokes from lib.utils import rot2D, princ_angle, sci_not @@ -416,7 +415,7 @@ def polarisation_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c X, Y = np.meshgrid(np.arange(stkI.data.shape[1]), np.arange(stkI.data.shape[0])) U, V = poldata*np.cos(np.pi/2.+pangdata*np.pi/180.), poldata*np.sin(np.pi/2.+pangdata*np.pi/180.) ax.quiver(X[::step_vec, ::step_vec], Y[::step_vec, ::step_vec], U[::step_vec, ::step_vec], V[::step_vec, ::step_vec], units='xy', angles='uv', - scale=1./vec_scale, scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.1, linewidth=0.5, color='w', edgecolor='k') + scale=1./vec_scale, scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.2, linewidth=0.3, color='w', edgecolor='k') pol_sc = AnchoredSizeBar(ax.transData, vec_scale, r"$P$= 100 %", 4, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color='w') ax.add_artist(pol_sc) @@ -448,7 +447,7 @@ def polarisation_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c if savename is not None: if savename[-4:] not in ['.png', '.jpg', '.pdf']: savename += '.pdf' - fig.savefig(path_join(plots_folder, savename), bbox_inches='tight', dpi=300) + fig.savefig(path_join(plots_folder, savename), bbox_inches='tight', dpi=200) plt.show() return fig, ax diff --git a/src/lib/reduction.py b/src/lib/reduction.py index 6019b98..582b0a4 100755 --- a/src/lib/reduction.py +++ b/src/lib/reduction.py @@ -1101,16 +1101,16 @@ def compute_Stokes(data_array, error_array, data_mask, headers, FWHM=None, scale same_filt3 = np.array([filt3 == header['filtnam3'] for header in headers]).all() same_filt4 = np.array([filt4 == header['filtnam4'] for header in headers]).all() if (same_filt2 and same_filt3 and same_filt4): - transmit2, transmit3, transmit4 = globals["trans2"][filt2.lower()], globals["trans3"][filt3.lower()], globals["trans4"][filt4.lower()] + transmit2, transmit3, transmit4 = globals()["trans2"][filt2.lower()], globals()["trans3"][filt3.lower()], globals()["trans4"][filt4.lower()] else: print("WARNING : All images in data_array are not from the same \ band filter, the limiting transmittance will be taken.") - transmit2 = np.min([globals["trans2"][header['filtnam2'].lower()] for header in headers]) - transmit3 = np.min([globals["trans3"][header['filtnam3'].lower()] for header in headers]) - transmit4 = np.min([globals["trans4"][header['filtnam4'].lower()] for header in headers]) + transmit2 = np.min([globals()["trans2"][header['filtnam2'].lower()] for header in headers]) + transmit3 = np.min([globals()["trans3"][header['filtnam3'].lower()] for header in headers]) + transmit4 = np.min([globals()["trans4"][header['filtnam4'].lower()] for header in headers]) if transmitcorr: transmit *= transmit2*transmit3*transmit4 - pol_eff = np.array([globals["pol_efficiency"]['pol0'], globals["pol_efficiency"]['pol60'], globals["pol_efficiency"]['pol120']]) + pol_eff = np.array([globals()["pol_efficiency"]['pol0'], globals()["pol_efficiency"]['pol60'], globals()["pol_efficiency"]['pol120']]) # Calculating correction factor corr = np.array([1.0*h['photflam']/h['exptime'] for h in pol_headers])*pol_headers[0]['exptime']/pol_headers[0]['photflam'] @@ -1122,11 +1122,11 @@ def compute_Stokes(data_array, error_array, data_mask, headers, FWHM=None, scale coeff_stokes = np.zeros((3, 3)) # Coefficients linking each polarizer flux to each Stokes parameter for i in range(3): - coeff_stokes[0, i] = pol_eff[(i+1) % 3]*pol_eff[(i+2) % 3]*np.sin(-2.*globals["theta"][(i+1) % 3]+2.*globals["theta"][(i+2) % 3])*2./transmit[i] - coeff_stokes[1, i] = (-pol_eff[(i+1) % 3]*np.sin(2.*globals["theta"][(i+1) % 3]) + - pol_eff[(i+2) % 3]*np.sin(2.*globals["theta"][(i+2) % 3]))*2./transmit[i] - coeff_stokes[2, i] = (pol_eff[(i+1) % 3]*np.cos(2.*globals["theta"][(i+1) % 3]) - - pol_eff[(i+2) % 3]*np.cos(2.*globals["theta"][(i+2) % 3]))*2./transmit[i] + coeff_stokes[0, i] = pol_eff[(i+1) % 3]*pol_eff[(i+2) % 3]*np.sin(-2.*globals()["theta"][(i+1) % 3]+2.*globals()["theta"][(i+2) % 3])*2./transmit[i] + coeff_stokes[1, i] = (-pol_eff[(i+1) % 3]*np.sin(2.*globals()["theta"][(i+1) % 3]) + + pol_eff[(i+2) % 3]*np.sin(2.*globals()["theta"][(i+2) % 3]))*2./transmit[i] + coeff_stokes[2, i] = (pol_eff[(i+1) % 3]*np.cos(2.*globals()["theta"][(i+1) % 3]) - + pol_eff[(i+2) % 3]*np.cos(2.*globals()["theta"][(i+2) % 3]))*2./transmit[i] # Normalization parameter for Stokes parameters computation A = (coeff_stokes[0, :]*transmit/2.).sum()