overplot IC5063 with Radio and IR

src/lib/plots.py

@@ -449,21 +449,21 @@ class align_maps(object):
         self.map = map1
         self.other_map = other_map
 
-        self.wcs_map = WCS(self.map[0]).deepcopy()
-        if self.wcs_map.naxis == 4:
-            self.wcs_map = WCS(self.map[0],naxis=[1,2]).deepcopy()
-            self.map[0].data = self.map[0].data[0,0]
-        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_map = deepcopy(WCS(self.map[0])).celestial
+#        if self.wcs_map.naxis == 4:
+#            self.wcs_map = WCS(self.map[0],naxis=[1,2]).deepcopy()
+#            self.map[0].data = self.map[0].data[0,0]
+#        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()
-            self.other_map[0].data = self.other_map[0].data[0,0]
-        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]
+        self.wcs_other = deepcopy(WCS(self.other_map[0])).celestial
+#        if self.wcs_other.naxis == 4:
+#            self.wcs_other = WCS(self.other_map[0],naxis=[1,2]).deepcopy()
+#            self.other_map[0].data = self.other_map[0].data[0,0]
+#        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']
@@ -773,7 +773,7 @@ class overplot_pol(align_maps):
         self.fig2.canvas.draw()
 
     def plot(self, SNRp_cut=3., SNRi_cut=30., savename=None, **kwargs) -> None:
-        while not self.aligned():
+        while not self.aligned:
             self.align()
         self.overplot(SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, savename=savename, **kwargs)
         plt.show(block=True)
@@ -1739,7 +1739,8 @@ class pol_map(object):
             label = r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]"
         elif self.display_selection.lower() in ['pol_flux']:
             self.data = self.I*self.convert_flux*self.P
-            vmin, vmax = 0., np.max(self.data[self.data > 0.])
+            vmin, vmax = np.min(self.I[self.cut]*self.convert_flux)/10., np.max(self.I[self.data > 0.]*self.convert_flux)
+            norm = LogNorm(vmin, vmax)
             label = r"$F_{\lambda} \cdot P$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]"
         elif self.display_selection.lower() in ['pol_deg']:
             self.data = self.P*100.
@@ -1894,12 +1895,14 @@ class pol_map(object):
                 ax = self.ax
             if hasattr(self, 'an_int'):
                 self.an_int.remove()
+            #self.an_int = ax.annotate(r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(self.pivot_wav,sci_not(I_reg*self.convert_flux,I_reg_err*self.convert_flux,2))+"\n"+r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_reg*100.,np.ceil(P_reg_err*1000.)/10.)+"\n"+r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_reg,np.ceil(PA_reg_err*10.)/10.), color='white', fontsize=12, xy=(0.01, 0.93), xycoords='axes fraction')
             self.an_int = ax.annotate(r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(self.pivot_wav,sci_not(I_reg*self.convert_flux,I_reg_err*self.convert_flux,2))+"\n"+r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_reg*100.,np.ceil(P_reg_err*1000.)/10.)+"\n"+r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_reg,np.ceil(PA_reg_err*10.)/10.)+"\n"+r"$P^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_cut*100.,np.ceil(P_cut_err*1000.)/10.)+"\n"+r"$\theta_{{P}}^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_cut,np.ceil(PA_cut_err*10.)/10.), color='white', fontsize=12, xy=(0.01, 0.85), xycoords='axes fraction')
             if not self.region is None:
                 self.cont = ax.contour(self.region.astype(float),levels=[0.5], colors='white', linewidths=0.8)
             fig.canvas.draw_idle()
             return self.an_int
         else:
+            #ax.annotate(r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(self.pivot_wav,sci_not(I_reg*self.convert_flux,I_reg_err*self.convert_flux,2))+"\n"+r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_reg*100.,np.ceil(P_reg_err*1000.)/10.)+"\n"+r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_reg,np.ceil(PA_reg_err*10.)/10.), color='white', fontsize=12, xy=(0.01, 0.94), xycoords='axes fraction')
             ax.annotate(r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(self.pivot_wav,sci_not(I_reg*self.convert_flux,I_reg_err*self.convert_flux,2))+"\n"+r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_reg*100.,np.ceil(P_reg_err*1000.)/10.)+"\n"+r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_reg,np.ceil(PA_reg_err*10.)/10.)+"\n"+r"$P^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_cut*100.,np.ceil(P_cut_err*1000.)/10.)+"\n"+r"$\theta_{{P}}^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_cut,np.ceil(PA_cut_err*10.)/10.), color='white', fontsize=12, xy=(0.01, 0.90), xycoords='axes fraction')
             if not self.region is None:
                 ax.contour(self.region.astype(float),levels=[0.5], colors='white', linewidths=0.8)

src/lib/reduction.py

@@ -246,10 +246,10 @@ def crop_array(data_array, headers, error_array=None, data_mask=None, step=5,
     Returns:
     cropped_array : numpy.ndarray
         Array containing the observationnal data homogeneously cropped.
-    headers : header list
-        Updated headers associated with the images in data_array.
     cropped_error : numpy.ndarray
         Array containing the error on the observationnal data homogeneously cropped.
+    headers : header list
+        Updated headers associated with the images in data_array.
     """
     if error_array is None:
         error_array = np.zeros(data_array.shape)
@@ -284,7 +284,7 @@ def crop_array(data_array, headers, error_array=None, data_mask=None, step=5,
         crop_error_array[i] = error_array[i][v_array[0]:v_array[1],v_array[2]:v_array[3]]
         #Update CRPIX value in the associated header
         curr_wcs = deepcopy(WCS(crop_headers[i]))
-        curr_wcs.wcs.crpix = curr_wcs.wcs.crpix - np.array([v_array[2], v_array[0]])
+        curr_wcs.wcs.crpix[:2] = curr_wcs.wcs.crpix[:2] - np.array([v_array[2], v_array[0]])
         crop_headers[i].update(curr_wcs.to_header())
         crop_headers[i]['naxis1'], crop_headers[i]['naxis2'] = crop_array[i].shape
 
@@ -1606,11 +1606,11 @@ def rotate_data(data_array, error_array, data_mask, headers, ang):
         Updated array containing the rotated images.
     new_error_array : numpy.ndarray
         Updated array containing the rotated errors.
+    new_data_mask : numpy.ndarray
+        Updated 2D boolean array delimiting the data to work on.
     new_headers : header list
         Updated list of headers corresponding to the reduced images accounting
         for the new orientation angle.
-    new_data_mask : numpy.ndarray
-        Updated 2D boolean array delimiting the data to work on.
     """
     #Rotate I_stokes, Q_stokes, U_stokes using rotation matrix
     alpha = ang*np.pi/180.
@@ -1649,8 +1649,8 @@ def rotate_data(data_array, error_array, data_mask, headers, ang):
 
         new_wcs = WCS(header).deepcopy()
 
-        new_wcs.wcs.pc = np.dot(mrot, new_wcs.wcs.pc)
-        new_wcs.wcs.crpix = np.dot(mrot, new_wcs.wcs.crpix - old_center[::-1]) + new_center[::-1]
+        new_wcs.wcs.pc[:2,:2] = np.dot(mrot, new_wcs.wcs.pc[:2,:2])
+        new_wcs.wcs.crpix[:2] = np.dot(mrot, new_wcs.wcs.crpix[:2] - old_center[::-1]) + new_center[::-1]
         new_wcs.wcs.set()
         for key, val in new_wcs.to_header().items():
             new_header[key] = val

src/overplot.py

@@ -26,47 +26,47 @@ from matplotlib.colors import LogNorm
 #levels24GHz = levelsMorganti*0.46*1e-3
 #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')
-#
+
 #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.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.]))
 #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')
-#
+
 #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.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.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 = overplot_pol(Stokes_UV, Stokes_IR, cmap='inferno')
 #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')
 
-data_folder1 = "../data/M87/POS1/"
-plots_folder1 = "../plots/M87/POS1/"
-basename1 = "M87_020_log_core"
-M87_1_95 = fits.open(data_folder1+"M87_POS1_1995_FOC_combine_FWHM020.fits")
-M87_1_96 = fits.open(data_folder1+"M87_POS1_1996_FOC_combine_FWHM020.fits")
-M87_1_97 = fits.open(data_folder1+"M87_POS1_1997_FOC_combine_FWHM020.fits")
-M87_1_98 = fits.open(data_folder1+"M87_POS1_1998_FOC_combine_FWHM020.fits")
-M87_1_99 = fits.open(data_folder1+"M87_POS1_1999_FOC_combine_FWHM020.fits")
+#data_folder1 = "../data/M87/POS1/"
+#plots_folder1 = "../plots/M87/POS1/"
+#basename1 = "M87_020_log"
+#M87_1_95 = fits.open(data_folder1+"M87_POS1_1995_FOC_combine_FWHM020.fits")
+#M87_1_96 = fits.open(data_folder1+"M87_POS1_1996_FOC_combine_FWHM020.fits")
+#M87_1_97 = fits.open(data_folder1+"M87_POS1_1997_FOC_combine_FWHM020.fits")
+#M87_1_98 = fits.open(data_folder1+"M87_POS1_1998_FOC_combine_FWHM020.fits")
+#M87_1_99 = fits.open(data_folder1+"M87_POS1_1999_FOC_combine_FWHM020.fits")
 
-H = align_pol(np.array([M87_1_95,M87_1_96,M87_1_97,M87_1_98,M87_1_99]), norm=LogNorm())
-H.plot(SNRp_cut=5.0, SNRi_cut=50.0, savename=plots_folder1+'animated_loop/'+basename1, norm=LogNorm())
-command("convert -delay 50 -loop 0 {0:s}animated_loop/{1:s}*.png {0:s}animated_loop/{1:s}.gif".format(plots_folder1, basename1))
+#H = align_pol(np.array([M87_1_95,M87_1_96,M87_1_97,M87_1_98,M87_1_99]), norm=LogNorm())
+#H.plot(SNRp_cut=5.0, SNRi_cut=50.0, savename=plots_folder1+'animated_loop/'+basename1, norm=LogNorm())
+#command("convert -delay 50 -loop 0 {0:s}animated_loop/{1:s}*.png {0:s}animated_loop/{1:s}.gif".format(plots_folder1, basename1))
 
-data_folder3 = "../data/M87/POS3/"
-plots_folder3 = "../plots/M87/POS3/"
-basename3 = "M87_020_log_star"
-M87_3_95 = fits.open(data_folder3+"M87_POS3_1995_FOC_combine_FWHM020.fits")
-M87_3_96 = fits.open(data_folder3+"M87_POS3_1996_FOC_combine_FWHM020.fits")
-M87_3_97 = fits.open(data_folder3+"M87_POS3_1997_FOC_combine_FWHM020.fits")
-M87_3_98 = fits.open(data_folder3+"M87_POS3_1998_FOC_combine_FWHM020.fits")
-M87_3_99 = fits.open(data_folder3+"M87_POS3_1999_FOC_combine_FWHM020.fits")
+#data_folder3 = "../data/M87/POS3/"
+#plots_folder3 = "../plots/M87/POS3/"
+#basename3 = "M87_020_log"
+#M87_3_95 = fits.open(data_folder3+"M87_POS3_1995_FOC_combine_FWHM020.fits")
+#M87_3_96 = fits.open(data_folder3+"M87_POS3_1996_FOC_combine_FWHM020.fits")
+#M87_3_97 = fits.open(data_folder3+"M87_POS3_1997_FOC_combine_FWHM020.fits")
+#M87_3_98 = fits.open(data_folder3+"M87_POS3_1998_FOC_combine_FWHM020.fits")
+#M87_3_99 = fits.open(data_folder3+"M87_POS3_1999_FOC_combine_FWHM020.fits")
 
-I = align_pol(np.array([M87_3_95,M87_3_96,M87_3_97,M87_3_98,M87_3_99]), norm=LogNorm())
-I.plot(SNRp_cut=5.0, SNRi_cut=50.0, savename=plots_folder3+'animated_loop/'+basename3, norm=LogNorm())
-command("convert -delay 20 -loop 0 {0:s}animated_loop/{1:s}*.png {0:s}animated_loop/{1:s}.gif".format(plots_folder3, basename3))
+#I = align_pol(np.array([M87_3_95,M87_3_96,M87_3_97,M87_3_98,M87_3_99]), norm=LogNorm())
+#I.plot(SNRp_cut=5.0, SNRi_cut=50.0, savename=plots_folder3+'animated_loop/'+basename3, norm=LogNorm())
+#command("convert -delay 20 -loop 0 {0:s}animated_loop/{1:s}*.png {0:s}animated_loop/{1:s}.gif".format(plots_folder3, basename3))