correction on error propagation in compute_Stokes

src/FOC_reduction.py

@@ -103,7 +103,7 @@ def main():
     rebin = True
     if rebin:
         pxsize = 0.10
-        px_scale = 'arcsec'         #pixel or arcsec
+        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
@@ -117,10 +117,11 @@ def main():
     rotate_data = False              #rotation to North convention can give erroneous results
     # Polarization map output
     figname = 'IC5063_FOC'         #target/intrument name
-    figtype = '_combine_FWHM020'    #additionnal informations
+    figtype = '_combine_FWHM020_pol'    #additionnal informations
-    SNRp_cut = 7.    #P measurments with SNR>3
+    SNRp_cut = 3.    #P measurments with SNR>3
-    SNRi_cut = 180.   #I measurments with SNR>30, which implies an uncertainty in P of 4.7%.
+    SNRi_cut = 70.   #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
+    step_vec = 0    #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
 
     ##### Pipeline start
     ## Step 1:
@@ -150,12 +151,17 @@ def main():
         if (data < 0.).any():
             print("ETAPE 4 : ", data)
     # Align and rescale images with oversampling.
+    data_mask = np.zeros(data_array.shape[1:]).astype(bool)
+    if px_scale.lower() not in ['full','integrate']:
         data_array, error_array, headers, data_mask = proj_red.align_data(data_array, headers, error_array, upsample_factor=int(Dxy.min()), ref_center=align_center, return_shifts=False)
     for data in data_array:
         if (data < 0.).any():
             print("ETAPE 5 : ", data)
 
+    if px_scale.lower() not in ['full','integrate']:
         vertex = image_hull((1.-data_mask),step=5,null_val=0.,inside=True)
+    else:
+        vertex = np.array([0.,0.,data_array.shape[2],data_array.shape[2]])
     shape = np.array([vertex[1]-vertex[0],vertex[3]-vertex[2]])
     rectangle = [vertex[2], vertex[0], shape[1], shape[0], 0., 'w']
 
@@ -206,6 +212,7 @@ def main():
     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, figname+figtype, data_folder=data_folder, return_hdul=True)
 
     # Plot polarization map (Background is either total Flux, Polarization degree or Polarization degree error).
+    if px_scale.lower() not in ['full','integrate']:
         proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, rectangle=None, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype, plots_folder=plots_folder, display=None)
         proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, rectangle=None, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype+"_P_flux", plots_folder=plots_folder, display='Pol_Flux')
         proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, rectangle=None, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype+"_P", plots_folder=plots_folder, display='Pol_deg')
@@ -213,6 +220,8 @@ def main():
         proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, rectangle=None, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype+"_P_err", plots_folder=plots_folder, display='Pol_deg_err')
         proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, rectangle=None, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype+"_SNRi", plots_folder=plots_folder, display='SNRi')
         proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, rectangle=None, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype+"_SNRp", plots_folder=plots_folder, display='SNRp')
+    else:
+        proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, rectangle=None, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype, plots_folder=plots_folder, display='default')
 
     return 0
 

src/lib/plots.py

@@ -328,7 +328,7 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
         data_mask = np.ones(stkI.shape).astype(bool)
 
     #Plot Stokes parameters map
-    if display is None:
+    if display is None or display.lower() == 'default':
         plot_Stokes(Stokes, savename=savename, plots_folder=plots_folder)
 
     #Compute SNR and apply cuts
@@ -406,19 +406,20 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
         cont = ax.contour(SNRp, levels=levelsP, colors='grey', linewidths=0.5)
     else:
         # Defaults to intensity map
-        vmin, vmax = 0., np.max(stkI.data[stkI.data > 0.]*convert_flux)
+        vmin, vmax = 0., np.max(stkI.data[stkI.data > 0.]*convert_flux*2.)
         im = ax.imshow(stkI.data*convert_flux, vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.)
         cbar = plt.colorbar(im, cax=cbar_ax, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]")
-        levelsI = np.linspace(SNRi_cut, SNRi.max(), 10)
-        cont = ax.contour(SNRi, levels=levelsI, colors='grey', linewidths=0.5)
 
+    if (display is None) or not(display.lower() in ['default']):
         fontprops = fm.FontProperties(size=16)
         px_size = wcs.wcs.get_cdelt()[0]*3600.
         px_sc = AnchoredSizeBar(ax.transData, 1./px_size, '1 arcsec', 3, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color='w', fontproperties=fontprops)
         ax.add_artist(px_sc)
 
