Save only data in "data_mask"

src/FOC_reduction.py

@@ -116,7 +116,7 @@ def main():
     rotate_stokes = True           #rotation to North convention can give erroneous results
     rotate_data = False              #rotation to North convention can give erroneous results
     # Final crop
-    crop = True        #Crop to desired ROI
+    crop = False        #Crop to desired ROI
     # Polarization map output
     figname = 'NGC1068_FOC'         #target/intrument name
     figtype = '_combine_FWHM020'    #additionnal informations
@@ -188,6 +188,7 @@ def main():
     ## Step 4:
     # Save image to FITS.
     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, data_mask, figname+figtype, data_folder=data_folder, return_hdul=True)
+    data_mask = Stokes_test[-1].data.astype(bool)
 
     ## Step 5:
     # crop to desired region of interest (roi)

src/lib/fits.py

@@ -12,6 +12,7 @@ prototypes :
 import numpy as np
 from astropy.io import fits
 from astropy import wcs
+from lib.convex_hull import image_hull
 
 
 def get_obs_data(infiles, data_folder="", compute_flux=False):
@@ -77,6 +78,33 @@ def get_obs_data(infiles, data_folder="", compute_flux=False):
     return data_array, headers
 
 
+def clean_ROI(image):
+    H,J = [],[]
+
+    shape = np.array(image.shape)
+    row, col = np.indices(shape)
+
+    for i in range(0,shape[0]):
+        r = row[i,:][image[i,:]>0.]
+        c = col[i,:][image[i,:]>0.]
+        if len(r)>1 and len(c)>1:
+            H.append((r[0],c[0]))
+            H.append((r[-1],c[-1]))
+    H = np.array(H)
+    for j in range(0,shape[1]):
+        r = row[:,j][image[:,j]>0.]
+        c = col[:,j][image[:,j]>0.]
+        if len(r)>1 and len(c)>1:
+            J.append((r[0],c[0]))
+            J.append((r[-1],c[-1]))
+    J = np.array(J)
+    xmin = np.min([H[:,1].min(),J[:,1].min()])-1
+    xmax = np.max([H[:,1].max(),J[:,1].max()])+1
+    ymin = np.min([H[:,0].min(),J[:,0].min()])-1
+    ymax = np.max([H[:,0].max(),J[:,0].max()])+1
+    return np.array([xmin,xmax,ymin,ymax])
+
+
 def 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, data_mask, filename, data_folder="",
         return_hdul=False):
@@ -116,6 +144,11 @@ def save_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, P, debiased_P, s_P,
     ref_header = headers[0]
     exp_tot = np.array([header['exptime'] for header in headers]).sum()
     new_wcs = wcs.WCS(ref_header).deepcopy()
+    
+    vertex = clean_ROI(1.-data_mask)
+    shape = vertex[1::2]-vertex[0::2]
+    new_wcs.array_shape = shape
+    new_wcs.wcs.crpix -= vertex[0::2]
 
     header = new_wcs.to_header()
     header['instrume'] = (ref_header['instrume'], 'Instrument from which data is reduced')
@@ -131,6 +164,27 @@ def save_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, P, debiased_P, s_P,
     header['PA_int'] = (ref_header['PA_int'], 'Integrated polarization angle')
     header['PA_int_err'] = (ref_header['PA_int_err'], 'Integrated polarization angle error')
 
+    #Crop Data to mask
+    I_stokes = I_stokes[vertex[0]:vertex[1],vertex[2]:vertex[3]]
+    Q_stokes = Q_stokes[vertex[0]:vertex[1],vertex[2]:vertex[3]]
+    U_stokes = U_stokes[vertex[0]:vertex[1],vertex[2]:vertex[3]]
+    P = P[vertex[0]:vertex[1],vertex[2]:vertex[3]]
+    debiased_P = debiased_P[vertex[0]:vertex[1],vertex[2]:vertex[3]]
+    s_P = s_P[vertex[0]:vertex[1],vertex[2]:vertex[3]]
+    s_P_P = s_P_P[vertex[0]:vertex[1],vertex[2]:vertex[3]]
+    PA = PA[vertex[0]:vertex[1],vertex[2]:vertex[3]]
+    s_PA = s_PA[vertex[0]:vertex[1],vertex[2]:vertex[3]]
+    s_PA_P = s_PA_P[vertex[0]:vertex[1],vertex[2]:vertex[3]]
+    
+    new_Stokes_cov = np.zeros((3,3,shape[0],shape[1]))
+    for i in range(3):
+        for j in range(3):
+            Stokes_cov[i,j][data_mask] = 0.
+            new_Stokes_cov[i,j] = Stokes_cov[i,j][vertex[0]:vertex[1],vertex[2]:vertex[3]]
+    Stokes_cov = new_Stokes_cov
+    
+    data_mask = data_mask[vertex[0]:vertex[1],vertex[2]:vertex[3]].astype(float, copy=False)
+
     #Create HDUList object
     hdul = fits.HDUList([])
 
