diff --git a/.gitignore b/.gitignore index 643fa47..e530053 100755 --- a/.gitignore +++ b/.gitignore @@ -12,3 +12,4 @@ ipython_config.py *.fits *.npy *.par +data/IC5063_x3nl030/POLARIZATION_COMPARISON/POLARIZATION_COMPARISON-20220329T133819Z-001.zip diff --git a/plots/IC5063_x3nl030/18GHz_overplot.png b/plots/IC5063_x3nl030/18GHz_overplot.png index 6832971..7d3face 100644 Binary files a/plots/IC5063_x3nl030/18GHz_overplot.png and b/plots/IC5063_x3nl030/18GHz_overplot.png differ diff --git a/plots/IC5063_x3nl030/18GHz_overplot_forced.png b/plots/IC5063_x3nl030/18GHz_overplot_forced.png index 1a6e4f1..d8a5b50 100644 Binary files a/plots/IC5063_x3nl030/18GHz_overplot_forced.png and b/plots/IC5063_x3nl030/18GHz_overplot_forced.png differ diff --git a/plots/IC5063_x3nl030/24GHz_overplot.png b/plots/IC5063_x3nl030/24GHz_overplot.png index 3e289d4..ee0b941 100644 Binary files a/plots/IC5063_x3nl030/24GHz_overplot.png and b/plots/IC5063_x3nl030/24GHz_overplot.png differ diff --git a/plots/IC5063_x3nl030/24GHz_overplot_forced.png b/plots/IC5063_x3nl030/24GHz_overplot_forced.png index d6ad8bd..98b6211 100644 Binary files a/plots/IC5063_x3nl030/24GHz_overplot_forced.png and b/plots/IC5063_x3nl030/24GHz_overplot_forced.png differ diff --git a/plots/IC5063_x3nl030/S2_overplot_forced.png b/plots/IC5063_x3nl030/S2_overplot_forced.png new file mode 100644 index 0000000..4b4019b Binary files /dev/null and b/plots/IC5063_x3nl030/S2_overplot_forced.png differ diff --git a/plots/IC5063_x3nl030/S2_sqrt_overplot_forced.png b/plots/IC5063_x3nl030/S2_sqrt_overplot_forced.png new file mode 100644 index 0000000..1833669 Binary files /dev/null and b/plots/IC5063_x3nl030/S2_sqrt_overplot_forced.png differ diff --git a/plots/IC5063_x3nl030/axes_24GHz_overplot.png b/plots/IC5063_x3nl030/axes_24GHz_overplot.png new file mode 100644 index 0000000..77fdc7b Binary files /dev/null and b/plots/IC5063_x3nl030/axes_24GHz_overplot.png differ diff --git a/src/FOC_reduction.py b/src/FOC_reduction.py index ba805d9..e41ed67 100755 --- a/src/FOC_reduction.py +++ b/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) diff --git a/src/lib/plots.py b/src/lib/plots.py index 1799a95..41aa0af 100755 --- a/src/lib/plots.py +++ b/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 \ No newline at end of file diff --git a/src/lib/reduction.py b/src/lib/reduction.py index 6e4f38b..88c2adc 100755 --- a/src/lib/reduction.py +++ b/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()) diff --git a/src/lib/test_overplot.py b/src/lib/test_overplot.py index 508e5e1..0e23a0b 100755 --- a/src/lib/test_overplot.py +++ b/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')