add overplot and crop classes in lib.plots

.gitignore

@@ -12,3 +12,4 @@ ipython_config.py
 *.fits
 *.npy
 *.par
+data/IC5063_x3nl030/POLARIZATION_COMPARISON/POLARIZATION_COMPARISON-20220329T133819Z-001.zip

src/FOC_reduction.py

@@ -184,15 +184,16 @@ def main():
     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)
 
     ## Step 4:
-    # crop to desired region of interest (roi)
-#    stokescrop = proj_plots.crop_map(deepcopy(stokes_test), deepcopy(data_mask), snrp_cut=snrp_cut, snri_cut=snri_cut)
-#    stokescrop.run()
-#    stokes_crop, data_mask = stokescrop.crop()
-
-    ## Step 5:
     # 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)
 
+    ## Step 5:
+    # crop to desired region of interest (roi)
+    stokescrop = proj_plots.crop_Stokes(deepcopy(Stokes_test))
+    stokescrop.crop()
+    stokescrop.writeto(data_folder+figname+figtype+"_crop.fits")
+    stokes_crop, data_mask = stokescrop.hdul_crop, stokescrop.data_mask
+
     # Plot polarization map (Background is either total Flux, Polarization degree or Polarization degree error).
     if px_scale.lower() not in ['full','integrate']:
         proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, rectangle=None, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype, plots_folder=plots_folder, display=None)

src/lib/plots.py

@@ -9,7 +9,7 @@ prototypes :
         Plots polarization map of polarimetric parameters saved in an HDUList
 """
 
-import copy
+from copy import deepcopy
 import numpy as np
 import matplotlib.pyplot as plt
 from matplotlib.patches import Rectangle
@@ -18,6 +18,7 @@ from matplotlib.transforms import Bbox
 import matplotlib.font_manager as fm
 from mpl_toolkits.axes_grid1.anchored_artists import AnchoredSizeBar, AnchoredDirectionArrows
 from astropy.wcs import WCS
+from astropy.io import fits
 
 
 def princ_angle(ang):
@@ -295,7 +296,7 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
     elif display.lower() in ['s_p','pol_err','pol_deg_err']:
         # Display polarization degree error map
         vmin, vmax = 0., 10.
-        p_err = pol_err.data.copy()
+        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.)
         cbar = plt.colorbar(im, cax=cbar_ax, label=r"$\sigma_P$ [%]")
@@ -496,7 +497,7 @@ class align_maps(object):
         plt.show(block=True)
         return self.get_aligned_wcs()
 
-class overplot_maps(align_maps):
+class overplot_radio(align_maps):
     """
     Class to overplot maps from different observations.
     Inherit from class align_maps in order to get the same WCS on both maps.
@@ -511,8 +512,8 @@ class overplot_maps(align_maps):
         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_unit = self.other_map[0].header['bunit']
         other_convert = 1.
+        other_unit = self.other_map[0].header['bunit']
         if other_unit.lower() == 'jy/beam':
             other_unit = r"mJy/Beam"
             other_convert = 1e3
@@ -554,6 +555,12 @@ class overplot_maps(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
+        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)
+        
         if not(savename is None):
             self.fig2.savefig(savename,bbox_inches='tight',dpi=200)
 
@@ -565,21 +572,27 @@ class overplot_maps(align_maps):
             self.overplot(other_levels=levels, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, savename=savename)
             plt.show(block=True)
 
-
-class crop_map(object):
+class overplot_pol(align_maps):
     """
-    Class to interactively crop I_stokes map to desired Region of Interest
+    Class to overplot maps from different observations.
+    Inherit from class align_maps in order to get the same WCS on both maps.
     """
-    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()
+    def overplot(self, SNRp_cut=3., SNRi_cut=30., savename=None):
+        #Get Data
+        obj = self.Stokes_UV[0].header['targname']
+        stkI = self.Stokes_UV[np.argmax([self.Stokes_UV[i].header['datatype']=='I_stokes' for i in range(len(self.Stokes_UV))])]
+        stk_cov = self.Stokes_UV[np.argmax([self.Stokes_UV[i].header['datatype']=='IQU_cov_matrix' for i in range(len(self.Stokes_UV))])]
+        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']
+        except KeyError:
+            other_convert = 1.
 
         #Compute SNR and apply cuts
         pol.data[pol.data == 0.] = np.nan
@@ -590,72 +603,221 @@ class crop_map(object):
         SNRi[np.isnan(SNRi)] = 0.
         pol.data[SNRi < SNRi_cut] = np.nan
 
-        convert_flux = Stokes[0].header['photflam']
+        plt.rcParams.update({'font.size': 16})
+        self.fig2 = plt.figure(figsize=(15,15))
+        self.ax = self.fig2.add_subplot(111, projection=self.wcs_UV)
+        self.ax.set_facecolor('k')
+        self.fig2.subplots_adjust(hspace=0, wspace=0, right=0.9)
+
+        #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')
+        self.ax.clabel(cont_stkI, inline=True, fontsize=8)
+        
+        self.ax.autoscale(False)
+
+        #Display full size polarization vectors
+        pol.data[np.isfinite(pol.data)] = 1./2.
+        step_vec = 1
+        X, Y = np.meshgrid(np.linspace(0,stkI.data.shape[0],stkI.data.shape[0]), np.linspace(0,stkI.data.shape[1],stkI.data.shape[1]))
+        U, V = pol.data*np.cos(np.pi/2.+pang.data*np.pi/180.), pol.data*np.sin(np.pi/2.+pang.data*np.pi/180.)
+        Q = self.ax.quiver(X[::step_vec,::step_vec],Y[::step_vec,::step_vec],U[::step_vec,::step_vec],V[::step_vec,::step_vec],units='xy',angles='uv',scale=0.5,scale_units='xy',pivot='mid',headwidth=0.,headlength=0.,headaxislength=0.,width=0.1,color='w')
+
+        #Display "other" intensity map
+        vmin, vmax = 0., np.max(other_data[other_data > 0.]*other_convert)
+        im = self.ax.imshow(other_data*other_convert, vmin=vmin, vmax=vmax, transform=self.ax.get_transform(self.wcs_other), 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}$]")
+
+        #Display pixel scale
+        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)
+        
+        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) -> None:
+        self.align()
+        if self.aligned:
+            self.overplot(SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, savename=savename)
+            plt.show(block=True)
+
+
+class crop_map(object):
+    """
+    Class to interactively crop a map to desired Region of Interest
+    """
+    def __init__(self, hdul):
+        #Get data
+        self.hdul = hdul
+        self.header = deepcopy(self.hdul[0].header)
+        self.wcs = WCS(self.header).deepcopy()
+        
+        self.data = deepcopy(self.hdul[0].data)
+        try:
+            convert_flux = self.header['photflam']
+        except KeyError:
+            convert_flux = 1.
 
         #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 = self.fig.add_subplot(111, projection=self.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.]
+        self.ax.plot(*self.wcs.wcs.crpix, 'r+')
+        self.extent = np.array([0.,self.data.shape[0],0., self.data.shape[1]])
+        self.center = np.array(self.data.shape)/2
+        self.RSextent = deepcopy(self.extent)
+        self.RScenter = deepcopy(self.center)
 
-        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.)
+        vmin, vmax = 0., np.max(self.data[self.data > 0.]*convert_flux)
+        im = self.ax.imshow(self.data*convert_flux, vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1., origin='lower')
         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)
+        #Selection button
+        self.axapply = self.fig.add_axes([0.80, 0.01, 0.1, 0.04])
+        self.bapply = Button(self.axapply, 'Apply crop')
+        self.axreset = self.fig.add_axes([0.60, 0.01, 0.1, 0.04])
+        self.breset = Button(self.axreset, 'Reset')
 
-        self.ax.set_title(
-        "Click and drag to crop to desired Region of Interest.\n"
-        "Press 'c' to toggle the selector on and off")
+        self.ax.set_title("Click and drag to crop to desired Region of Interest.")
+    
+    def crpix_in_RS(self):
+        crpix = self.wcs.wcs.crpix
+        x_lim, y_lim = self.RSextent[:2], self.RSextent[2:]
+        if (crpix[0] > x_lim[0] and crpix[0] < x_lim[1]):
+            if (crpix[1] > y_lim[0] and crpix[1] < y_lim[1]):
+                return True
+        return False
+
+    def reset_crop(self, event):
+        self.RSextent = deepcopy(self.extent)
+        self.RScenter = deepcopy(self.center)
+        self.ax.set_xlim(*self.extent[:2])
+        self.ax.set_ylim(*self.extent[2:])
+        self.rect_selector.clear()
+        self.fig.canvas.draw_idle()
 
     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)]
+        self.RSextent = np.array(self.rect_selector.extents)
+        self.RScenter = np.array(self.rect_selector.center)
 
-        # 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 apply_crop(self, event):
+        vertex = self.RSextent.astype(int)
+        shape = vertex[1::2] - vertex[0::2]
+        #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
+        if self.crpix_in_RS():
+            self.wcs_crop.wcs.crpix -= self.RSextent[::2]
+        else:
+            self.wcs_crop.wcs.crval = self.wcs.wcs_pix2world([self.RScenter],1)[0]
+            self.wcs_crop.wcs.crpix = self.RScenter-self.RSextent[::2]
+         
+        # Crop dataset
+        self.data_crop = self.data[vertex[2]:vertex[3], vertex[0]:vertex[1]]
 
-    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)
+        #Write cropped map to new HDUList
+        self.header_crop = deepcopy(self.header)
+        self.header_crop.update(self.wcs_crop.to_header())
+        self.hdul_crop = fits.HDUList([fits.PrimaryHDU(self.data_crop,self.header_crop)])
+
+        try:
+            convert_flux = self.header_crop['photflam']
+        except KeyError:
+            convert_flux = 1.
+
+        self.fig.clear()
+        self.ax = self.fig.add_subplot(111,projection=self.wcs_crop)
+        vmin, vmax = 0., np.max(self.data_crop[self.data_crop > 0.]*convert_flux)
+        im = self.ax.imshow(self.data_crop*convert_flux, vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1., origin='lower')
+        self.ax.plot(*self.wcs_crop.wcs.crpix, 'r+')
+        cbar_ax = self.fig.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}$]")
+        xlim, ylim = self.RSextent[1::2]-self.RSextent[0::2]
+        self.ax.set_xlim(0,xlim)
+        self.ax.set_ylim(0,ylim)
+        self.rect_selector.clear()
+
+        self.fig.canvas.draw_idle()
+
+    def crop(self) -> None:
+        if self.fig.canvas.manager.toolbar.mode == '':
+            self.rect_selector = RectangleSelector(self.ax, self.onselect_crop,
+                    drawtype='box', button=[1], interactive=True)
+        self.bapply.on_clicked(self.apply_crop)
+        self.breset.on_clicked(self.reset_crop)
         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
+    def writeto(self, filename):
+        self.hdul_crop.writeto(filename,overwrite=True)
 
-        return Stokes_crop, data_mask
+class crop_Stokes(crop_map):
+    """
+    Class to interactively crop a polarization map to desired Region of Interest.
+    Inherit from crop_map.
+    """
+    def apply_crop(self,event):
+        """
+        Redefine apply_crop method for the Stokes HDUList.
+        """
+        self.hdul_crop = deepcopy(self.hdul)
+        vertex = self.RSextent.astype(int)
+        shape = vertex[1::2] - vertex[0::2]
+        #Update WCS and header in new cropped image
+        crpix = np.array(self.wcs.wcs.crpix)
+        self.wcs_crop = self.wcs.deepcopy()
+        self.wcs_crop.array_shape = shape
+        if self.crpix_in_RS():
+            self.wcs_crop.wcs.crpix -= self.RSextent[::2]
+        else:
+            self.wcs_crop.wcs.crval = self.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':
+                stokes_cov = np.zeros((3,3,shape[1],shape[0]))
+                for i in range(3):
+                    for j in range(3):
+                        stokes_cov[i,j] = dataset.data[i,j][vertex[2]:vertex[3], vertex[0]:vertex[1]]
+                dataset.data = stokes_cov
+            else:
+                dataset.data = dataset.data[vertex[2]:vertex[3], vertex[0]:vertex[1]]
+            dataset.header.update(self.wcs_crop.to_header())
+        self.data_mask = self.hdul_crop[-1].data
 
+        try:
+            convert_flux = self.hdul_crop[0].header['photflam']
+        except KeyError:
+            convert_flux = 1.
+
+        data_crop = self.hdul_crop[0].data
+        self.fig.clear()
+        self.ax = self.fig.add_subplot(111,projection=self.wcs_crop)
+        vmin, vmax = 0., np.max(data_crop[data_crop > 0.]*convert_flux)
+        im = self.ax.imshow(data_crop*convert_flux, vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1., origin='lower')
+        self.ax.plot(*self.wcs_crop.wcs.crpix, 'r+')
+        cbar_ax = self.fig.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}$]")
+        xlim, ylim = self.RSextent[1::2]-self.RSextent[0::2]
+        self.ax.set_xlim(0,xlim)
+        self.ax.set_ylim(0,ylim)
+        self.rect_selector.clear()
+
+        self.fig.canvas.draw_idle()
+    
+    @property
+    def data_mask(self):
+        return self.data_mask

src/lib/reduction.py

@@ -1389,7 +1389,8 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers,
             new_Stokes_cov[i,j] = sc_rotate(new_Stokes_cov[i,j], ang,
                     reshape=False, cval=0.)
         new_Stokes_cov[i,i] = np.abs(new_Stokes_cov[i,i])
-    center = np.array(new_I_stokes.shape)/2
+    old_center = np.array(I_stokes.shape)/2
+    new_center = np.array(new_I_stokes.shape)/2
 
     #Update headers to new angle
     new_headers = []
@@ -1417,7 +1418,7 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers,
         elif new_wcs.wcs.has_pc():      # PC matrix + CDELT
             newpc = np.dot(mrot, new_wcs.wcs.get_pc())
             new_wcs.wcs.pc = newpc
-        new_wcs.wcs.crpix = np.dot(mrot, new_wcs.wcs.crpix - center) + center
+        new_wcs.wcs.crpix = np.dot(mrot, new_wcs.wcs.crpix - old_center) + new_center
         new_wcs.wcs.set()
         new_header.update(new_wcs.to_header())
 

src/lib/test_overplot.py

@@ -1,7 +1,7 @@
 #!/usr/bin/python3
 from astropy.io import fits
 import numpy as np
-from plots import overplot_maps
+from plots import overplot_radio
 
 Stokes_UV = fits.open("../../data/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_pol.fits")
 Stokes_18GHz = fits.open("../../data/IC5063_x3nl030/radio/IC5063.18GHz.fits")
@@ -11,10 +11,10 @@ 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=10.0, SNRi_cut=100.0, savename='../../plots/IC5063_x3nl030/18GHz_overplot.png')
+A = overplot_radio(Stokes_UV, Stokes_18GHz)
+A.plot(levels=levels18GHz, SNRp_cut=10.0, SNRi_cut=100.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 = levelsMorganti*0.46*1e-3
-B = overplot_maps(Stokes_UV, Stokes_24GHz)
-B.plot(levels=levels24GHz, SNRp_cut=10.0, SNRi_cut=100.0, savename='../../plots/IC5063_x3nl030/24GHz_overplot.png')
+B = overplot_radio(Stokes_UV, Stokes_24GHz)
+B.plot(levels=levels24GHz, SNRp_cut=10.0, SNRi_cut=100.0, savename='../../plots/IC5063_x3nl030/24GHz_overplot_forced.png')