@@ -139,8 +193,6 @@ def save_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, P, debiased_P, s_P,
     primary_hdu = fits.PrimaryHDU(data=I_stokes, header=header)
     hdul.append(primary_hdu)
 
-    data_mask = data_mask.astype(float, copy=False)
-
     #Add Q, U, Stokes_cov, P, s_P, PA, s_PA to the HDUList
     for data, name in [[Q_stokes,'Q_stokes'],[U_stokes,'U_stokes'],
             [Stokes_cov,'IQU_cov_matrix'],[P, 'Pol_deg'],

src/lib/plots.py

@@ -242,6 +242,7 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
 
     #Compute SNR and apply cuts
     pol.data[pol.data == 0.] = np.nan
+    pol_err.data[pol_err.data == 0.] = np.nan
     SNRp = pol.data/pol_err.data
     SNRp[np.isnan(SNRp)] = 0.
     pol.data[SNRp < SNRp_cut] = np.nan
@@ -386,39 +387,55 @@ class align_maps(object):
     """
     Class to interactively align maps with different WCS.
     """
-    def __init__(self, Stokes, other_map):
+    def __init__(self, map1, other_map):
         self.aligned = False
-        self.Stokes_UV = Stokes
+        self.map = map1
         self.other_map = other_map
 
-        self.wcs_UV = WCS(self.Stokes_UV[0]).deepcopy()
-        convert_flux = self.Stokes_UV[0].header['photflam']
+        self.wcs_map = WCS(self.map[0]).deepcopy()
+        if self.wcs_map.naxis > 2:
+            self.wcs_map = WCS(self.map[0],naxis=[1,2]).deepcopy()
+            self.map[0].data = self.map[0].data[0,0]
+        
         self.wcs_other = WCS(self.other_map[0]).deepcopy()
         if self.wcs_other.naxis > 2:
             self.wcs_other = WCS(self.other_map[0],naxis=[1,2]).deepcopy()
             self.other_map[0].data = self.other_map[0].data[0,0]
+        
+        try:
+            convert_flux = self.map[0].header['photflam']
+        except KeyError:
+            convert_flux = 1.
+        try:
+            other_convert = self.other_map[0].header['photflam']
+        except KeyError:
+            other_convert = 1.
+        
         #Get data
-        stkI = self.Stokes_UV[np.argmax([self.Stokes_UV[i].header['datatype']=='I_stokes' for i in range(len(self.Stokes_UV))])]
+        data = self.map[0].data
         other_data = self.other_map[0].data
 
         plt.rcParams.update({'font.size': 16})
         self.fig = plt.figure(figsize=(25,15))
         #Plot the UV map
-        self.ax1 = self.fig.add_subplot(121, projection=self.wcs_UV)
+        self.ax1 = self.fig.add_subplot(121, projection=self.wcs_map)
         self.ax1.set_facecolor('k')
 
-        vmin, vmax = 0., np.max(stkI.data[stkI.data > 0.]*convert_flux)
-        im1 = self.ax1.imshow(stkI.data*convert_flux, vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.)
+        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.)
 
         fontprops = fm.FontProperties(size=16)
-        px_size = self.wcs_UV.wcs.get_cdelt()[0]*3600.
+        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)
         
-        north_dir1 = AnchoredDirectionArrows(self.ax1.transAxes, "E", "N", length=-0.08, fontsize=0.03, loc=1, aspect_ratio=-1, sep_y=0.01, sep_x=0.01, angle=-Stokes[0].header['orientat'], color='w', arrow_props={'ec': 'w', 'fc': 'w', 'alpha': 1,'lw': 2})
-        self.ax1.add_artist(north_dir1)
+        try:
+            north_dir1 = AnchoredDirectionArrows(self.ax1.transAxes, "E", "N", length=-0.08, fontsize=0.03, loc=1, aspect_ratio=-1, sep_y=0.01, sep_x=0.01, angle=-self.map[0].header['orientat'], color='w', arrow_props={'ec': 'w', 'fc': 'w', 'alpha': 1,'lw': 2})
+            self.ax1.add_artist(north_dir1)
+        except KeyError:
+            pass
 
-        self.cr_UV, = self.ax1.plot(*self.wcs_UV.wcs.crpix, 'r+')
+        self.cr_map, = self.ax1.plot(*self.wcs_map.wcs.crpix, 'r+')
 
         self.ax1.set(xlabel="Right Ascension (J2000)", ylabel="Declination (J2000)", title="Click on selected point of reference.")
 
@@ -426,14 +443,20 @@ class align_maps(object):
         self.ax2 = self.fig.add_subplot(122, projection=self.wcs_other)
         self.ax2.set_facecolor('k')
 
-        vmin, vmax = 0., np.max(other_data[other_data > 0.])
-        im2 = self.ax2.imshow(other_data, vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.)
+        vmin, vmax = 0., np.max(other_data[other_data > 0.]*other_convert)
+        im2 = self.ax2.imshow(other_data*other_convert, vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.)
 
         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': 'w', 'fc': 'w', 'alpha': 1,'lw': 2})
+            self.ax2.add_artist(north_dir2)
+        except KeyError:
+            pass
+
         self.cr_other, = self.ax2.plot(*self.wcs_other.wcs.crpix, 'r+')
 
         self.ax2.set(xlabel="Right Ascension (J2000)", ylabel="Declination (J2000)", title="Click on selected point of reference.")
@@ -445,7 +468,7 @@ class align_maps(object):
         self.breset = Button(self.axreset, 'Leave as is')
      
     def get_aligned_wcs(self):
-        return self.wcs_UV, self.wcs_other
+        return self.wcs_map, self.wcs_other
 
     def onclick_ref(self, event) -> None:
         if self.fig.canvas.manager.toolbar.mode == '':
@@ -453,7 +476,7 @@ class align_maps(object):
                 x = event.xdata
                 y = event.ydata
 
-                self.cr_UV.set(data=[x,y])
+                self.cr_map.set(data=[x,y])
                 self.fig.canvas.draw_idle()
             
             if (event.inaxes is not None) and (event.inaxes == self.ax2):
@@ -464,7 +487,7 @@ class align_maps(object):
                 self.fig.canvas.draw_idle()
     
     def reset_align(self, event):
-        self.wcs_UV.wcs.crpix = WCS(self.Stokes_UV[0].header).wcs.crpix[:2]
+        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]
         self.fig.canvas.draw_idle()
 
@@ -474,9 +497,9 @@ class align_maps(object):
         self.aligned = True
 
     def apply_align(self, event):
-        self.wcs_UV.wcs.crpix = np.array(self.cr_UV.get_data())
+        self.wcs_map.wcs.crpix = np.array(self.cr_map.get_data())
         self.wcs_other.wcs.crpix = np.array(self.cr_other.get_data())
-        self.wcs_other.wcs.crval = self.wcs_UV.wcs.crval
+        self.wcs_other.wcs.crval = self.wcs_map.wcs.crval
         self.fig.canvas.draw_idle()
 
         if self.aligned:
@@ -497,6 +520,7 @@ class align_maps(object):
         plt.show(block=True)
         return self.get_aligned_wcs()
 
+
 class overplot_radio(align_maps):
     """
     Class to overplot maps from different observations.
@@ -555,12 +579,15 @@ class overplot_radio(align_maps):
 
         self.ax.set(xlabel="Right Ascension (J2000)", ylabel="Declination (J2000)", title="HST/FOC UV polarization map of {0:s} overplotted with {1:.2f}GHz map in {2:s}.".format(obj, other_freq*1e-9, other_unit))
 
-        #Display pixel scale
+        #Display pixel scale and North direction
         fontprops = fm.FontProperties(size=16)
         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': 'w', '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)
 
@@ -629,11 +656,14 @@ class overplot_pol(align_maps):
         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
+        #Display pixel scale and North direction
         fontprops = fm.FontProperties(size=16)
         px_size = self.wcs_other.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': 'w', '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))
 
@@ -718,7 +748,7 @@ class crop_map(object):
         #Update WCS and header in new cropped image
         crpix = np.array(self.wcs.wcs.crpix)
         self.wcs_crop = self.wcs.deepcopy()
-        self.wcs_crop.wcs.array_shape = shape
+        self.wcs_crop.array_shape = shape
         if self.crpix_in_RS():
             self.wcs_crop.wcs.crpix -= self.RSextent[::2]
         else: