general improvements to pipeline

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
+    pxsize = 0.10
-    px_scale = 'pixel'         #pixel, arcsec or full
+    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_function = 'weighted_gaussian'  #gaussian_after, weighted_gaussian_after, gaussian, weighted_gaussian or combine
-    smoothing_FWHM = None           #If None, no smoothing is done
+    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,7 +173,6 @@ 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)
 
     # Align and rescale images with oversampling.
@@ -181,8 +180,6 @@ def main():
 
     # 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)

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')

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

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))
+        output[0] += r" $\pm$ {0}".format(round(err*10**power,rnd))
-    return output+r")e{0}".format(-power)
+        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)