diff --git a/src/FOC_reduction.py b/src/FOC_reduction.py index 30e8874..e6e4e38 100755 --- a/src/FOC_reduction.py +++ b/src/FOC_reduction.py @@ -12,9 +12,10 @@ import lib.reduction as proj_red #Functions used in reduction pipeline import lib.plots as proj_plots #Functions for plotting data from lib.deconvolve import from_file_psf from lib.query import retrieve_products, path_exists, system +from matplotlib.colors import LogNorm -def main(target=None, proposal_id=None, infiles=None, output_dir="./data"): +def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=0, interactive=0): ## Reduction parameters # Deconvolution deconvolve = False @@ -53,18 +54,18 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data"): smoothing_scale = 'arcsec' #pixel or arcsec # Rotation - rotate_stokes = True rotate_data = False #rotation to North convention can give erroneous results + rotate_stokes = True # Final crop - crop = False #Crop to desired ROI - final_display = True #Whether to display all polarization map outputs + #crop = False #Crop to desired ROI + #interactive = False #Whether to output to intercative analysis tool # Polarization map output SNRp_cut = 3. #P measurments with SNR>3 SNRi_cut = 30. #I measurments with SNR>30, which implies an uncertainty in P of 4.7%. - flux_lim = [5e-19,5e-14] #lowest and highest flux displayed on plot, defaults to bkg and maximum in cut if None - vec_scale = 2.0 + flux_lim = None #lowest and highest flux displayed on plot, defaults to bkg and maximum in cut if None + vec_scale = 5 step_vec = 1 #plot all vectors in the array. if step_vec = 2, then every other vector will be plotted # if step_vec = 0 then all vectors are displayed at full length @@ -168,21 +169,21 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data"): ## Step 5: # crop to desired region of interest (roi) if crop: - figtype += "crop" - stokescrop = proj_plots.crop_Stokes(deepcopy(Stokes_test)) + figtype += "_crop" + stokescrop = proj_plots.crop_Stokes(deepcopy(Stokes_test),norm=LogNorm()) stokescrop.crop() stokescrop.writeto("/".join([data_folder,"_".join([figname,figtype+".fits"])])) - Stokes_test, data_mask = stokescrop.hdul_crop, stokescrop.data_mask + Stokes_test, data_mask, headers = stokescrop.hdul_crop, stokescrop.data_mask, [dataset.header for dataset in stokescrop.hdul_crop] print("F_int({0:.0f} Angs) = ({1} ± {2})e{3} ergs.cm^-2.s^-1.Angs^-1".format(headers[0]['photplam'],*proj_plots.sci_not(Stokes_test[0].data[data_mask].sum()*headers[0]['photflam'],np.sqrt(Stokes_test[3].data[0,0][data_mask].sum())*headers[0]['photflam'],2,out=int))) print("P_int = {0:.1f} ± {1:.1f} %".format(headers[0]['p_int']*100.,np.ceil(headers[0]['p_int_err']*1000.)/10.)) - print("PA_int = {0:.1f} ± {1:.1f} °".format(headers[0]['pa_int'],np.ceil(headers[0]['pa_int_err']*10.)/10.)) + print("PA_int = {0:.1f} ±t {1:.1f} °".format(headers[0]['pa_int'],np.ceil(headers[0]['pa_int_err']*10.)/10.)) # Background values print("F_bkg({0:.0f} Angs) = ({1} ± {2})e{3} ergs.cm^-2.s^-1.Angs^-1".format(headers[0]['photplam'],*proj_plots.sci_not(I_bkg[0,0]*headers[0]['photflam'],np.sqrt(S_cov_bkg[0,0][0,0])*headers[0]['photflam'],2,out=int))) print("P_bkg = {0:.1f} ± {1:.1f} %".format(debiased_P_bkg[0,0]*100.,np.ceil(s_P_bkg[0,0]*1000.)/10.)) print("PA_bkg = {0:.1f} ± {1:.1f} °".format(PA_bkg[0,0],np.ceil(s_PA_bkg[0,0]*10.)/10.)) # Plot polarization map (Background is either total Flux, Polarization degree or Polarization degree error). - if px_scale.lower() not in ['full','integrate'] and final_display: + if px_scale.lower() not in ['full','integrate'] and not interactive: proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype]), plots_folder=plots_folder) proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"I"]), plots_folder=plots_folder, display='Intensity') proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"P_flux"]), plots_folder=plots_folder, display='Pol_Flux') @@ -192,7 +193,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data"): proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"P_err"]), plots_folder=plots_folder, display='Pol_deg_err') proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"SNRi"]), plots_folder=plots_folder, display='SNRi') proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"SNRp"]), plots_folder=plots_folder, display='SNRp') - elif final_display: + elif not interactive: proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, savename="_".join([figname,figtype]), plots_folder=plots_folder, display='integrate') elif px_scale.lower() not in ['full', 'integrate']: pol_map = proj_plots.pol_map(Stokes_test, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim) @@ -212,6 +213,10 @@ if __name__ == "__main__": help='the full or relative path to the data products', default=None) parser.add_argument('-o','--output_dir', metavar='directory_path', required=False, help='output directory path for the data products', type=str, default="./data") + parser.add_argument('-c','--crop', metavar='crop_boolean', required=False, + help='whether to crop the analysis region', type=int, default=0) + parser.add_argument('-i','--interactive', metavar='interactive_boolean', required=False, + help='whether to output to the interactive analysis tool', type=int, default=0) args = parser.parse_args() - exitcode = main(target=args.target, proposal_id=args.proposal_id, infiles=args.files, output_dir=args.output_dir) + exitcode = main(target=args.target, proposal_id=args.proposal_id, infiles=args.files, output_dir=args.output_dir, crop=args.crop, interactive=args.interactive) print("Finished with ExitCode: ",exitcode) diff --git a/src/lib/plots.py b/src/lib/plots.py index 8b4b8a0..4811b13 100755 --- a/src/lib/plots.py +++ b/src/lib/plots.py @@ -1068,7 +1068,7 @@ class crop_map(object): """ Class to interactively crop a map to desired Region of Interest """ - def __init__(self, hdul, fig=None, ax=None): + def __init__(self, hdul, fig=None, ax=None, **kwargs): #Get data self.cropped=False self.hdul = hdul @@ -1080,10 +1080,13 @@ class crop_map(object): self.convert_flux = self.header['photflam'] except KeyError: self.convert_flux = 1. + try: + self.kwargs = kwargs + except AttributeError: + self.kwargs = {} #Plot the map plt.rcParams.update({'font.size': 12}) - plt.ioff() if fig is None: self.fig = plt.figure(figsize=(15,15)) self.fig.suptitle("Click and drag to crop to desired Region of Interest.") @@ -1104,28 +1107,36 @@ class crop_map(object): self.rect_selector = RectangleSelector(self.ax, self.onselect_crop, button=[1]) self.embedded = True - self.display() - plt.ion() + self.display(self.data, self.wcs, self.convert_flux, **self.kwargs) 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) - plt.show() - def display(self, data=None, wcs=None, convert_flux=None): + def display(self, data=None, wcs=None, convert_flux=None, **kwargs): if data is None: data = self.data if wcs is None: wcs = self.wcs if convert_flux is None: convert_flux = self.convert_flux + if kwargs is None: + kwargs = self.kwargs + else: + kwargs = {**self.kwargs, **kwargs} - vmin, vmax = 0., np.max(data[data > 0.]*convert_flux) + vmin, vmax = np.min(data[data > 0.]*convert_flux), np.max(data[data > 0.]*convert_flux) + for key, value in [["cmap",[["cmap","inferno"]]], ["origin",[["origin","lower"]]], ["aspect",[["aspect","equal"]]], ["alpha",[["alpha",self.mask_alpha]]], ["norm",[["vmin",vmin],["vmax",vmax]]]]: + try: + test = kwargs[key] + except KeyError: + for key_i, val_i in value: + kwargs[key_i] = val_i if hasattr(self, 'im'): self.im.remove() - self.im = self.ax.imshow(data*convert_flux, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=self.mask_alpha, origin='lower') + self.im = self.ax.imshow(data*convert_flux, **kwargs) if hasattr(self, 'cr'): self.cr[0].set_data(*wcs.wcs.crpix) else: @@ -1199,14 +1210,9 @@ class crop_map(object): 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.rect_selector.clear() self.ax.reset_wcs(self.wcs_crop) - self.display(data=self.data_crop, wcs=self.wcs_crop, convert_flux=convert_flux) + self.display(data=self.data_crop, wcs=self.wcs_crop) xlim, ylim = self.RSextent[1::2]-self.RSextent[0::2] self.ax.set_xlim(0,xlim) @@ -1215,6 +1221,7 @@ class crop_map(object): if self.fig.canvas.manager.toolbar.mode == '': self.rect_selector = RectangleSelector(self.ax, self.onselect_crop, button=[1]) + self.fig.canvas.draw_idle() def on_close(self, event) -> None: @@ -1283,16 +1290,11 @@ class crop_Stokes(crop_map): dataset.data = deepcopy(dataset.data[vertex[2]:vertex[3], vertex[0]:vertex[1]]) dataset.header.update(self.wcs_crop.to_header()) - try: - convert_flux = self.hdul_crop[0].header['photflam'] - except KeyError: - convert_flux = 1. - self.data_crop = self.hdul_crop[0].data self.rect_selector.clear() if not self.embedded: self.ax.reset_wcs(self.wcs_crop) - self.display(data=self.data_crop, wcs=self.wcs_crop, convert_flux=convert_flux) + self.display(data=self.data_crop, wcs=self.wcs_crop) xlim, ylim = self.RSextent[1::2]-self.RSextent[0::2] self.ax.set_xlim(0,xlim) @@ -1303,11 +1305,34 @@ class crop_Stokes(crop_map): if self.fig.canvas.manager.toolbar.mode == '': self.rect_selector = RectangleSelector(self.ax, self.onselect_crop, button=[1]) + # Update integrated values + mask = np.logical_and(self.hdul_crop[-1].data.astype(bool), self.hdul_crop[0].data >0) + I_diluted = self.hdul_crop[0].data[mask].sum() + Q_diluted = self.hdul_crop[1].data[mask].sum() + U_diluted = self.hdul_crop[2].data[mask].sum() + I_diluted_err = np.sqrt(np.sum(self.hdul_crop[3].data[0,0][mask])) + Q_diluted_err = np.sqrt(np.sum(self.hdul_crop[3].data[1,1][mask])) + U_diluted_err = np.sqrt(np.sum(self.hdul_crop[3].data[2,2][mask])) + IQ_diluted_err = np.sqrt(np.sum(self.hdul_crop[3].data[0,1][mask]**2)) + IU_diluted_err = np.sqrt(np.sum(self.hdul_crop[3].data[0,2][mask]**2)) + QU_diluted_err = np.sqrt(np.sum(self.hdul_crop[3].data[1,2][mask]**2)) + + P_diluted = np.sqrt(Q_diluted**2+U_diluted**2)/I_diluted + P_diluted_err = (1./I_diluted)*np.sqrt((Q_diluted**2*Q_diluted_err**2 + U_diluted**2*U_diluted_err**2 + 2.*Q_diluted*U_diluted*QU_diluted_err)/(Q_diluted**2 + U_diluted**2) + ((Q_diluted/I_diluted)**2 + (U_diluted/I_diluted)**2)*I_diluted_err**2 - 2.*(Q_diluted/I_diluted)*IQ_diluted_err - 2.*(U_diluted/I_diluted)*IU_diluted_err) + + PA_diluted = princ_angle((90./np.pi)*np.arctan2(U_diluted,Q_diluted)) + PA_diluted_err = (90./(np.pi*(Q_diluted**2 + U_diluted**2)))*np.sqrt(U_diluted**2*Q_diluted_err**2 + Q_diluted**2*U_diluted_err**2 - 2.*Q_diluted*U_diluted*QU_diluted_err) + + for dataset in self.hdul_crop: + dataset.header['P_int'] = (P_diluted, 'Integrated polarization degree') + dataset.header['P_int_err'] = (np.ceil(P_diluted_err*1000.)/1000., 'Integrated polarization degree error') + dataset.header['PA_int'] = (PA_diluted, 'Integrated polarization angle') + dataset.header['PA_int_err'] = (np.ceil(PA_diluted_err*10.)/10., 'Integrated polarization angle error') self.fig.canvas.draw_idle() @property def data_mask(self): - return self.hdul_crop[-1].data + return self.hdul_crop[-1].data.astype(int) class image_lasso_selector(object): diff --git a/src/overplot.py b/src/overplot_IC5063.py similarity index 100% rename from src/overplot.py rename to src/overplot_IC5063.py diff --git a/src/overplot_MRK463E.py b/src/overplot_MRK463E.py new file mode 100755 index 0000000..133aa7e --- /dev/null +++ b/src/overplot_MRK463E.py @@ -0,0 +1,25 @@ +#!/usr/bin/python3 +from os import system as command +from astropy.io import fits +import numpy as np +from copy import deepcopy +from lib.plots import overplot_chandra, overplot_pol, align_pol +from matplotlib.colors import LogNorm + +Stokes_UV = fits.open("./data/MRK463E/5960/MRK463E_FOC_b0.05arcsec_c0.10arcsec.fits") +Stokes_IR = fits.open("./data/MRK463E/WFPC2/IR_rot_crop.fits") +Stokes_Xr = fits.open("./data/MRK463E/Chandra/4913/primary/acisf04913N004_cntr_img2.fits") + +levels = np.geomspace(1.,99.,10) + +A = overplot_chandra(Stokes_UV, Stokes_Xr) +A.plot(levels=levels, SNRp_cut=3.0, SNRi_cut=20.0, zoom=1, savename='./plots/MRK463E/Chandra_overplot.pdf') + +B = overplot_chandra(Stokes_UV, Stokes_Xr, norm=LogNorm()) +B.plot(levels=levels, SNRp_cut=3.0, SNRi_cut=20.0, zoom=1, savename='./plots/MRK463E/Chandra_overplot_forced.pdf') + +C = overplot_pol(Stokes_UV, Stokes_IR) +C.plot(SNRp_cut=3.0, SNRi_cut=20.0, savename='./plots/MRK463E/IR_overplot.pdf') + +D = overplot_pol(Stokes_UV, Stokes_IR, norm=LogNorm()) +D.plot(SNRp_cut=3.0, SNRi_cut=30.0, vec_scale=2, norm=LogNorm(1e-18,1e-15), savename='./plots/MRK463E/IR_overplot_forced.pdf')