implement axis error in compute_pol, redo some plots

src/FOC_reduction.py

@@ -17,11 +17,11 @@ from lib.convex_hull import image_hull
 def main():
     ##### User inputs
     ## Input and output locations
-#    globals()['data_folder'] = "../data/NGC1068_x274020/"
+    globals()['data_folder'] = "../data/NGC1068_x274020/"
-#    infiles = ['x274020at.c0f.fits','x274020bt.c0f.fits','x274020ct.c0f.fits',
+    infiles = ['x274020at.c0f.fits','x274020bt.c0f.fits','x274020ct.c0f.fits',
-#            'x274020dt.c0f.fits','x274020et.c0f.fits','x274020ft.c0f.fits',
+            'x274020dt.c0f.fits','x274020et.c0f.fits','x274020ft.c0f.fits',
-#            'x274020gt.c0f.fits','x274020ht.c0f.fits','x274020it.c0f.fits']
+            'x274020gt.c0f.fits','x274020ht.c0f.fits','x274020it.c0f.fits']
-#    globals()['plots_folder'] = "../plots/NGC1068_x274020/"
+    globals()['plots_folder'] = "../plots/NGC1068_x274020/"
 
 #    globals()['data_folder'] = "../data/NGC1068_x14w010/"
 #    infiles = ['x14w0101t_c0f.fits','x14w0102t_c0f.fits','x14w0103t_c0f.fits',
@@ -60,9 +60,9 @@ def main():
 #            'x3995202r_c0f.fits','x3995206r_c0f.fits']
 #    globals()['plots_folder'] = "../plots/PG1630+377_x39510/"
 
-    globals()['data_folder'] = "../data/IC5063_x3nl030/"
+#    globals()['data_folder'] = "../data/IC5063_x3nl030/"
-    infiles = ['x3nl0301r_c0f.fits','x3nl0302r_c0f.fits','x3nl0303r_c0f.fits']
+#    infiles = ['x3nl0301r_c0f.fits','x3nl0302r_c0f.fits','x3nl0303r_c0f.fits']
-    globals()['plots_folder'] = "../plots/IC5063_x3nl030/"
+#    globals()['plots_folder'] = "../plots/IC5063_x3nl030/"
 
 #    globals()['data_folder'] = "../data/MKN3_x3nl010/"
 #    infiles = ['x3nl0101r_c0f.fits','x3nl0102r_c0f.fits','x3nl0103r_c0f.fits']
@@ -77,10 +77,6 @@ def main():
 #    infiles = ['x3nl0201r_c0f.fits','x3nl0202r_c0f.fits','x3nl0203r_c0f.fits']
 #    globals()['plots_folder'] = "../plots/MKN78_x3nl020/"
 
-#    globals()['data_folder'] = "../data/PictorA_x25d040/"
-#    infiles = ['x25d0401t_c0f.fits','x25d0402t_c0f.fits','x25d0403t_c0f.fits']
-#    globals()['plots_folder'] = "../plots/PictorA_x25d040/"
-
 #    globals()['data_folder'] = "../data/3C273_x0u20/"
 #    infiles = ['x0u20101t_c0f.fits','x0u20102t_c0f.fits','x0u20103t_c0f.fits','x0u20104t_c0f.fits','x0u20105t_c0f.fits','x0u20106t_c0f.fits','x0u20201t_c0f.fits','x0u20202t_c0f.fits','x0u20203t_c0f.fits','x0u20204t_c0f.fits','x0u20205t_c0f.fits','x0u20206t_c0f.fits','x0u20301t_c0f.fits','x0u20302t_c0f.fits','x0u20303t_c0f.fits','x0u20304t_c0f.fits','x0u20305t_c0f.fits','x0u20306t_c0f.fits']
 #    globals()['plots_folder'] = "../plots/3C273_x0u20/"
@@ -116,11 +112,11 @@ def main():
     rotate_stokes = True           #rotation to North convention can give erroneous results
     rotate_data = False              #rotation to North convention can give erroneous results
     # Polarization map output
-    figname = 'IC5063_FOC'         #target/intrument name
+    figname = 'NGC1068_FOC'         #target/intrument name
-    figtype = '_combine_FWHM020_pol'    #additionnal informations
+    figtype = '_combine_FWHM020_wae'    #additionnal informations
     SNRp_cut = 3.    #P measurments with SNR>3
-    SNRi_cut = 70.   #I measurments with SNR>30, which implies an uncertainty in P of 4.7%.
+    SNRi_cut = 30.   #I measurments with SNR>30, which implies an uncertainty in P of 4.7%.
-    step_vec = 0    #plot all vectors in the array. if step_vec = 2, then every other vector will be plotted
+    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
 
     ##### Pipeline start
@@ -163,7 +159,7 @@ def main():
     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']
+    rectangle = [vertex[2], vertex[0], shape[1], shape[0], 0., 'g']
 
     # Rotate data to have North up
     ref_header = deepcopy(headers[0])
@@ -186,7 +182,7 @@ def main():
     # FWHM of FOC have been estimated at about 0.03" across 1500-5000 Angstrom band, which is about 2 detector pixels wide
     # see Jedrzejewski, R.; Nota, A.; Hack, W. J., A Comparison Between FOC and WFPC2
     # Bibcode : 1995chst.conf...10J
-    I_stokes, Q_stokes, U_stokes, Stokes_cov = proj_red.compute_Stokes(data_array, error_array, data_mask, headers, FWHM=smoothing_FWHM, scale=smoothing_scale, smoothing=smoothing_function)
+    I_stokes, Q_stokes, U_stokes, Stokes_cov, dP_dtheta, dPA_dtheta = proj_red.compute_Stokes(data_array, error_array, data_mask, headers, FWHM=smoothing_FWHM, scale=smoothing_scale, smoothing=smoothing_function)
 
     ## Step 3:
     # Rotate images to have North up
@@ -199,7 +195,7 @@ def main():
         rectangle[4] = alpha
         I_stokes, Q_stokes, U_stokes, Stokes_cov, headers, data_mask = proj_red.rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers, -ref_header['orientat'], SNRi_cut=None)
     # Compute polarimetric parameters (polarization degree and angle).
-    P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P = proj_red.compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers)
+    P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P = proj_red.compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, dP_dtheta, dPA_dtheta, headers)
 
     ## Step 4:
     # crop to desired region of interest (roi)

src/lib/plots.py

@@ -93,16 +93,16 @@ def plot_obs(data_array, headers, shape=None, vmin=0., vmax=6., rectangle=None,
         exptime = headers[i]['exptime']
         filt = headers[i]['filtnam1']
         #plots
-        im = ax.imshow(data, vmin=vmin, vmax=vmax, origin='lower')
+        im = ax.imshow(data, vmin=vmin, vmax=vmax, origin='lower', cmap='gray')
         if not(rectangle is None):
             x, y, width, height, angle, color = rectangle[i]
             ax.add_patch(Rectangle((x, y), width, height, angle=angle,
                 edgecolor=color, fill=False))
         #position of centroid
-        ax.plot([data.shape[1]/2, data.shape[1]/2], [0,data.shape[0]-1], lw=1,
+        ax.plot([data.shape[1]/2, data.shape[1]/2], [0,data.shape[0]-1], '--', lw=1,
-                color='black')
+                color='grey', alpha=0.5)
-        ax.plot([0,data.shape[1]-1], [data.shape[1]/2, data.shape[1]/2], lw=1,
+        ax.plot([0,data.shape[1]-1], [data.shape[1]/2, data.shape[1]/2], '--', lw=1,
-                color='black')
+                color='grey', alpha=0.5)
         ax.annotate(instr+":"+rootname,color='white',fontsize=5,xy=(0.02, 0.95),
                 xycoords='axes fraction')
         ax.annotate(filt,color='white',fontsize=10,xy=(0.02, 0.02),
@@ -110,12 +110,12 @@ def plot_obs(data_array, headers, shape=None, vmin=0., vmax=6., rectangle=None,
         ax.annotate(exptime,color='white',fontsize=5,xy=(0.80, 0.02),
                 xycoords='axes fraction')
 
-    fig.subplots_adjust(hspace=0, wspace=0, right=0.85)
+    fig.subplots_adjust(hspace=0.01, wspace=0.01, right=0.85)
     cbar_ax = fig.add_axes([0.9, 0.12, 0.02, 0.75])
-    fig.colorbar(im, cax=cbar_ax)
+    fig.colorbar(im, cax=cbar_ax, label=r'$Counts \cdot s^{-1}$')
 
     if not (savename is None):
-        fig.suptitle(savename)
+        #fig.suptitle(savename)
         fig.savefig(plots_folder+savename+".png",bbox_inches='tight')
     plt.show()
     return 0
@@ -149,121 +149,27 @@ def plot_Stokes(Stokes, savename=None, plots_folder=""):
     fig = plt.figure(figsize=(30,10))
 
     ax = fig.add_subplot(131, projection=wcs)
-    im = ax.imshow(stkI, origin='lower')
+    im = ax.imshow(stkI, origin='lower', cmap='inferno')
     plt.colorbar(im)
     ax.set(xlabel="RA", ylabel="DEC", title=r"$I_{stokes}$")
 
     ax = fig.add_subplot(132, projection=wcs)
-    im = ax.imshow(stkQ, origin='lower')
+    im = ax.imshow(stkQ, origin='lower', cmap='inferno')
     plt.colorbar(im)
     ax.set(xlabel="RA", ylabel="DEC", title=r"$Q_{stokes}$")
 
     ax = fig.add_subplot(133, projection=wcs)
-    im = ax.imshow(stkU, origin='lower')
+    im = ax.imshow(stkU, origin='lower', cmap='inferno')
     plt.colorbar(im)
     ax.set(xlabel="RA", ylabel="DEC", title=r"$U_{stokes}$")
 
     if not (savename is None):
-        fig.suptitle(savename+"_IQU")
+        #fig.suptitle(savename+"_IQU")
         fig.savefig(plots_folder+savename+"_IQU.png",bbox_inches='tight')
     plt.show()
     return 0
 
 
-class crop_map(object):
-    """
-    Class to interactively crop I_stokes map to desired Region of Interest
-    """
-    def __init__(self, Stokes, data_mask, SNRp_cut=3., SNRi_cut=30.):
-        #Get data
-        stkI = Stokes[np.argmax([Stokes[i].header['datatype']=='I_stokes' for i in range(len(Stokes))])]
-        stk_cov = Stokes[np.argmax([Stokes[i].header['datatype']=='IQU_cov_matrix' for i in range(len(Stokes))])]
-        pol = Stokes[np.argmax([Stokes[i].header['datatype']=='Pol_deg_debiased' for i in range(len(Stokes))])]
-        pol_err = Stokes[np.argmax([Stokes[i].header['datatype']=='Pol_deg_err' for i in range(len(Stokes))])]
-        self.Stokes = Stokes
-        self.data_mask = data_mask
-
-        wcs = WCS(Stokes[0]).deepcopy()
-
-        #Compute SNR and apply cuts
-        pol.data[pol.data == 0.] = np.nan
-        SNRp = pol.data/pol_err.data
-        SNRp[np.isnan(SNRp)] = 0.
-        pol.data[SNRp < SNRp_cut] = np.nan
-        SNRi = stkI.data/np.sqrt(stk_cov.data[0,0])
-        SNRi[np.isnan(SNRi)] = 0.
-        pol.data[SNRi < SNRi_cut] = np.nan
-
-        convert_flux = Stokes[0].header['photflam']
-
-        #Plot the map
-        plt.rcParams.update({'font.size': 16})
-        self.fig = plt.figure(figsize=(15,15))
-        self.ax = self.fig.add_subplot(111, projection=wcs)
-        self.ax.set_facecolor('k')
-        self.fig.subplots_adjust(hspace=0, wspace=0, right=0.9)
-        cbar_ax = self.fig.add_axes([0.95, 0.12, 0.01, 0.75])
-
-        self.extent = [-stkI.data.shape[1]/2.,stkI.data.shape[1]/2.,-stkI.data.shape[0]/2.,stkI.data.shape[0]/2.]
-        self.center = [stkI.data.shape[1]/2.,stkI.data.shape[0]/2.]
-
-        vmin, vmax = 0., np.max(stkI.data[stkI.data > 0.]*convert_flux)
-        im = self.ax.imshow(stkI.data*convert_flux,extent=self.extent, 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^{-1}$]")
-        levelsI = np.linspace(SNRi_cut, np.max(SNRi[SNRi > 0.]), 10)
-        cont = self.ax.contour(SNRi, extent=self.extent, levels=levelsI, colors='grey', linewidths=0.5)
-
-        fontprops = fm.FontProperties(size=16)
-        px_size = wcs.wcs.get_cdelt()[0]
-        px_sc = AnchoredSizeBar(self.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)
-        self.ax.add_artist(px_sc)
-
-        self.ax.set_title(
-        "Click and drag to crop to desired Region of Interest.\n"
-        "Press 'c' to toggle the selector on and off")
-
-    def onselect_crop(self, eclick, erelease) -> None:
-        # Obtain (xmin, xmax, ymin, ymax) values
-        self.RSextent = self.rect_selector.extents
-        self.RScenter = [self.center[i]+self.rect_selector.center[i] for i in range(2)]
-
-        # Zoom to selection
-        print("CROP TO : ",self.RSextent)
-        print("CENTER : ",self.RScenter)
-        #self.ax.set_xlim(self.extent[0], self.extent[1])
-        #self.ax.set_ylim(self.extent[2], self.extent[3])
-
-    def run(self) -> None:
-        self.rect_selector = RectangleSelector(self.ax, self.onselect_crop,
-                drawtype='box', button=[1], interactive=True)
-        #self.fig.canvas.mpl_connect('key_press_event', self.toggle_selector)
-        plt.show()
-
-    def crop(self):
-        Stokes_crop = copy.deepcopy(self.Stokes)
-        # Data sets to crop
-        stkI = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='I_stokes' for i in range(len(Stokes_crop))])]
-        stkQ = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='Q_stokes' for i in range(len(Stokes_crop))])]
-        stkU = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='U_stokes' for i in range(len(Stokes_crop))])]
-        stk_cov = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='IQU_cov_matrix' for i in range(len(Stokes_crop))])]
-        pol = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='Pol_deg_debiased' for i in range(len(Stokes_crop))])]
-        pol_err = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='Pol_deg_err' for i in range(len(Stokes_crop))])]
-        pang = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='Pol_ang' for i in range(len(Stokes_crop))])]
-        pang_err = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='Pol_ang_err' for i in range(len(Stokes_crop))])]
-        # Crop datasets
-        vertex = [int(self.RScenter[0]+self.RSextent[0]), int(self.RScenter[0]+self.RSextent[1]), int(self.RScenter[1]+self.RSextent[2]), int(self.RScenter[1]+self.RSextent[3])]
-        for dataset in [stkI, stkQ, stkU, pol, pol_err, pang, pang_err]:
-            dataset.data = dataset.data[vertex[2]:vertex[3], vertex[0]:vertex[1]]
-        data_mask = self.data_mask[vertex[2]:vertex[3], vertex[0]:vertex[1]]
-        new_stk_cov = np.zeros((3,3,stkI.data.shape[0],stkI.data.shape[1]))
-        for i in range(3):
-            for j in range(3):
-                new_stk_cov[i,j] = stk_cov.data[i,j][vertex[2]:vertex[3], vertex[0]:vertex[1]]
-        stk_cov.data = new_stk_cov
-
-        return Stokes_crop, data_mask
-
-
 def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_cut=30.,
         step_vec=1, savename=None, plots_folder="", display=None):
     """
@@ -366,7 +272,9 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
         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^{-1}$]")
         levelsI = np.linspace(SNRi_cut, np.max(SNRi[SNRi > 0.]), 10)
+        print("SNRi contour levels : ", levelsI)
         cont = ax.contour(SNRi, levels=levelsI, colors='grey', linewidths=0.5)
+        #ax.clabel(cont,inline=True,fontsize=6)
     elif display.lower() in ['pol_flux']:
         # Display polarisation flux
         pf_mask = (stkI.data > 0.) * (pol.data > 0.)
@@ -674,3 +582,98 @@ class overplot_maps(align_maps):
         if self.aligned:
             self.overplot(other_levels=levels, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, savename=savename)
             plt.show(block=True)
+
+
+class crop_map(object):
+    """
+    Class to interactively crop I_stokes map to desired Region of Interest
+    """
+    def __init__(self, Stokes, data_mask, SNRp_cut=3., SNRi_cut=30.):
+        #Get data
+        stkI = Stokes[np.argmax([Stokes[i].header['datatype']=='I_stokes' for i in range(len(Stokes))])]
+        stk_cov = Stokes[np.argmax([Stokes[i].header['datatype']=='IQU_cov_matrix' for i in range(len(Stokes))])]
+        pol = Stokes[np.argmax([Stokes[i].header['datatype']=='Pol_deg_debiased' for i in range(len(Stokes))])]
+        pol_err = Stokes[np.argmax([Stokes[i].header['datatype']=='Pol_deg_err' for i in range(len(Stokes))])]
+        self.Stokes = Stokes
+        self.data_mask = data_mask
+
+        wcs = WCS(Stokes[0]).deepcopy()
+
+        #Compute SNR and apply cuts
+        pol.data[pol.data == 0.] = np.nan
+        SNRp = pol.data/pol_err.data
+        SNRp[np.isnan(SNRp)] = 0.
+        pol.data[SNRp < SNRp_cut] = np.nan
+        SNRi = stkI.data/np.sqrt(stk_cov.data[0,0])
+        SNRi[np.isnan(SNRi)] = 0.
+        pol.data[SNRi < SNRi_cut] = np.nan
+
+        convert_flux = Stokes[0].header['photflam']
+
+        #Plot the map
+        plt.rcParams.update({'font.size': 16})
+        self.fig = plt.figure(figsize=(15,15))
+        self.ax = self.fig.add_subplot(111, projection=wcs)
+        self.ax.set_facecolor('k')
+        self.fig.subplots_adjust(hspace=0, wspace=0, right=0.9)
+        cbar_ax = self.fig.add_axes([0.95, 0.12, 0.01, 0.75])
+
+        self.extent = [-stkI.data.shape[1]/2.,stkI.data.shape[1]/2.,-stkI.data.shape[0]/2.,stkI.data.shape[0]/2.]
+        self.center = [stkI.data.shape[1]/2.,stkI.data.shape[0]/2.]
+
+        vmin, vmax = 0., np.max(stkI.data[stkI.data > 0.]*convert_flux)
+        im = self.ax.imshow(stkI.data*convert_flux,extent=self.extent, 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^{-1}$]")
+        levelsI = np.linspace(SNRi_cut, np.max(SNRi[SNRi > 0.]), 10)
+        cont = self.ax.contour(SNRi, extent=self.extent, levels=levelsI, colors='grey', linewidths=0.5)
+
+        fontprops = fm.FontProperties(size=16)
+        px_size = wcs.wcs.get_cdelt()[0]
+        px_sc = AnchoredSizeBar(self.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)
+        self.ax.add_artist(px_sc)
+
+        self.ax.set_title(
+        "Click and drag to crop to desired Region of Interest.\n"
+        "Press 'c' to toggle the selector on and off")
+
+    def onselect_crop(self, eclick, erelease) -> None:
+        # Obtain (xmin, xmax, ymin, ymax) values
+        self.RSextent = self.rect_selector.extents
+        self.RScenter = [self.center[i]+self.rect_selector.center[i] for i in range(2)]
+
+        # Zoom to selection
+        print("CROP TO : ",self.RSextent)
+        print("CENTER : ",self.RScenter)
+        #self.ax.set_xlim(self.extent[0], self.extent[1])
+        #self.ax.set_ylim(self.extent[2], self.extent[3])
+
+    def run(self) -> None:
+        self.rect_selector = RectangleSelector(self.ax, self.onselect_crop,
+                drawtype='box', button=[1], interactive=True)
+        #self.fig.canvas.mpl_connect('key_press_event', self.toggle_selector)
+        plt.show()
+
+    def crop(self):
+        Stokes_crop = copy.deepcopy(self.Stokes)
+        # Data sets to crop
+        stkI = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='I_stokes' for i in range(len(Stokes_crop))])]
+        stkQ = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='Q_stokes' for i in range(len(Stokes_crop))])]
+        stkU = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='U_stokes' for i in range(len(Stokes_crop))])]
+        stk_cov = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='IQU_cov_matrix' for i in range(len(Stokes_crop))])]
+        pol = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='Pol_deg_debiased' for i in range(len(Stokes_crop))])]
+        pol_err = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='Pol_deg_err' for i in range(len(Stokes_crop))])]
+        pang = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='Pol_ang' for i in range(len(Stokes_crop))])]
+        pang_err = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='Pol_ang_err' for i in range(len(Stokes_crop))])]
+        # Crop datasets
+        vertex = [int(self.RScenter[0]+self.RSextent[0]), int(self.RScenter[0]+self.RSextent[1]), int(self.RScenter[1]+self.RSextent[2]), int(self.RScenter[1]+self.RSextent[3])]
+        for dataset in [stkI, stkQ, stkU, pol, pol_err, pang, pang_err]:
+            dataset.data = dataset.data[vertex[2]:vertex[3], vertex[0]:vertex[1]]
+        data_mask = self.data_mask[vertex[2]:vertex[3], vertex[0]:vertex[1]]
+        new_stk_cov = np.zeros((3,3,stkI.data.shape[0],stkI.data.shape[1]))
+        for i in range(3):
+            for j in range(3):
+                new_stk_cov[i,j] = stk_cov.data[i,j][vertex[2]:vertex[3], vertex[0]:vertex[1]]
+        stk_cov.data = new_stk_cov
+
+        return Stokes_crop, data_mask
+

src/lib/reduction.py

@@ -246,21 +246,21 @@ def crop_array(data_array, headers, error_array=None, step=5, null_val=None,
     new_shape = np.array([v_array[1]-v_array[0],v_array[3]-v_array[2]])
     rectangle = [v_array[2], v_array[0], new_shape[1], new_shape[0], 0., 'b']
     if display:
-        fig, ax = plt.subplots()
+        fig, ax = plt.subplots(figsize=(10,10))
         data = data_array[0]
         instr = headers[0]['instrume']
         rootname = headers[0]['rootname']
         exptime = headers[0]['exptime']
         filt = headers[0]['filtnam1']
         #plots
-        im = ax.imshow(data, vmin=data.min(), vmax=data.max(), origin='lower')
+        im = ax.imshow(data, vmin=data.min(), vmax=data.max(), origin='lower', cmap='gray')
         x, y, width, height, angle, color = rectangle
         ax.add_patch(Rectangle((x, y),width,height,edgecolor=color,fill=False))
         #position of centroid
-        ax.plot([data.shape[1]/2, data.shape[1]/2], [0,data.shape[0]-1], lw=1,
+        ax.plot([data.shape[1]/2, data.shape[1]/2], [0,data.shape[0]-1], '--', lw=1,
-                color='black')
+                color='grey', alpha=0.3)
-        ax.plot([0,data.shape[1]-1], [data.shape[1]/2, data.shape[1]/2], lw=1,
+        ax.plot([0,data.shape[1]-1], [data.shape[1]/2, data.shape[1]/2], '--', lw=1,
-                color='black')
+                color='grey', alpha=0.3)
         ax.annotate(instr+":"+rootname, color='white', fontsize=5,
                 xy=(0.02, 0.95), xycoords='axes fraction')
         ax.annotate(filt, color='white', fontsize=10, xy=(0.02, 0.02),
@@ -273,10 +273,10 @@ def crop_array(data_array, headers, error_array=None, step=5, null_val=None,
 
         fig.subplots_adjust(hspace=0, wspace=0, right=0.85)
         cbar_ax = fig.add_axes([0.9, 0.12, 0.02, 0.75])
-        fig.colorbar(im, cax=cbar_ax)
+        fig.colorbar(im, cax=cbar_ax, label=r'$Counts \cdot s^{-1}$')
 
         if not(savename is None):
-            fig.suptitle(savename+'_'+filt+'_crop_region')
+            #fig.suptitle(savename+'_'+filt+'_crop_region')
             fig.savefig(plots_folder+savename+'_'+filt+'_crop_region.png',
                     bbox_inches='tight')
             plot_obs(data_array, headers, vmin=data_array.min(),
@@ -486,7 +486,7 @@ def get_error(data_array, headers, sub_shape=(15,15), display=False,
         plt.legend()
 
         if not(savename is None):
-            fig.suptitle(savename+"_background_flux")
+            #fig.suptitle(savename+"_background_flux")
             fig.savefig(plots_folder+savename+"_background_flux.png",
                     bbox_inches='tight')
             vmin = np.min(np.log10(data[data > 0.]))
@@ -584,6 +584,10 @@ def rebin_array(data_array, error_array, headers, pxsize, scale,
             # Propagate error
             rms_image = np.sqrt(bin_ndarray(image**2, new_shape=new_shape,
                 operation='average'))
+            if operation.lower() in ["mean", "average", "avg"]:
+                new_error = 1./np.sqrt(Dxy[0]*Dxy[1])*bin_ndarray(error,
+                    new_shape=new_shape, operation='average')
+            else:
                 new_error = np.sqrt(Dxy[0]*Dxy[1])*bin_ndarray(error,
                     new_shape=new_shape, operation='average')
             rebinned_error.append(np.sqrt(rms_image**2 + new_error**2))
@@ -1083,11 +1087,12 @@ def compute_Stokes(data_array, error_array, data_mask, headers,
         pol_eff[1] = pol_efficiency['pol60']
         pol_eff[2] = pol_efficiency['pol120']
 
-        # Orientation and error for each polarizer  ## THIS IS WHERE WE IMPLEMENT THE ERROR THAT IS GOING WRONG
+        # Orientation and error for each polarizer
         # POL0 = 0deg, POL60 = 60deg, POL120=120deg
         theta = np.array([180.*np.pi/180., 60.*np.pi/180., 120.*np.pi/180.])
         # 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.])
+        fmax = np.finfo(np.float64).max
         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])])
 
@@ -1130,19 +1135,42 @@ def compute_Stokes(data_array, error_array, data_mask, headers,
         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_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_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_dtheta = np.array([dU_dtheta1, dU_dtheta2, dU_dtheta3])
+        
+        # Compute the derivative of each polarization component with respect to the polarizer orientation
+        mask = I_stokes > 0.
+        dP_dtheta1 = np.ones(I_stokes.shape)*fmax
+        dP_dtheta2 = np.ones(I_stokes.shape)*fmax
+        dP_dtheta3 = np.ones(I_stokes.shape)*fmax
+        dP_dtheta1[mask] = 1./(I_stokes[mask]*np.sqrt(Q_stokes[mask]**2 + U_stokes[mask]**2))*(Q_stokes[mask]*dQ_dtheta1[mask] + U_stokes[mask]*dQ_dtheta1[mask]) - 1./I_stokes[mask]**2 * dI_dtheta1[mask]
+        dP_dtheta2[mask] = 1./(I_stokes[mask]*np.sqrt(Q_stokes[mask]**2 + U_stokes[mask]**2))*(Q_stokes[mask]*dQ_dtheta2[mask] + U_stokes[mask]*dQ_dtheta2[mask]) - 1./I_stokes[mask]**2 * dI_dtheta2[mask]
+        dP_dtheta3[mask] = 1./(I_stokes[mask]*np.sqrt(Q_stokes[mask]**2 + U_stokes[mask]**2))*(Q_stokes[mask]*dQ_dtheta3[mask] + U_stokes[mask]*dQ_dtheta3[mask]) - 1./I_stokes[mask]**2 * dI_dtheta3[mask]
+        dP_dtheta = np.array([dP_dtheta1, dP_dtheta2, dP_dtheta3])
+
+        dPA_dtheta1 = np.ones(I_stokes.shape)*fmax
+        dPA_dtheta2 = np.ones(I_stokes.shape)*fmax
+        dPA_dtheta3 = np.ones(I_stokes.shape)*fmax
+        dPA_dtheta1[mask] = 90./np.pi*1./(U_stokes[mask]**2 + Q_stokes[mask]**2)*(Q_stokes[mask]*dU_dtheta1[mask] - U_stokes[mask]*dQ_dtheta1[mask])
+        dPA_dtheta2[mask] = 90./np.pi*1./(U_stokes[mask]**2 + Q_stokes[mask]**2)*(Q_stokes[mask]*dU_dtheta2[mask] - U_stokes[mask]*dQ_dtheta2[mask])
+        dPA_dtheta3[mask] = 90./np.pi*1./(U_stokes[mask]**2 + Q_stokes[mask]**2)*(Q_stokes[mask]*dU_dtheta3[mask] - U_stokes[mask]*dQ_dtheta3[mask])
+        dPA_dtheta = np.array([dPA_dtheta1, dPA_dtheta2, dPA_dtheta3])
 
         # Compute the uncertainty associated with the polarizers' orientation (see Kishimoto 1999)
-        s_I2_axis = (dI_dtheta1**2*sigma_theta[0]**2 + dI_dtheta2**2*sigma_theta[1]**2 + dI_dtheta3**2*sigma_theta[2]**2)
+        s_I2_axis = np.sum([dI_dtheta[i]**2 * sigma_theta[i]**2 for i in range(len(sigma_theta))],axis=0) #(dI_dtheta1**2*sigma_theta[0]**2 + dI_dtheta2**2*sigma_theta[1]**2 + dI_dtheta3**2*sigma_theta[2]**2)
-        s_Q2_axis = (dQ_dtheta1**2*sigma_theta[0]**2 + dQ_dtheta2**2*sigma_theta[1]**2 + dQ_dtheta3**2*sigma_theta[2]**2)
+        s_Q2_axis = np.sum([dQ_dtheta[i]**2 * sigma_theta[i]**2 for i in range(len(sigma_theta))],axis=0) #(dQ_dtheta1**2*sigma_theta[0]**2 + dQ_dtheta2**2*sigma_theta[1]**2 + dQ_dtheta3**2*sigma_theta[2]**2)
-        s_U2_axis = (dU_dtheta1**2*sigma_theta[0]**2 + dU_dtheta2**2*sigma_theta[1]**2 + dU_dtheta3**2*sigma_theta[2]**2)
+        s_U2_axis = np.sum([dU_dtheta[i]**2 * sigma_theta[i]**2 for i in range(len(sigma_theta))],axis=0) #(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
+        # Add quadratically the uncertainty to the Stokes covariance matrix
         Stokes_cov[0,0] += s_I2_axis
         Stokes_cov[1,1] += s_Q2_axis
         Stokes_cov[2,2] += s_U2_axis
@@ -1158,10 +1186,10 @@ def compute_Stokes(data_array, error_array, data_mask, headers,
             I_stokes, Q_stokes, U_stokes = Stokes_array
             Stokes_cov[0,0], Stokes_cov[1,1], Stokes_cov[2,2] = Stokes_error**2
 
-    return I_stokes, Q_stokes, U_stokes, Stokes_cov
+    return I_stokes, Q_stokes, U_stokes, Stokes_cov, dP_dtheta, dPA_dtheta
 
 
-def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers):
+def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, dP_dtheta, dPA_dtheta, headers):
     """
     Compute the polarization degree (in %) and angle (in deg) and their
     respective errors from given Stokes parameters.
@@ -1215,13 +1243,24 @@ def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers):
         print("WARNING : found {0:d} pixels for which P > 1".format(P[P>1.].size))
 
     #Associated errors
+    sigma_theta = np.array([3.*np.pi/180., 3.*np.pi/180., 3.*np.pi/180.])
     fmax = np.finfo(np.float64).max
-
     s_P = np.ones(I_stokes.shape)*fmax
-    s_P[mask] = (1/I_stokes[mask])*np.sqrt((Q_stokes[mask]**2*Stokes_cov[1,1][mask] + U_stokes[mask]**2*Stokes_cov[2,2][mask] + 2.*Q_stokes[mask]*U_stokes[mask]*Stokes_cov[1,2][mask])/(Q_stokes[mask]**2 + U_stokes[mask]**2) + ((Q_stokes[mask]/I_stokes[mask])**2 + (U_stokes[mask]/I_stokes[mask])**2)*Stokes_cov[0,0][mask] - 2.*(Q_stokes[mask]/I_stokes[mask])*Stokes_cov[0,1][mask] - 2.*(U_stokes[mask]/I_stokes[mask])*Stokes_cov[0,2][mask])
-    s_P[np.isnan(s_P)] = fmax
     s_PA = np.ones(I_stokes.shape)*fmax
+
+    # Propagate previously computed errors
+    s_P[mask] = (1/I_stokes[mask])*np.sqrt((Q_stokes[mask]**2*Stokes_cov[1,1][mask] + U_stokes[mask]**2*Stokes_cov[2,2][mask] + 2.*Q_stokes[mask]*U_stokes[mask]*Stokes_cov[1,2][mask])/(Q_stokes[mask]**2 + U_stokes[mask]**2) + ((Q_stokes[mask]/I_stokes[mask])**2 + (U_stokes[mask]/I_stokes[mask])**2)*Stokes_cov[0,0][mask] - 2.*(Q_stokes[mask]/I_stokes[mask])*Stokes_cov[0,1][mask] - 2.*(U_stokes[mask]/I_stokes[mask])*Stokes_cov[0,2][mask])
     s_PA[mask] = (90./(np.pi*(Q_stokes[mask]**2 + U_stokes[mask]**2)))*np.sqrt(U_stokes[mask]**2*Stokes_cov[1,1][mask] + Q_stokes[mask]**2*Stokes_cov[2,2][mask] - 2.*Q_stokes[mask]*U_stokes[mask]*Stokes_cov[1,2][mask])
+    # Compute the uncertainty associated with the polarizers' orientation (see Kishimoto 1999)
+    s_P_axis = np.sqrt(np.sum([dP_dtheta[i]**2 * sigma_theta[i]**2 for i in range(len(sigma_theta))], axis=0))
+    s_PA_axis = np.sqrt(np.sum([dPA_dtheta[i]**2 * sigma_theta[i]**2 for i in range(len(sigma_theta))], axis=0))
+    # Add axis uncertainties quadratically
+#    with warnings.catch_warnings(record=True) as w:
+#        s_P[mask] = np.sqrt(s_P[mask]**2 + s_P_axis[mask]**2)
+#        s_PA[mask] = np.sqrt(s_PA[mask]**2 + s_PA_axis[mask]**2)
+
+    s_P[np.isnan(s_P)] = fmax
+    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:

src/lib/test_overplot.py

@@ -12,9 +12,9 @@ levelsMorganti = np.array([1.,2.,3.,8.,16.,32.,64.,128.])
 #levels18GHz = np.array([0.6, 1.5, 3, 6, 12, 24, 48, 96])/100.*Stokes_18GHz[0].data.max()
 levels18GHz = levelsMorganti*0.28*1e-3
 A = overplot_maps(Stokes_UV, Stokes_18GHz)
-A.plot(levels=levels18GHz, SNRp_cut=6.0, SNRi_cut=180.0, savename='../../plots/IC5063_x3nl030/18GHz_overplot.png')
+A.plot(levels=levels18GHz, SNRp_cut=3.0, SNRi_cut=50.0, savename='../../plots/IC5063_x3nl030/18GHz_overplot.png')
 
 #levels24GHz = np.array([1.,1.5, 3, 6, 12, 24, 48, 96])/100.*Stokes_24GHz[0].data.max()
 levels24GHz = levelsMorganti*0.46*1e-3
 B = overplot_maps(Stokes_UV, Stokes_24GHz)
-B.plot(levels=levels24GHz, SNRp_cut=6.0, SNRi_cut=180.0, savename='../../plots/IC5063_x3nl030/24GHz_overplot.png')
+B.plot(levels=levels24GHz, SNRp_cut=3.0, SNRi_cut=50.0, savename='../../plots/IC5063_x3nl030/24GHz_overplot.png')