From fc80bafc665747212dbdf0c9946786b24615aee0 Mon Sep 17 00:00:00 2001 From: Tibeuleu Date: Tue, 31 Jan 2023 17:39:25 +0100 Subject: [PATCH] general improvements to pipeline --- src/FOC_reduction.py | 20 ++++++++++---------- src/lib/background.py | 4 ++-- src/lib/fits.py | 2 ++ src/lib/plots.py | 35 +++++++++++++++++++++-------------- 4 files changed, 35 insertions(+), 26 deletions(-) diff --git a/src/FOC_reduction.py b/src/FOC_reduction.py index 6476085..47fe4ae 100755 --- a/src/FOC_reduction.py +++ b/src/FOC_reduction.py @@ -131,19 +131,19 @@ def main(): display_crop = False # Error estimation error_sub_type = 'freedman-diaconis' #sqrt, sturges, rice, scott, freedman-diaconis (default) or shape (example (15,15)) - subtract_error = False + subtract_error = True display_error = False # Data binning rebin = True - pxsize = 10 - px_scale = 'pixel' #pixel, arcsec or full + pxsize = 0.10 + px_scale = 'arcsec' #pixel, arcsec or full rebin_operation = 'sum' #sum or average # Alignement align_center = 'image' #If None will align image to image center display_data = False # Smoothing - smoothing_function = 'combine' #gaussian_after, weighted_gaussian_after, gaussian, weighted_gaussian or combine - smoothing_FWHM = None #If None, no smoothing is done + smoothing_function = 'weighted_gaussian' #gaussian_after, weighted_gaussian_after, gaussian, weighted_gaussian or combine + smoothing_FWHM = 0.20 #If None, no smoothing is done smoothing_scale = 'arcsec' #pixel or arcsec # Rotation rotate_stokes = True @@ -153,7 +153,7 @@ def main(): final_display = True # Polarization map output figname = 'NGC1068_FOC' #target/intrument name - figtype = '_bin10px' #additionnal informations + figtype = '_wg_020' #additionnal informations SNRp_cut = 5. #P measurments with SNR>3 SNRi_cut = 50. #I measurments with SNR>30, which implies an uncertainty in P of 4.7%. step_vec = 1 #plot all vectors in the array. if step_vec = 2, then every other vector will be plotted @@ -173,16 +173,13 @@ def main(): # Estimate error from data background, estimated from sub-image of desired sub_shape. background = None - if px_scale.lower() not in ['full','integrate']: - data_array, error_array, headers, background = proj_red.get_error(data_array, headers, error_array, sub_type=error_sub_type, subtract_error=subtract_error, display=display_error, savename=figname+"_errors", plots_folder=plots_folder, return_background=True) + data_array, error_array, headers, background = proj_red.get_error(data_array, headers, error_array, sub_type=error_sub_type, subtract_error=subtract_error, display=display_error, savename=figname+"_errors", plots_folder=plots_folder, return_background=True) # Align and rescale images with oversampling. data_array, error_array, headers, data_mask = proj_red.align_data(data_array, headers, error_array=error_array, background=background, upsample_factor=10, ref_center=align_center, return_shifts=False) # Rebin data to desired pixel size. if rebin: - if px_scale.lower() in ['full','integrate']: - data_array, error_array, headers = proj_red.get_error(data_array, headers, error_array, data_mask, sub_shape=error_sub_shape, display=display_error, savename=figname+"_errors", plots_folder=plots_folder) data_array, error_array, headers, Dxy, data_mask = proj_red.rebin_array(data_array, error_array, headers, pxsize=pxsize, scale=px_scale, operation=rebin_operation, data_mask=data_mask) # Rotate data to have North up @@ -229,6 +226,9 @@ def main(): stokescrop.writeto(data_folder+figname+figtype+".fits") Stokes_test, data_mask = stokescrop.hdul_crop, stokescrop.data_mask + 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(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.)) # 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, savename=figname+figtype, plots_folder=plots_folder) diff --git a/src/lib/background.py b/src/lib/background.py index 436ffdc..6491e03 100644 --- a/src/lib/background.py +++ b/src/lib/background.py @@ -89,10 +89,10 @@ def display_bkg(data, background, std_bkg, headers, histograms=None, binning=Non filt = headers[0]['filtnam1'] #plots im = ax2.imshow(data0, norm=LogNorm(data0[data0>0.].mean()/10.,data0.max()), origin='lower', cmap='gray') - bkg_im = ax2.imshow(bkg_data0, origin='lower', cmap='Reds', alpha=0.7) + bkg_im = ax2.imshow(bkg_data0, origin='lower', cmap='Reds', alpha=0.5) if not(rectangle is None): x, y, width, height, angle, color = rectangle[0] - ax2.add_patch(Rectangle((x, y),width,height,edgecolor=color,fill=False)) + ax2.add_patch(Rectangle((x, y),width,height,edgecolor=color,fill=False,lw=2)) ax2.annotate(instr+":"+rootname, color='white', fontsize=10, xy=(0.02, 0.95), xycoords='axes fraction') ax2.annotate(filt, color='white', fontsize=14, xy=(0.02, 0.02), xycoords='axes fraction') ax2.annotate(str(exptime)+" s", color='white', fontsize=10, xy=(0.80, 0.02), xycoords='axes fraction') diff --git a/src/lib/fits.py b/src/lib/fits.py index f45c18e..6d6146f 100755 --- a/src/lib/fits.py +++ b/src/lib/fits.py @@ -169,6 +169,7 @@ def save_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, P, debiased_P, s_P, header['datatype'] = ('I_stokes', 'type of data stored in the HDU') I_stokes[(1-data_mask).astype(bool)] = 0. primary_hdu = fits.PrimaryHDU(data=I_stokes, header=header) + primary_hdu.name = 'I_stokes' hdul.append(primary_hdu) #Add Q, U, Stokes_cov, P, s_P, PA, s_PA to the HDUList @@ -183,6 +184,7 @@ def save_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, P, debiased_P, s_P, if not name == 'IQU_cov_matrix': data[(1-data_mask).astype(bool)] = 0. hdu = fits.ImageHDU(data=data,header=hdu_header) + hdu.name = name hdul.append(hdu) #Save fits file to designated filepath diff --git a/src/lib/plots.py b/src/lib/plots.py index edfaf47..893a4f3 100755 --- a/src/lib/plots.py +++ b/src/lib/plots.py @@ -67,18 +67,23 @@ def princ_angle(ang): return A[0] -def sci_not(v,err,rnd=1): +def sci_not(v,err,rnd=1,out=str): """ Return the scientifque error notation as a string. """ power = - int(('%E' % v)[-3:])+1 - output = r"({0}".format(round(v*10**power,rnd)) + output = [r"({0}".format(round(v*10**power,rnd)),round(v*10**power,rnd)] if type(err) == list: for error in err: - output += r" $\pm$ {0}".format(round(error*10**power,rnd)) + output[0] += r" $\pm$ {0}".format(round(error*10**power,rnd)) + output.append(round(error*10**power,rnd)) else: - output += r" $\pm$ {0}".format(round(err*10**power,rnd)) - return output+r")e{0}".format(-power) + output[0] += r" $\pm$ {0}".format(round(err*10**power,rnd)) + output.append(round(err*10**power,rnd)) + if out==str: + return output[0]+r")e{0}".format(-power) + else: + return *output[1:],-power def plot_obs(data_array, headers, shape=None, vmin=0., vmax=6., rectangle=None, @@ -309,7 +314,7 @@ 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' - vmin, vmax = np.max(np.sqrt(stk_cov.data[0,0][mask])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) + vmin, vmax = 3.*np.mean(np.sqrt(stk_cov.data[0,0][mask])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) 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) @@ -320,7 +325,7 @@ 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.) - vmin, vmax = np.max(np.sqrt(stk_cov.data[0,0][mask])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) + vmin, vmax = 3.*np.mean(np.sqrt(stk_cov.data[0,0][mask])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) 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) @@ -350,13 +355,13 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c elif display.lower() in ['s_i','i_err']: # Display intensity error map display='s_i' - vmin, vmax = 0., np.max(np.sqrt(stk_cov.data[0,0][stk_cov.data[0,0] > 0.])*convert_flux) + vmin, vmax = np.min(np.sqrt(stk_cov.data[0,0][stk_cov.data[0,0] > 0.])*convert_flux), np.max(np.sqrt(stk_cov.data[0,0][stk_cov.data[0,0] > 0.])*convert_flux) im = ax.imshow(np.sqrt(stk_cov.data[0,0])*convert_flux, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) cbar = plt.colorbar(im, cax=cbar_ax, label=r"$\sigma_I$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") elif display.lower() in ['snr','snri']: # Display I_stokes signal-to-noise map display='snri' - vmin, vmax = 0., np.max(SNRi[SNRi > 0.]) + vmin, vmax = 0., np.max(SNRi[np.isfinite(SNRi)]) im = ax.imshow(SNRi, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) cbar = plt.colorbar(im, cax=cbar_ax, label=r"$I_{Stokes}/\sigma_{I}$") levelsSNRi = np.linspace(SNRi_cut, vmax*0.99, 10) @@ -366,7 +371,7 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c elif display.lower() in ['snrp']: # Display polarization degree signal-to-noise map display='snrp' - vmin, vmax = SNRp_cut, np.max(SNRp[SNRp > 0.]) + vmin, vmax = 0., np.max(SNRp[np.isfinite(SNRp)]) im = ax.imshow(SNRp, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) cbar = plt.colorbar(im, cax=cbar_ax, label=r"$P/\sigma_{P}$") levelsSNRp = np.linspace(SNRp_cut, vmax*0.99, 10) @@ -375,7 +380,7 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c #ax.clabel(cont,inline=True,fontsize=6) else: # Defaults to intensity map - vmin, vmax = np.min(stkI.data[SNRi > SNRi_cut]*convert_flux)/5., np.max(stkI.data[SNRi > SNRi_cut]*convert_flux) + vmin, vmax = 3.*np.mean(np.sqrt(stk_cov.data[0,0][mask])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) #im = ax.imshow(stkI.data*convert_flux, vmin=vmin, vmax=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$]") im = ax.imshow(stkI.data*convert_flux, norm=LogNorm(vmin,vmax), aspect='equal', cmap='inferno', alpha=1.) @@ -1545,6 +1550,7 @@ class pol_map(object): def submit_save(expression): ax_text_save.set(visible=False) if expression != '': + plt.rcParams.update({'font.size': 15}) save_fig = plt.figure(figsize=(15,15)) save_ax = save_fig.add_subplot(111, projection=self.wcs) self.ax_cosmetics(ax=save_ax) @@ -1560,6 +1566,7 @@ class pol_map(object): ax_snr_reset.set(visible=True) ax_save.set(visible=True) ax_dump.set(visible=True) + plt.rcParams.update({'font.size': 10}) self.fig.canvas.draw_idle() text_save.on_submit(submit_save) @@ -1736,17 +1743,17 @@ class pol_map(object): self.display_selection = "total_flux" if self.display_selection.lower() in ['total_flux']: self.data = self.I*self.convert_flux - vmin, vmax = np.max(np.sqrt(self.IQU_cov[0,0][self.cut])*self.convert_flux), np.max(self.data[self.data > 0.]) + vmin, vmax = .5*np.mean(np.sqrt(self.IQU_cov[0,0][self.IQU_cov[0,0] < 3.*self.I])*self.convert_flux), np.max(self.data[self.data > 0.]) 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 = np.max(np.sqrt(self.IQU_cov[0,0][self.cut])*self.convert_flux), np.max(self.I[self.data > 0.]*self.convert_flux) + vmin, vmax = .5*np.mean(np.sqrt(self.IQU_cov[0,0][self.IQU_cov[0,0] < 3.*self.I])*self.convert_flux), np.max(self.I[self.data > 0.]*self.convert_flux) 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']: self.data = self.P*100. - vmin, vmax = 0., 100. #np.max(self.data[self.data > 0.]) + vmin, vmax = 0., np.max(self.data[self.P > self.s_P]) label = r"$P$ [%]" elif self.display_selection.lower() in ['pol_ang']: self.data = princ_angle(self.PA)