From d8365e984d9ae835be17e85934ea5a83dbb3889a Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Mon, 8 Jul 2024 14:48:18 +0200 Subject: [PATCH] fix plots aspect ratio, fits saving of reduction parameters --- package/lib/fits.py | 42 ++++++++++++++++++------------------- package/lib/plots.py | 49 +++++++++++++++++++++++++------------------- 2 files changed, 48 insertions(+), 43 deletions(-) diff --git a/package/lib/fits.py b/package/lib/fits.py index 9113f21..1506a29 100755 --- a/package/lib/fits.py +++ b/package/lib/fits.py @@ -16,6 +16,7 @@ from astropy.io import fits from astropy.wcs import WCS from .convex_hull import clean_ROI +from .utils import wcs_PA def get_obs_data(infiles, data_folder="", compute_flux=False): @@ -57,22 +58,20 @@ def get_obs_data(infiles, data_folder="", compute_flux=False): if new_wcs.wcs.has_cd() or (new_wcs.wcs.cdelt[:2] == np.array([1.0, 1.0])).all(): # Update WCS with relevant information if new_wcs.wcs.has_cd(): - old_cd = new_wcs.wcs.cd del new_wcs.wcs.cd keys = list(new_wcs.to_header().keys()) + ["CD1_1", "CD1_2", "CD1_3", "CD2_1", "CD2_2", "CD2_3", "CD3_1", "CD3_2", "CD3_3"] for key in keys: header.remove(key, ignore_missing=True) - new_cdelt = np.linalg.eig(old_cd)[0] - elif (new_wcs.wcs.cdelt == np.array([1.0, 1.0])).all() and (new_wcs.array_shape in [(512, 512), (1024, 512), (512, 1024), (1024, 1024)]): - old_cd = new_wcs.wcs.pc - new_wcs.wcs.pc = np.dot(old_cd, np.diag(1.0 / new_cdelt)) + new_cdelt = np.linalg.eigvals(wcs.wcs.cd) + new_cdelt.sort() + new_wcs.wcs.pc = wcs.wcs.cd.dot(np.diag(1.0 / new_cdelt)) new_wcs.wcs.cdelt = new_cdelt for key, val in new_wcs.to_header().items(): header[key] = val try: _ = header["ORIENTAT"] except KeyError: - header["ORIENTAT"] = -np.arccos(new_wcs.wcs.pc[0, 0]) * 180.0 / np.pi + header["ORIENTAT"] = wcs_PA(new_wcs.wcs.pc[1, 0], np.diag(new_wcs.wcs.pc).mean()) # force WCS for POL60 to have same pixel size as POL0 and POL120 is_pol60 = np.array([head["filtnam1"].lower() == "pol60" for head in headers], dtype=bool) @@ -130,7 +129,6 @@ def save_Stokes( Only returned if return_hdul is True. """ # Create new WCS object given the modified images - exp_tot = header_stokes['exptime'] new_wcs = WCS(header_stokes).deepcopy() if data_mask.shape != (1, 1): @@ -140,23 +138,23 @@ def save_Stokes( new_wcs.wcs.crpix = np.array(new_wcs.wcs.crpix) - vertex[0::-2] header = new_wcs.to_header() - header["TELESCOP"] = (header_stokes["telescop"] if "TELESCOP" in list(header_stokes.keys()) else "HST", "telescope used to acquire data") - header["INSTRUME"] = (header_stokes["instrume"] if "INSTRUME" in list(header_stokes.keys()) else "FOC", "identifier for instrument used to acuire data") - header["PHOTPLAM"] = (header_stokes["photplam"], "Pivot Wavelength") - header["PHOTFLAM"] = (header_stokes["photflam"], "Inverse Sensitivity in DN/sec/cm**2/Angst") - header["EXPTIME"] = (exp_tot, "Total exposure time in sec") - header["PROPOSID"] = (header_stokes["proposid"], "PEP proposal identifier for observation") - header["TARGNAME"] = (header_stokes["targname"], "Target name") - header["ORIENTAT"] = (np.arccos(new_wcs.wcs.pc[0, 0]) * 180.0 / np.pi, "Angle between North and the y-axis of the image") - header["FILENAME"] = (filename, "Original filename") + header["TELESCOP"] = (header_stokes["TELESCOP"] if "TELESCOP" in list(header_stokes.keys()) else "HST", "telescope used to acquire data") + header["INSTRUME"] = (header_stokes["INSTRUME"] if "INSTRUME" in list(header_stokes.keys()) else "FOC", "identifier for instrument used to acuire data") + header["PHOTPLAM"] = (header_stokes["PHOTPLAM"], "Pivot Wavelength") + header["PHOTFLAM"] = (header_stokes["PHOTFLAM"], "Inverse Sensitivity in DN/sec/cm**2/Angst") + header["EXPTIME"] = (header_stokes["EXPTIME"], "Total exposure time in sec") + header["PROPOSID"] = (header_stokes["PROPOSID"], "PEP proposal identifier for observation") + header["TARGNAME"] = (header_stokes["TARGNAME"], "Target name") + header["ORIENTAT"] = (header_stokes["ORIENTAT"], "Angle between North and the y-axis of the image") + header["FILENAME"] = (filename, "ORIGINAL FILENAME") header["BKG_TYPE"] = (header_stokes["BKG_TYPE"], "Bkg estimation method used during reduction") header["BKG_SUB"] = (header_stokes["BKG_SUB"], "Amount of bkg subtracted from images") - header["SMOOTH"] = (header_stokes["SMOOTH"], "Smoothing method used during reduction") - header["SAMPLING"] = (header_stokes["SAMPLING"], "Resampling performed during reduction") - header["P_INT"] = (header_stokes["P_int"], "Integrated polarization degree") - header["sP_INT"] = (header_stokes["sP_int"], "Integrated polarization degree error") - header["PA_INT"] = (header_stokes["PA_int"], "Integrated polarization angle") - header["sPA_INT"] = (header_stokes["sPA_int"], "Integrated polarization angle error") + header["SMOOTH"] = (header_stokes["SMOOTH"] if "SMOOTH" in list(header_stokes.keys()) else "None", "Smoothing method used during reduction") + header["SAMPLING"] = (header_stokes["SAMPLING"] if "SAMPLING" in list(header_stokes.keys()) else "None", "Resampling performed during reduction") + header["P_INT"] = (header_stokes["P_INT"], "Integrated polarization degree") + header["sP_INT"] = (header_stokes["sP_INT"], "Integrated polarization degree error") + header["PA_INT"] = (header_stokes["PA_INT"], "Integrated polarization angle") + header["sPA_INT"] = (header_stokes["sPA_INT"], "Integrated polarization angle error") # Crop Data to mask if data_mask.shape != (1, 1): diff --git a/package/lib/plots.py b/package/lib/plots.py index ecbebd3..a7b326e 100755 --- a/package/lib/plots.py +++ b/package/lib/plots.py @@ -182,9 +182,11 @@ def plot_Stokes(Stokes, savename=None, plots_folder=""): wcs = WCS(Stokes[0]).deepcopy() # Plot figure - plt.rcParams.update({"font.size": 10}) - fig, (axI, axQ, axU) = plt.subplots(ncols=3, figsize=(20, 6), subplot_kw=dict(projection=wcs)) - fig.subplots_adjust(hspace=0, wspace=0.75, bottom=0.01, top=0.99, left=0.08, right=0.95) + plt.rcParams.update({"font.size": 12}) + ratiox = max(int(stkI.shape[1]/stkI.shape[0]),1) + ratioy = max(int(stkI.shape[0]/stkI.shape[1]),1) + fig, (axI, axQ, axU) = plt.subplots(ncols=3, figsize=(20*ratiox, 8*ratioy), subplot_kw=dict(projection=wcs)) + fig.subplots_adjust(hspace=0, wspace=0.50, bottom=0.01, top=0.99, left=0.07, right=0.97) fig.suptitle("I, Q, U Stokes parameters") imI = axI.imshow(stkI, origin="lower", cmap="inferno") @@ -320,9 +322,11 @@ def polarization_map( print("No pixel with polarization information above requested SNR.") # Plot the map - plt.rcParams.update({"font.size": 10}) + plt.rcParams.update({"font.size": 12}) plt.rcdefaults() - fig, ax = plt.subplots(figsize=(10, 10), layout="constrained", subplot_kw=dict(projection=wcs)) + ratiox = max(int(stkI.shape[1]/stkI.shape[0]),1) + ratioy = max(int(stkI.shape[0]/stkI.shape[1]),1) + fig, ax = plt.subplots(figsize=(10*ratiox, 10*ratioy), layout="constrained", subplot_kw=dict(projection=wcs)) ax.set(aspect="equal", fc="k") # fig.subplots_adjust(hspace=0, wspace=0, left=0.102, right=1.02) @@ -439,17 +443,17 @@ def polarization_map( ax.transAxes, "E", "N", - length=-0.08, - fontsize=0.025, + length=-0.05, + fontsize=0.02, loc=1, - aspect_ratio=-1, + aspect_ratio=-(stkI.shape[1]/stkI.shape[0]), sep_y=0.01, sep_x=0.01, back_length=0.0, head_length=10.0, head_width=10.0, angle=-Stokes[0].header["orientat"], - text_props={"ec": "k", "fc": "w", "alpha": 1, "lw": -0.2}, + text_props={"ec": "k", "fc": "w", "alpha": 1, "lw": 0.4}, arrow_props={"ec": "k", "fc": "w", "alpha": 1, "lw": 1}, ) @@ -666,7 +670,7 @@ class align_maps(object): length=-0.08, fontsize=0.03, loc=1, - aspect_ratio=-1, + aspect_ratio=-(self.map_data.shape[1]/self.map_data.shape[0]), sep_y=0.01, sep_x=0.01, angle=-self.map_header["orientat"], @@ -724,7 +728,7 @@ class align_maps(object): length=-0.08, fontsize=0.03, loc=1, - aspect_ratio=-1, + aspect_ratio=-(self.other_data.shape[1]/self.other_data.shape[0]), sep_y=0.01, sep_x=0.01, angle=-self.other_header["orientat"], @@ -988,7 +992,7 @@ class overplot_radio(align_maps): length=-0.08, fontsize=0.03, loc=1, - aspect_ratio=-1, + aspect_ratio=-(stkI.shape[1]/stkI.shape[0]), sep_y=0.01, sep_x=0.01, angle=-self.Stokes_UV[0].header["orientat"], @@ -1190,7 +1194,7 @@ class overplot_chandra(align_maps): length=-0.08, fontsize=0.03, loc=1, - aspect_ratio=-1, + aspect_ratio=-(stkI.shape[1]/stkI.shape[0]), sep_y=0.01, sep_x=0.01, angle=-self.Stokes_UV[0].header["orientat"], @@ -1329,7 +1333,6 @@ class overplot_pol(align_maps): else: self.scale_vec = scale_vec step_vec = 1 - px_scale = np.abs(self.wcs_UV.wcs.get_cdelt()[0] / self.other_wcs.wcs.get_cdelt()[0]) self.X, self.Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0])) self.U, self.V = pol * np.cos(np.pi / 2.0 + pang * np.pi / 180.0), pol * np.sin(np.pi / 2.0 + pang * np.pi / 180.0) self.Q = self.ax_overplot.quiver( @@ -1339,7 +1342,7 @@ class overplot_pol(align_maps): self.V[::step_vec, ::step_vec], units="xy", angles="uv", - scale=px_scale / self.scale_vec, + scale=1. / self.scale_vec, scale_units="xy", pivot="mid", headwidth=0.0, @@ -1385,7 +1388,7 @@ class overplot_pol(align_maps): length=-0.08, fontsize=0.03, loc=1, - aspect_ratio=-1, + aspect_ratio=-(stkI.shape[1]/stkI.shape[0]), sep_y=0.01, sep_x=0.01, angle=-self.Stokes_UV[0].header["orientat"], @@ -1395,7 +1398,7 @@ class overplot_pol(align_maps): self.ax_overplot.add_artist(north_dir) pol_sc = AnchoredSizeBar( self.ax_overplot.transData, - self.scale_vec / px_scale, + self.scale_vec, r"$P$= 100%", 4, pad=0.5, @@ -1550,7 +1553,7 @@ class align_pol(object): length=-0.08, fontsize=0.025, loc=1, - aspect_ratio=-1, + aspect_ratio=-(stkI.shape[1]/stkI.shape[0]), sep_y=0.01, sep_x=0.01, back_length=0.0, @@ -1814,6 +1817,8 @@ class crop_map(object): # Write cropped map to new HDUList self.header_crop = deepcopy(header) self.header_crop.update(self.wcs_crop.to_header()) + if self.header_crop["FILENAME"][-4:] != "crop": + self.header_crop["FILENAME"] += "_crop" self.hdul_crop = fits.HDUList([fits.PrimaryHDU(self.data_crop, self.header_crop)]) self.rect_selector.clear() @@ -1936,6 +1941,8 @@ class crop_Stokes(crop_map): ) for dataset in self.hdul_crop: + if dataset.header["FILENAME"][-4:] != "crop": + dataset.header["FILENAME"] += "_crop" dataset.header["P_int"] = (P_diluted, "Integrated polarization degree") dataset.header["sP_int"] = (np.ceil(P_diluted_err * 1000.0) / 1000.0, "Integrated polarization degree error") dataset.header["PA_int"] = (PA_diluted, "Integrated polarization angle") @@ -2797,10 +2804,10 @@ class pol_map(object): ax.transAxes, "E", "N", - length=-0.08, - fontsize=0.025, + length=-0.05, + fontsize=0.02, loc=1, - aspect_ratio=-1, + aspect_ratio=-(self.I.shape[1]/self.I.shape[0]), sep_y=0.01, sep_x=0.01, back_length=0.0,