diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020.png index b4b8330..ca027ef 100644 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020.png and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_IQU.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_IQU.png index e48eff7..26f643a 100644 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_IQU.png and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_IQU.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_I_err.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_I_err.png index d8611f8..56a2684 100644 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_I_err.png and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_I_err.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P.png index 61ee811..89f8c62 100644 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P.png and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_err.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_err.png index ecb1251..c73bea0 100644 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_err.png and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_err.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_flux.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_flux.png index f730101..20d6e8b 100644 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_flux.png and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_flux.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRi.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRi.png index f4ad036..d99e87e 100644 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRi.png and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRi.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRp.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRp.png index 3432db0..bf0205e 100644 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRp.png and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRp.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020.png index d3b9beb..2f2efc2 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_IQU.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_IQU.png index 22aaa7e..1022652 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_IQU.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_IQU.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_I_err.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_I_err.png index 8f471a1..5e0699f 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_I_err.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_I_err.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P.png index 1eaefcd..000d933 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P_err.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P_err.png index 94f4e05..e2cd4cd 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P_err.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P_err.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P_flux.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P_flux.png index 26669d9..b4220e6 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P_flux.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P_flux.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_SNRi.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_SNRi.png index d47e24c..a5badbb 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_SNRi.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_SNRi.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_SNRp.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_SNRp.png index 17d6590..bcc043d 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_SNRp.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_SNRp.png differ diff --git a/src/FOC_reduction.py b/src/FOC_reduction.py index d87e35d..993d08e 100755 --- a/src/FOC_reduction.py +++ b/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) diff --git a/src/lib/fits.py b/src/lib/fits.py index 68c3c76..54b96a1 100755 --- a/src/lib/fits.py +++ b/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'], diff --git a/src/lib/plots.py b/src/lib/plots.py index 38b8dd3..db91ffe 100755 --- a/src/lib/plots.py +++ b/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: