revamp imagerie

src/FOC_reduction.py

@@ -20,17 +20,17 @@ from astropy.wcs import WCS
 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']
-    psf_file = 'NGC1068_f253m00.fits'
+##    psf_file = 'NGC1068_f253m00.fits'
-    globals()['plots_folder'] = "../plots/NGC1068_x274020/"
+#    globals()['plots_folder'] = "../plots/NGC1068_x274020/"
 
-#    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']
 #    psf_file = 'IC5063_f502m00.fits'
-#    globals()['plots_folder'] = "../plots/IC5063_x3nl030/"
+    globals()['plots_folder'] = "../plots/IC5063_x3nl030/"
 
 #    globals()['data_folder'] = "../data/NGC1068_x14w010/"
 #    infiles = ['x14w0101t_c0f.fits','x14w0102t_c0f.fits','x14w0103t_c0f.fits',
@@ -122,10 +122,10 @@ def main():
     crop = False                    #Crop to desired ROI
     final_display = False
     # Polarization map output
-    figname = 'NGC1068_FOC'         #target/intrument name
+    figname = 'IC5063_FOC'         #target/intrument name
     figtype = '_combine_FWHM020'    #additionnal informations
-    SNRp_cut = 5.    #P measurments with SNR>3
+    SNRp_cut = 3.    #P measurments with SNR>3
-    SNRi_cut = 50.   #I measurments with SNR>30, which implies an uncertainty in P of 4.7%.
+    SNRi_cut = 60.   #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
                     # if step_vec = 0 then all vectors are displayed at full length
 
@@ -176,13 +176,13 @@ def main():
     # 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)
-    
+
     ## Step 3:
     # Rotate images to have North up
     if rotate_stokes:
         alpha = headers[0]['orientat']
         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, -alpha, 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)
 

src/lib/plots.py

@@ -271,7 +271,7 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
     if display is None:
         # If no display selected, show intensity map
         vmin, vmax = 0., np.max(stkI.data[stkI.data > 0.]*convert_flux)
-        im = ax.imshow(stkI.data*convert_flux, vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.)
+        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^{-1}$]")
         levelsI = np.linspace(vmax*0.01, vmax*0.99, 10)
         print("Total flux contour levels : ", levelsI)
@@ -281,7 +281,7 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
         # Display polarisation flux
         pf_mask = (stkI.data > 0.) * (pol.data > 0.)
         vmin, vmax = 0., np.max(stkI.data[pf_mask]*convert_flux*pol.data[pf_mask])
-        im = ax.imshow(stkI.data*convert_flux*pol.data, vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.)
+        im = ax.imshow(stkI.data*convert_flux*pol.data, vmin=vmin, vmax=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)
         print("Polarized flux contour levels : ", levelsPf)
@@ -290,24 +290,24 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
     elif display.lower() in ['p','pol','pol_deg']:
         # Display polarization degree map
         vmin, vmax = 0., 100.
-        im = ax.imshow(pol.data*100., vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.)
+        im = ax.imshow(pol.data*100., vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.)
         cbar = plt.colorbar(im, cax=cbar_ax, label=r"$P$ [%]")
     elif display.lower() in ['s_p','pol_err','pol_deg_err']:
         # Display polarization degree error map
         vmin, vmax = 0., 10.
         p_err = deepcopy(pol_err.data)
         p_err[p_err > vmax/100.] = np.nan
-        im = ax.imshow(p_err*100., vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.)
+        im = ax.imshow(p_err*100., vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.)
         cbar = plt.colorbar(im, cax=cbar_ax, label=r"$\sigma_P$ [%]")
     elif display.lower() in ['s_i','i_err']:
         # Display intensity error map
         vmin, vmax = 0., 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='auto', cmap='inferno', alpha=1.)
+        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
         vmin, vmax = 0., np.max(SNRi[SNRi > 0.])
-        im = ax.imshow(SNRi, vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.)
+        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)
         print("SNRi contour levels : ", levelsSNRi)
@@ -316,7 +316,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
         vmin, vmax = SNRp_cut, np.max(SNRp[SNRp > 0.])
-        im = ax.imshow(SNRp, vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.)
+        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)
         print("SNRp contour levels : ", levelsSNRp)
@@ -325,7 +325,7 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
     else:
         # Defaults to intensity map
         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.)