-    #pol.data[np.isfinite(pol.data)] = 1./2.
+        if step_vec == 0:
-    X, Y = np.meshgrid(np.linspace(0,stkI.data.shape[0],stkI.data.shape[0]), np.linspace(0,stkI.data.shape[1],stkI.data.shape[1]))
+            pol.data[np.isfinite(pol.data)] = 1./2.
+            step_vec = 1
+        X, Y = np.meshgrid(np.linspace(0,stkI.data.shape[0],stkI.data.shape[0])-0.5, np.linspace(0,stkI.data.shape[1],stkI.data.shape[1])-0.5)
         U, V = pol.data*np.cos(np.pi/2.+pang.data*np.pi/180.), pol.data*np.sin(np.pi/2.+pang.data*np.pi/180.)
         Q = 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=0.5,scale_units='xy',pivot='mid',headwidth=0.,headlength=0.,headaxislength=0.,width=0.1,color='w')
         pol_sc = AnchoredSizeBar(ax.transData, 2., r"$P$= 100 %", 4, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color='w', fontproperties=fontprops)

src/lib/reduction.py

@@ -568,11 +568,12 @@ def rebin_array(data_array, error_array, headers, pxsize, scale,
                 Dxy = np.array([pxsize,]*2)
             elif scale.lower() in ['arcsec','arcseconds']:
                 Dxy = np.floor(pxsize/w.wcs.cdelt/3600.).astype(int)
+            elif scale.lower() in ['full','integrate']:
+                Dxy = np.floor(image.shape).astype(int)
             else:
                 raise ValueError("'{0:s}' invalid scale for binning.".format(scale))
 
             if (Dxy <= 1.).any():
-                print(Dxy, pxsize, w.wcs.cdelt*3600.)
                 raise ValueError("Requested pixel size is below resolution.")
             new_shape = (image.shape//Dxy).astype(int)
 
@@ -1049,7 +1050,7 @@ def compute_Stokes(data_array, error_array, data_mask, headers,
     # Routine for the FOC instrument
     if instr == 'FOC':
         # Get image from each polarizer and covariance matrix
-        pol_array, pol_cov = polarizer_avg(data_array, error_array, data_mask,
+        pol_array, pol_cov_m = polarizer_avg(data_array, error_array, data_mask,
                 headers, FWHM=FWHM, scale=scale, smoothing=smoothing)
         pol0, pol60, pol120 = pol_array
 
@@ -1088,6 +1089,7 @@ def compute_Stokes(data_array, error_array, data_mask, headers,
         # Uncertainties on the orientation of the polarizers' axes taken to be 3deg (see Nota et. al 1996, p36; Robinson & Thomson 1995)
         sigma_theta = np.array([3.*np.pi/180., 3.*np.pi/180., 3.*np.pi/180.])
         pol_flux = 2.*np.array([pol0/transmit[0], pol60/transmit[1], pol120/transmit[2]])
+        pol_cov = np.array([[4.*pol_cov_m[i,j]/(transmit[i]*transmit[j]) for j in range(pol_cov_m.shape[1])] for i in range(pol_cov_m.shape[0])])
 
         # Normalization parameter for Stokes parameters computation
         A = pol_eff[1]*pol_eff[2]*np.sin(-2.*theta[1]+2.*theta[2]) \
@@ -1141,9 +1143,9 @@ def compute_Stokes(data_array, error_array, data_mask, headers,
         s_U2_axis = (dU_dtheta1**2*sigma_theta[0]**2 + dU_dtheta2**2*sigma_theta[1]**2 + dU_dtheta3**2*sigma_theta[2]**2)
 
         # Add quadratically the uncertainty to the Stokes covariance matrix ## THIS IS WHERE THE PROBLEMATIC UNCERTAINTY IS ADDED TO THE PIPELINE
-        #Stokes_cov[0,0] += s_I2_axis
+        Stokes_cov[0,0] += s_I2_axis
-        #Stokes_cov[1,1] += s_Q2_axis
+        Stokes_cov[1,1] += s_Q2_axis
-        #Stokes_cov[2,2] += s_U2_axis
+        Stokes_cov[2,2] += s_U2_axis
 
         if not(FWHM is None) and (smoothing.lower() in ['gaussian_after','gauss_after']):
             Stokes_array = np.array([I_stokes, Q_stokes, U_stokes])