diff --git a/package/FOC_reduction.py b/package/FOC_reduction.py index 55bd806..a5c1b0f 100755 --- a/package/FOC_reduction.py +++ b/package/FOC_reduction.py @@ -41,12 +41,12 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= # Background estimation error_sub_type = "scott" # sqrt, sturges, rice, scott, freedman-diaconis (default) or shape (example (51, 51)) - subtract_error = 2.0 + subtract_error = 1.50 display_bkg = True # Data binning - pxsize = 4 - pxscale = "px" # pixel, arcsec or full + pxsize = 0.10 + pxscale = "arcsec" # pixel, arcsec or full rebin_operation = "sum" # sum or average # Alignement @@ -59,8 +59,8 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= # Smoothing smoothing_function = "combine" # gaussian_after, weighted_gaussian_after, gaussian, weighted_gaussian or combine - smoothing_FWHM = 1.5 # If None, no smoothing is done - smoothing_scale = "px" # pixel or arcsec + smoothing_FWHM = 0.150 # If None, no smoothing is done + smoothing_scale = "arcsec" # pixel or arcsec # Rotation rotate_North = True @@ -69,7 +69,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= P_cut = 3 # if >=1.0 cut on the signal-to-noise else cut on the confidence level in Q, U SNRi_cut = 1.0 # I measurments with SNR>30, which implies an uncertainty in P of 4.7%. flux_lim = None # lowest and highest flux displayed on plot, defaults to bkg and maximum in cut if None - scale_vec = 3 + scale_vec = 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 # Pipeline start diff --git a/package/lib/fits.py b/package/lib/fits.py index 0c82afd..62425f4 100755 --- a/package/lib/fits.py +++ b/package/lib/fits.py @@ -65,24 +65,25 @@ def get_obs_data(infiles, data_folder="", compute_flux=False): for wcs, header in zip(wcs_array, headers): new_wcs = wcs.deepcopy() if new_wcs.wcs.has_cd(): + # Update WCS with relevant information 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) - cdelt, orient = wcs_CD_to_PC(wcs.wcs.cd) - new_wcs.wcs.pc = wcs.wcs.cd.dot(np.diag(1.0 / cdelt)) - new_wcs.wcs.cdelt = cdelt + new_cdelt = np.linalg.eigvals(wcs.wcs.cd) + new_wcs.wcs.pc = wcs.wcs.cd.dot(np.diag(1.0 / new_cdelt)) + new_wcs.wcs.cdelt = new_cdelt try: header["ORIENTAT"] = float(header["ORIENTAT"]) except KeyError: - header["ORIENTAT"] = -orient - elif new_wcs.wcs.has_pc() and (new_wcs.wcs.cdelt[:2] != np.array([1.0, 1.0])).all(): + header["ORIENTAT"] = -wcs_CD_to_PC(new_wcs.wcs.cd)[1] + elif (new_wcs.wcs.cdelt[:2] != np.array([1.0, 1.0])).all(): try: header["ORIENTAT"] = float(header["ORIENTAT"]) except KeyError: header["ORIENTAT"] = -wcs_PA(new_wcs.wcs.pc, new_wcs.wcs.cdelt) else: - print("Could not compute ORIENTAT or CDELT from WCS") + print("Couldn't compute ORIENTAT from WCS") for key, val in new_wcs.to_header().items(): header[key] = val diff --git a/package/lib/plots.py b/package/lib/plots.py index c5a92d2..e428d28 100755 --- a/package/lib/plots.py +++ b/package/lib/plots.py @@ -363,7 +363,7 @@ def polarization_map( fig = plt.figure(figsize=(7 * ratiox, 7 * ratioy), layout="constrained") if ax is None: ax = fig.add_subplot(111, projection=wcs) - ax.set(aspect="equal", fc="k") + ax.set(aspect="equal", fc="k") # , ylim=[-0.05 * stkI.shape[0], 1.05 * stkI.shape[0]]) # fig.subplots_adjust(hspace=0, wspace=0, left=0.102, right=1.02) # ax.coords.grid(True, color='white', ls='dotted', alpha=0.5) @@ -444,7 +444,7 @@ def polarization_map( elif display.lower() in ["p", "pol", "pol_deg"]: # Display polarization degree map display = "p" - vmin, vmax = 0.0, 100.0 + vmin, vmax = 0.0, min(pol[np.isfinite(pol)].max(), 1.0) * 100.0 im = ax.imshow(pol * 100.0, vmin=vmin, vmax=vmax, aspect="equal", cmap=kwargs["cmap"], alpha=1.0) fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$P$ [%]") elif display.lower() in ["pa", "pang", "pol_ang"]: @@ -529,7 +529,7 @@ def polarization_map( PA_diluted = Stokes[0].header["PA_int"] PA_diluted_err = Stokes[0].header["sPA_int"] - plt.rcParams.update({"font.size": 11}) + plt.rcParams.update({"font.size": 12}) px_size = wcs.wcs.get_cdelt()[0] * 3600.0 px_sc = AnchoredSizeBar(ax.transData, 1.0 / px_size, "1 arcsec", 3, pad=0.25, sep=5, borderpad=0.25, frameon=False, size_vertical=0.005, color=font_color) north_dir = AnchoredDirectionArrows( @@ -2783,12 +2783,12 @@ class pol_map(object): def submit_save(expression): ax_text_save.set(visible=False) if expression != "": - save_fig, save_ax = plt.subplots(figsize=(12, 10), layout="constrained", subplot_kw=dict(projection=self.wcs)) + save_fig, save_ax = plt.subplots(figsize=(8, 6), layout="constrained", subplot_kw=dict(projection=self.wcs)) self.ax_cosmetics(ax=save_ax) self.display(fig=save_fig, ax=save_ax) self.pol_vector(fig=save_fig, ax=save_ax) self.pol_int(fig=save_fig, ax=save_ax) - save_fig.suptitle(r"{0:s} with $SNR_{{p}} \geq$ {1:d} and $SNR_{{I}} \geq$ {2:d}".format(self.targ, int(self.SNRp), int(self.SNRi))) + # save_fig.suptitle(r"{0:s} with $SNR_{{p}} \geq$ {1:d} and $SNR_{{I}} \geq$ {2:d}".format(self.targ, int(self.SNRp), int(self.SNRi))) if expression[-4:] not in [".png", ".jpg", ".pdf"]: expression += ".pdf" save_fig.savefig(expression, bbox_inches="tight", dpi=150, facecolor="None") @@ -3044,8 +3044,8 @@ class pol_map(object): sep_y=0.01, sep_x=0.01, back_length=0.0, - head_length=10.0, - head_width=10.0, + head_length=7.5, + head_width=7.5, angle=-self.Stokes[0].header["orientat"], color="white", text_props={"ec": None, "fc": "w", "alpha": 1, "lw": 0.4}, @@ -3425,7 +3425,7 @@ class pol_map(object): fontsize=12, xy=(0.01, 1.00), xycoords="axes fraction", - path_effects=[pe.withStroke(linewidth=0.5, foreground="k")], + path_effects=[pe.withStroke(linewidth=1.0, foreground="k")], verticalalignment="top", horizontalalignment="left", ) @@ -3461,7 +3461,7 @@ class pol_map(object): fontsize=12, xy=(0.01, 1.00), xycoords="axes fraction", - path_effects=[pe.withStroke(linewidth=0.5, foreground="k")], + path_effects=[pe.withStroke(linewidth=1.0, foreground="k")], verticalalignment="top", horizontalalignment="left", )