+        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$]")
 
     if (display is None) or not(display.lower() in ['default']):
@@ -398,7 +398,7 @@ class align_maps(object):
         elif self.wcs_map.naxis == 3:
             self.wcs_map = WCS(self.map[0],naxis=[1,2]).deepcopy()
             self.map[0].data = self.map[0].data[1]
-        
+
         self.wcs_other = WCS(self.other_map[0]).deepcopy()
         if self.wcs_other.naxis == 4:
             self.wcs_other = WCS(self.other_map[0],naxis=[1,2]).deepcopy()
@@ -406,7 +406,7 @@ class align_maps(object):
         elif self.wcs_other.naxis == 3:
             self.wcs_other = WCS(self.other_map[0],naxis=[1,2]).deepcopy()
             self.other_map[0].data = self.other_map[0].data[1]
-        
+
         try:
             convert_flux = self.map[0].header['photflam']
         except KeyError:
@@ -415,7 +415,7 @@ class align_maps(object):
             other_convert = self.other_map[0].header['photflam']
         except KeyError:
             other_convert = 1.
-        
+
         #Get data
         data = self.map[0].data
         other_data = self.other_map[0].data
@@ -427,13 +427,13 @@ class align_maps(object):
         self.ax1.set_facecolor('k')
 
         vmin, vmax = 0., np.max(data[data > 0.]*convert_flux)
-        im1 = self.ax1.imshow(data*convert_flux, vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.)
+        im1 = self.ax1.imshow(data*convert_flux, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.)
 
         fontprops = fm.FontProperties(size=16)
         px_size = self.wcs_map.wcs.get_cdelt()[0]*3600.
         px_sc = AnchoredSizeBar(self.ax1.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.ax1.add_artist(px_sc)
-        
+
         try:
             north_dir1 = AnchoredDirectionArrows(self.ax1.transAxes, "E", "N", length=-0.08, fontsize=0.025, loc=1, aspect_ratio=-1, sep_y=0.01, sep_x=0.01, back_length=0., head_length=10., head_width=10., angle=-self.map[0].header['orientat'], color='white', text_props={'ec': None, 'fc': 'w', 'alpha': 1, 'lw': 0.4}, arrow_props={'ec': None,'fc':'w','alpha': 1,'lw': 1})
             self.ax1.add_artist(north_dir1)
@@ -454,14 +454,14 @@ class align_maps(object):
                 test = kwargs[key]
             except KeyError:
                 for key_i, val_i in value:
-                    kwargs[key_i] = val_i                    
+                    kwargs[key_i] = val_i
-        im2 = self.ax2.imshow(other_data*other_convert, aspect='auto', **kwargs)
+        im2 = self.ax2.imshow(other_data*other_convert, aspect='equal', **kwargs)
 
         fontprops = fm.FontProperties(size=16)
         px_size = self.wcs_other.wcs.get_cdelt()[0]*3600.
         px_sc = AnchoredSizeBar(self.ax2.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.ax2.add_artist(px_sc)
-        
+
         try:
             north_dir2 = AnchoredDirectionArrows(self.ax2.transAxes, "E", "N", length=-0.08, fontsize=0.03, loc=1, aspect_ratio=-1, sep_y=0.01, sep_x=0.01, angle=-self.other_map[0].header['orientat'], color='w', arrow_props={'ec': None, 'fc': 'w', 'alpha': 1,'lw': 2})
             self.ax2.add_artist(north_dir2)
@@ -477,7 +477,7 @@ class align_maps(object):
         self.bapply = Button(self.axapply, 'Apply reference')
         self.axreset = self.fig.add_axes([0.60, 0.01, 0.1, 0.04])
         self.breset = Button(self.axreset, 'Leave as is')
-     
+
     def get_aligned_wcs(self):
         return self.wcs_map, self.wcs_other
 
@@ -489,14 +489,14 @@ class align_maps(object):
 
                 self.cr_map.set(data=[x,y])
                 self.fig.canvas.draw_idle()
-            
+
             if (event.inaxes is not None) and (event.inaxes == self.ax2):
                 x = event.xdata
                 y = event.ydata
 
                 self.cr_other.set(data=[x,y])
                 self.fig.canvas.draw_idle()
-    
+
     def reset_align(self, event):
         self.wcs_map.wcs.crpix = WCS(self.map[0].header).wcs.crpix[:2]
         self.wcs_other.wcs.crpix = WCS(self.other_map[0].header).wcs.crpix[:2]
@@ -504,7 +504,7 @@ class align_maps(object):
 
         if self.aligned:
             plt.close()
-        
+
         self.aligned = True
 
     def apply_align(self, event):
@@ -515,13 +515,13 @@ class align_maps(object):
 
         if self.aligned:
             plt.close()
-        
+
         self.aligned = True
-    
+
     def on_close_align(self, event):
         self.aligned = True
         #print(self.get_aligned_wcs())
-    
+
     def align(self):
         self.fig.canvas.draw()
         self.fig.canvas.mpl_connect('button_press_event', self.onclick_ref)
@@ -547,7 +547,7 @@ class overplot_radio(align_maps):
         pol = self.Stokes_UV[np.argmax([self.Stokes_UV[i].header['datatype']=='Pol_deg_debiased' for i in range(len(self.Stokes_UV))])]
         pol_err = self.Stokes_UV[np.argmax([self.Stokes_UV[i].header['datatype']=='Pol_deg_err' for i in range(len(self.Stokes_UV))])]
         pang = self.Stokes_UV[np.argmax([self.Stokes_UV[i].header['datatype']=='Pol_ang' for i in range(len(self.Stokes_UV))])]
-        
+
         other_data = self.other_map[0].data
         other_convert = 1.
         other_unit = self.other_map[0].header['bunit']
@@ -555,7 +555,7 @@ class overplot_radio(align_maps):
             other_unit = r"mJy/Beam"
             other_convert = 1e3
         other_freq = self.other_map[0].header['crval3']
-        
+
         convert_flux = self.Stokes_UV[0].header['photflam']
 
         #Compute SNR and apply cuts
@@ -575,7 +575,7 @@ class overplot_radio(align_maps):
 
         #Display UV intensity map with polarization vectors
         vmin, vmax = 0., np.max(stkI.data[stkI.data > 0.]*convert_flux)
-        im = self.ax.imshow(stkI.data*convert_flux, vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.)
+        im = self.ax.imshow(stkI.data*convert_flux, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.)
         cbar_ax = self.fig2.add_axes([0.95, 0.12, 0.01, 0.75])
         cbar = plt.colorbar(im, cax=cbar_ax, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]")
 
@@ -600,12 +600,12 @@ class overplot_radio(align_maps):
         north_dir = AnchoredDirectionArrows(self.ax.transAxes, "E", "N", length=-0.08, fontsize=0.03, loc=1, aspect_ratio=-1, sep_y=0.01, sep_x=0.01, angle=-self.Stokes_UV[0].header['orientat'], color='w', arrow_props={'ec': None, 'fc': 'w', 'alpha': 1,'lw': 2})
         self.ax.add_artist(north_dir)
 
- 
+
         if not(savename is None):
             self.fig2.savefig(savename,bbox_inches='tight',dpi=200)
 
         self.fig2.canvas.draw()
-    
+
     def plot(self, levels, SNRp_cut=3., SNRi_cut=30., savename=None) -> None:
         self.align()
         if self.aligned:
@@ -628,9 +628,9 @@ class overplot_pol(align_maps):
         pol = self.Stokes_UV[np.argmax([self.Stokes_UV[i].header['datatype']=='Pol_deg_debiased' for i in range(len(self.Stokes_UV))])]
         pol_err = self.Stokes_UV[np.argmax([self.Stokes_UV[i].header['datatype']=='Pol_deg_err' for i in range(len(self.Stokes_UV))])]
         pang = self.Stokes_UV[np.argmax([self.Stokes_UV[i].header['datatype']=='Pol_ang' for i in range(len(self.Stokes_UV))])]
-        
+
         convert_flux = self.Stokes_UV[0].header['photflam']
-        
+
         other_data = self.other_map[0].data
         try:
             other_convert = self.other_map[0].header['photflam']
@@ -654,9 +654,9 @@ class overplot_pol(align_maps):
 
         #Display Stokes I as contours
         levels_stkI = np.rint(np.linspace(10,99,10))/100.*np.max(stkI.data[stkI.data > 0.]*convert_flux)
-        cont_stkI = self.ax.contour(stkI.data*convert_flux, transform=self.ax.get_transform(self.wcs_UV), levels=levels_stkI, colors='grey')
+        cont_stkI = self.ax.contour(stkI.data*convert_flux, transform=self.ax.get_transform(self.wcs_UV), levels=levels_stkI, colors='grey', alpha=0.)
-        self.ax.clabel(cont_stkI, inline=True, fontsize=8)
+        #self.ax.clabel(cont_stkI, inline=True, fontsize=8)
-        
+
         self.ax.autoscale(False)
 
         #Display full size polarization vectors
@@ -673,27 +673,27 @@ class overplot_pol(align_maps):
                 test = kwargs[key]
             except KeyError:
                 for key_i, val_i in value:
-                    kwargs[key_i] = val_i                    
+                    kwargs[key_i] = val_i
         im = self.ax.imshow(other_data*other_convert, transform=self.ax.get_transform(self.wcs_other), alpha=1., **kwargs)
         cbar_ax = self.fig2.add_axes([0.95, 0.12, 0.01, 0.75])
         cbar = plt.colorbar(im, cax=cbar_ax, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]")
 
         #Display pixel scale and North direction
         fontprops = fm.FontProperties(size=16)
-        px_size = self.wcs_other.wcs.get_cdelt()[0]*3600.
+        px_size = self.wcs_UV.wcs.get_cdelt()[0]*3600.
         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)
         north_dir = AnchoredDirectionArrows(self.ax.transAxes, "E", "N", length=-0.08, fontsize=0.03, loc=1, aspect_ratio=-1, sep_y=0.01, sep_x=0.01, angle=-self.Stokes_UV[0].header['orientat'], color='w', arrow_props={'ec': None, 'fc': 'w', 'alpha': 1,'lw': 2})
         self.ax.add_artist(north_dir)
 
-        
+
         self.ax.set(xlabel="Right Ascension (J2000)", ylabel="Declination (J2000)", title="{0:s} overplotted with polarization vectors and Stokes I contours from HST/FOC".format(obj))
 
         if not(savename is None):
             self.fig2.savefig(savename,bbox_inches='tight',dpi=200)
 
         self.fig2.canvas.draw()
-    
+
     def plot(self, SNRp_cut=3., SNRi_cut=30., savename=None, **kwargs) -> None:
         self.align()
         if self.aligned:
@@ -711,7 +711,7 @@ class crop_map(object):
         self.hdul = hdul
         self.header = deepcopy(self.hdul[0].header)
         self.wcs = WCS(self.header).deepcopy()
-        
+
         self.data = deepcopy(self.hdul[0].data)
         try:
             self.convert_flux = self.header['photflam']
@@ -758,16 +758,16 @@ class crop_map(object):
             wcs = self.wcs
         if convert_flux is None:
             convert_flux = self.convert_flux
-        
+
         vmin, vmax = 0., np.max(data[data > 0.]*convert_flux)
         if hasattr(self, 'im'):
             self.im.remove()
-        self.im = self.ax.imshow(data*convert_flux, vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=self.mask_alpha, origin='lower')
+        self.im = self.ax.imshow(data*convert_flux, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=self.mask_alpha, origin='lower')
         if hasattr(self, 'cr'):
             self.cr[0].set_data(*wcs.wcs.crpix)
         else:
             self.cr = self.ax.plot(*wcs.wcs.crpix, 'r+')
-        self.fig.canvas.draw_idle()        
+        self.fig.canvas.draw_idle()
         return self.im
 
     @property
@@ -811,10 +811,10 @@ class crop_map(object):
             header = self.header
             data = self.data
             wcs = self.wcs
-        
+
         vertex = self.RSextent.astype(int)
         shape = vertex[1::2] - vertex[0::2]
-        
+
         extent = np.array(self.im.get_extent())
         shape_im = extent[1::2] - extent[0::2]
         if (shape_im.astype(int) != shape).any() and (self.RSextent != self.extent).any():
@@ -827,7 +827,7 @@ class crop_map(object):
             else:
                 self.wcs_crop.wcs.crval = wcs.wcs_pix2world([self.RScenter],1)[0]
                 self.wcs_crop.wcs.crpix = self.RScenter-self.RSextent[::2]
-             
+
             # Crop dataset
             self.data_crop = deepcopy(data[vertex[2]:vertex[3], vertex[0]:vertex[1]])
 
@@ -853,7 +853,7 @@ class crop_map(object):
                 self.rect_selector = RectangleSelector(self.ax, self.onselect_crop,
                         drawtype='box', button=[1], interactive=True)
         self.fig.canvas.draw_idle()
-    
+
     def on_close(self, event) -> None:
         if not hasattr(self, 'hdul_crop'):
             self.hdul_crop = self.hdul
@@ -890,7 +890,7 @@ class crop_Stokes(crop_map):
             hdul = self.hdul
             data = self.data
             wcs = self.wcs
-        
+
         vertex = self.RSextent.astype(int)
         shape = vertex[1::2] - vertex[0::2]
 
@@ -907,7 +907,7 @@ class crop_Stokes(crop_map):
             else:
                 self.wcs_crop.wcs.crval = wcs.wcs_pix2world([self.RScenter],1)[0]
                 self.wcs_crop.wcs.crpix = self.RScenter-self.RSextent[::2]
-             
+
             # Crop dataset
             for dataset in self.hdul_crop:
                 if dataset.header['datatype']=='IQU_cov_matrix':
@@ -941,7 +941,7 @@ class crop_Stokes(crop_map):
                 self.rect_selector = RectangleSelector(self.ax, self.onselect_crop,
                         drawtype='box', button=[1], interactive=True)
         self.fig.canvas.draw_idle()
-    
+
     @property
     def data_mask(self):
         return self.hdul_crop[-1].data
@@ -967,13 +967,13 @@ class image_lasso_selector:
             self.ax = ax
             self.mask_alpha = 0.1
             self.embedded = True
-        self.displayed = self.ax.imshow(self.img, vmin=self.vmin, vmax=self.vmax, aspect='auto', cmap='inferno',alpha=self.mask_alpha)
+        self.displayed = self.ax.imshow(self.img, vmin=self.vmin, vmax=self.vmax, aspect='equal', cmap='inferno',alpha=self.mask_alpha)
         plt.ion()
-        
+
         lineprops = {'color': 'grey', 'linewidth': 1, 'alpha': 0.8}
         self.lasso = LassoSelector(self.ax, self.onselect,props=lineprops, useblit=False)
         self.lasso.set_visible(True)
-        
+
         pix_x = np.arange(self.img.shape[0])
         pix_y = np.arange(self.img.shape[1])
         xv, yv = np.meshgrid(pix_y,pix_x)
@@ -993,10 +993,10 @@ class image_lasso_selector:
         p = Path(verts)
         self.indices = p.contains_points(self.pix, radius=0).reshape(self.img.shape[:2])
         self.update_mask()
-        
+
     def update_mask(self):
         self.displayed.remove()
-        self.displayed = self.ax.imshow(self.img, vmin=self.vmin, vmax=self.vmax, aspect='auto', cmap='inferno',alpha=self.mask_alpha)
+        self.displayed = self.ax.imshow(self.img, vmin=self.vmin, vmax=self.vmax, aspect='equal', cmap='inferno',alpha=self.mask_alpha)
         array = self.displayed.get_array().data
 
         self.mask = np.zeros(self.img.shape[:2],dtype=bool)
@@ -1040,14 +1040,14 @@ class pol_map(object):
         self.ax = self.fig.add_subplot(111,projection=self.wcs)
         self.ax_cosmetics()
         self.cbar_ax = self.fig.add_axes([0.925, 0.12, 0.01, 0.75])
- 
+
         #Display selected data (Default to total flux)
         self.display()
         #Display polarization vectors in SNR_cut
         self.pol_vector()
         #Display integrated values in ROI
         self.pol_int()
-        
+
         #Set axes for sliders (SNRp_cut, SNRi_cut)
         ax_I_cut = self.fig.add_axes([0.125, 0.080, 0.35, 0.01])
         ax_P_cut = self.fig.add_axes([0.125, 0.055, 0.35, 0.01])
@@ -1121,7 +1121,7 @@ class pol_map(object):
         b_crop = Button(ax_crop,"Crop")
         self.cropped = False
         b_crop_reset = Button(ax_crop_reset,"Reset")
-        
+
         def crop(event):
             if self.cropped:
                 self.cropped = False
@@ -1284,7 +1284,7 @@ class pol_map(object):
         b_snrp.on_clicked(d_snrp)
 
         plt.show()
-    
+
     @property
     def wcs(self):
         return deepcopy(WCS(self.Stokes[0].header))
@@ -1312,7 +1312,7 @@ class pol_map(object):
     @property
     def data_mask(self):
         return self.Stokes[np.argmax([self.Stokes[i].header['datatype']=='Data_mask' for i in range(len(self.Stokes))])].data
-    
+
     def set_data_mask(self, mask):
         self.Stokes[np.argmax([self.Stokes[i].header['datatype']=='Data_mask' for i in range(len(self.Stokes))])].data = mask.astype(float)
 
@@ -1337,7 +1337,7 @@ class pol_map(object):
         ax.coords[1].set_axislabel_position('l')
         ax.coords[1].set_ticklabel_position('l')
         ax.axis('equal')
-        
+
         #Display scales and orientation
         fontprops = fm.FontProperties(size=14)
         px_size = self.wcs.wcs.cdelt[0]*3600.
@@ -1357,7 +1357,7 @@ class pol_map(object):
     def display(self, fig=None, ax=None):
         if self.display_selection is None:
             self.display_selection = "total_flux"
-        
+
         if self.display_selection.lower() in ['total_flux']:
             self.data = self.I*self.convert_flux
             vmin, vmax = 0., np.max(self.data[self.data > 0.])
@@ -1390,12 +1390,12 @@ class pol_map(object):
                 ax = self.ax
             if hasattr(self, 'im'):
                 self.im.remove()
-            self.im = ax.imshow(self.data, vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno')
+            self.im = ax.imshow(self.data, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno')
             self.cbar = plt.colorbar(self.im, cax=self.cbar_ax, label=label)
             fig.canvas.draw_idle()
             return self.im
         else:
-            im = ax.imshow(self.data, vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno')
+            im = ax.imshow(self.data, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno')
             ax.set_xlim(0,self.data.shape[1])
             ax.set_ylim(0,self.data.shape[0])
             plt.colorbar(im, pad=0.025, aspect=80, label=label)
@@ -1419,7 +1419,7 @@ class pol_map(object):
         else:
             ax.quiver(X, Y, XY_U, XY_V, units='xy', scale=0.5, scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.1, color='white')
             fig.canvas.draw_idle()
-    
+
     def pol_int(self, fig=None, ax=None):
         if self.region is None:
             n_pix = self.I.size
@@ -1475,4 +1475,3 @@ class pol_map(object):
             if not self.region is None:
                 ax.contour(self.region.astype(float),levels=[0.5], colors='white', linewidths=0.8)
             fig.canvas.draw_idle()
-   

src/lib/reduction.py

@@ -482,7 +482,7 @@ def get_error(data_array, headers, error_array=None, data_mask=None, sub_shape=N
         err_flat = data_array[i]*0.03
 
         error_array[i] = np.sqrt(error_array[i]**2 + error_bkg[i]**2 + err_wav**2 + err_psf**2 + err_flat**2)
-        
+
         background[i] = sub_image.sum()
         if (data_array[i] < 0.).any():
             print(data_array[i])
@@ -808,7 +808,7 @@ def align_data(data_array, headers, error_array=None, upsample_factor=1.,
 
         shifts.append(shift)
         errors.append(error)
-    
+
     shifts = np.array(shifts)
     errors = np.array(errors)
 
@@ -818,7 +818,7 @@ def align_data(data_array, headers, error_array=None, upsample_factor=1.,
     for i in range(len(headers_wcs)):
         headers_wcs[i].wcs.crpix = new_crpix[0]
         headers[i].update(headers_wcs[i].to_header())
-    
+
     data_mask = rescaled_mask.all(axis=0)
     data_array, error_array, data_mask, headers = crop_array(rescaled_image, headers, rescaled_error, data_mask)
 
@@ -926,7 +926,7 @@ def smooth_data(data_array, error_array, data_mask, headers, FWHM=1.,
 
     else:
         raise ValueError("{} is not a valid smoothing option".format(smoothing))
-    
+
     return smoothed, error
 
 
@@ -1039,7 +1039,7 @@ def polarizer_avg(data_array, error_array, data_mask, headers, FWHM=None,
             err60 = np.sqrt(np.sum(err60_array**2,axis=0))
             err120 = np.sqrt(np.sum(err120_array**2,axis=0))
             polerr_array = np.array([err0, err60, err120])
-            
+
             # Update headers
             for header in headers:
                 if header['filtnam1']=='POL0':
@@ -1186,12 +1186,12 @@ def compute_Stokes(data_array, error_array, data_mask, headers,
         Q_stokes = np.zeros(pol_array[0].shape)
         U_stokes = np.zeros(pol_array[0].shape)
         Stokes_cov = np.zeros((3,3,I_stokes.shape[0],I_stokes.shape[1]))
-        
+
         for i in range(I_stokes.shape[0]):
             for j in range(I_stokes.shape[1]):
                 I_stokes[i,j], Q_stokes[i,j], U_stokes[i,j] = np.dot(coeff_stokes, pol_flux[:,i,j]).T
                 Stokes_cov[:,:,i,j] = np.dot(coeff_stokes, np.dot(pol_cov[:,:,i,j], coeff_stokes.T))
-        
+
         mask = (Q_stokes**2 + U_stokes**2) > I_stokes**2
         if mask.any():
             print("WARNING : found {0:d} pixels for which I_pol > I_stokes".format(I_stokes[mask].size))
@@ -1201,17 +1201,17 @@ def compute_Stokes(data_array, error_array, data_mask, headers,
         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 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)
@@ -1221,7 +1221,7 @@ def compute_Stokes(data_array, error_array, data_mask, headers,
         Stokes_cov[0,0] += s_I2_axis
         Stokes_cov[1,1] += s_Q2_axis
         Stokes_cov[2,2] += s_U2_axis
-        
+
         if not(FWHM is None) and (smoothing.lower() in ['weighted_gaussian_after','weight_gauss_after','gaussian_after','gauss_after']):
             smoothing = smoothing.lower()[:-6]
             Stokes_array = np.array([I_stokes, Q_stokes, U_stokes])
@@ -1233,7 +1233,7 @@ 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
-        
+
         #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
@@ -1419,11 +1419,11 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers,
     mrot = np.array([[1., 0., 0.],
                     [0., np.cos(2.*alpha), np.sin(2.*alpha)],
                     [0, -np.sin(2.*alpha), np.cos(2.*alpha)]])
-    
+
     old_center = np.array(I_stokes.shape)/2
     shape = np.fix(np.array(I_stokes.shape)*np.sqrt(2.5)).astype(int)
     new_center = np.array(shape)/2
-    
+
     I_stokes = zeropad(I_stokes, shape)
     Q_stokes = zeropad(Q_stokes, shape)
     U_stokes = zeropad(U_stokes, shape)
@@ -1433,7 +1433,7 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers,
     new_Q_stokes = np.zeros(shape)
     new_U_stokes = np.zeros(shape)
     new_Stokes_cov = np.zeros((*Stokes_cov.shape[:-2],*shape))
-    
+
     for i in range(shape[0]):
         for j in range(shape[1]):
             new_I_stokes[i,j], new_Q_stokes[i,j], new_U_stokes[i,j] = np.dot(mrot, np.array([I_stokes[i,j], Q_stokes[i,j], U_stokes[i,j]])).T
@@ -1554,7 +1554,7 @@ def rotate_data(data_array, error_array, data_mask, headers, ang):
     old_center = np.array(data_array[0].shape)/2
     shape = np.fix(np.array(data_array[0].shape)*np.sqrt(2.5)).astype(int)
     new_center = np.array(shape)/2
-    
+
     data_array = zeropad(data_array, [data_array.shape[0],*shape])
     error_array = zeropad(error_array, [error_array.shape[0],*shape])
     data_mask = zeropad(data_mask, shape)

src/overplot.py

@@ -16,30 +16,30 @@ Stokes_IR = fits.open("../data/IC5063_x3nl030/IR/u2e65g01t_c0f_rot.fits")
 
 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 = np.array([0.6, 1.5, 3, 6, 12, 24, 48, 96])/100.*Stokes_18GHz[0].data.max()
-levels18GHz = levelsMorganti*0.28*1e-3
+#levels18GHz = levelsMorganti*0.28*1e-3
-A = overplot_radio(Stokes_UV, Stokes_18GHz)
+#A = overplot_radio(Stokes_UV, Stokes_18GHz)
-A.plot(levels=levels18GHz, SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/18GHz_overplot_forced.png')
+#A.plot(levels=levels18GHz, SNRp_cut=3.0, SNRi_cut=60.0, savename='../plots/IC5063_x3nl030/18GHz_overplot_forced.png')
 
-#levels24GHz = np.array([1.,1.5, 3, 6, 12, 24, 48, 96])/100.*Stokes_24GHz[0].data.max()
+##levels24GHz = np.array([1.,1.5, 3, 6, 12, 24, 48, 96])/100.*Stokes_24GHz[0].data.max()
-levels24GHz = levelsMorganti*0.46*1e-3
+#levels24GHz = levelsMorganti*0.46*1e-3
-B = overplot_radio(Stokes_UV, Stokes_24GHz)
+#B = overplot_radio(Stokes_UV, Stokes_24GHz)
-B.plot(levels=levels24GHz, SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/24GHz_overplot_forced.png')
+#B.plot(levels=levels24GHz, SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/24GHz_overplot_forced.png')
-
+#
-levels103GHz = np.linspace(1,99,11)/100.*np.max(deepcopy(Stokes_103GHz[0].data[Stokes_103GHz[0].data > 0.]))
+#levels103GHz = np.linspace(1,99,11)/100.*np.max(deepcopy(Stokes_103GHz[0].data[Stokes_103GHz[0].data > 0.]))
-C = overplot_radio(Stokes_UV, Stokes_103GHz)
+#C = overplot_radio(Stokes_UV, Stokes_103GHz)
-C.plot(levels=levels103GHz, SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/103GHz_overplot_forced.png')
+#C.plot(levels=levels103GHz, SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/103GHz_overplot_forced.png')
-
+#
-levels229GHz = np.linspace(1,99,11)/100.*np.max(deepcopy(Stokes_229GHz[0].data[Stokes_229GHz[0].data > 0.]))
+#levels229GHz = np.linspace(1,99,11)/100.*np.max(deepcopy(Stokes_229GHz[0].data[Stokes_229GHz[0].data > 0.]))
-D = overplot_radio(Stokes_UV, Stokes_229GHz)
+#D = overplot_radio(Stokes_UV, Stokes_229GHz)
-D.plot(levels=levels229GHz, SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/229GHz_overplot_forced.png')
+#D.plot(levels=levels229GHz, SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/229GHz_overplot_forced.png')
-
+#
-levels357GHz = np.linspace(1,99,11)/100.*np.max(deepcopy(Stokes_357GHz[0].data[Stokes_357GHz[0].data > 0.]))
+#levels357GHz = np.linspace(1,99,11)/100.*np.max(deepcopy(Stokes_357GHz[0].data[Stokes_357GHz[0].data > 0.]))
-E = overplot_radio(Stokes_UV, Stokes_357GHz)
+#E = overplot_radio(Stokes_UV, Stokes_357GHz)
-E.plot(levels=levels357GHz, SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/357GHz_overplot_forced.png')
+#E.plot(levels=levels357GHz, SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/357GHz_overplot_forced.png')
-
+#
-F = overplot_pol(Stokes_UV, Stokes_S2)
+#F = overplot_pol(Stokes_UV, Stokes_S2)
-F.plot(SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/S2_overplot_forced.png', norm=LogNorm(vmin=5e-20,vmax=5e-18))
+#F.plot(SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/S2_overplot_forced.png', norm=LogNorm(vmin=5e-20,vmax=5e-18))
 
 G = overplot_pol(Stokes_UV, Stokes_IR, norm=LogNorm(vmin=1e-17,vmax=5e-15), cmap='inferno_r')
-G.plot(SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/IR_overplot_forced.png', norm=LogNorm(vmin=1e-17,vmax=5e-15), cmap='inferno_r')
+G.plot(SNRp_cut=3.0, SNRi_cut=60.0, savename='../plots/IC5063_x3nl030/IR_overplot_forced.png', norm=LogNorm(vmin=1e-17,vmax=5e-15), cmap='inferno_r')