From 6b52b85002f9654818ecec42e6873dcbb7338356 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Fri, 7 Mar 2025 10:48:55 +0100 Subject: [PATCH 1/7] small improvments to displays resolve addition of display improvments --- package/lib/background.py | 13 +- package/lib/plots.py | 211 +++++++-------------------------- package/lib/utils.py | 2 +- package/src/emission_center.py | 6 +- 4 files changed, 51 insertions(+), 181 deletions(-) diff --git a/package/lib/background.py b/package/lib/background.py index ed9638c..d120ae6 100755 --- a/package/lib/background.py +++ b/package/lib/background.py @@ -86,20 +86,19 @@ def display_bkg(data, background, std_bkg, headers, histograms=None, binning=Non ) ax_h.plot([background[i] * convert_flux[i], background[i] * convert_flux[i]], [hist.min(), hist.max()], "x--", color="C{0:d}".format(i), alpha=0.8) if i == 0: - xmin, xmax = np.min(np.array(bins)[np.array(hist) > 1e2]) * convert_flux[0], np.max(np.array(bins)[np.array(hist) > 1e2]) * convert_flux[0] + xmin, xmax = np.min(np.array(bins)[np.array(hist) > 5e1]) * convert_flux[0], np.max(np.array(bins)[np.array(hist) > 5e1]) * convert_flux[0] else: xmin, xmax = ( - min(xmin, np.min(np.array(bins)[np.array(hist) > 1e2]) * convert_flux[0]), - max(xmax, np.max(np.array(bins)[np.array(hist) > 1e2]) * convert_flux[0]), + min(xmin, np.min(np.array(bins)[np.array(hist) > 5e1]) * convert_flux[0]), + max(xmax, np.max(np.array(bins)[np.array(hist) > 5e1]) * convert_flux[0]), ) if coeff is not None: # ax_h.plot(bins*convert_flux[i], gausspol(bins, *coeff[i]), '--', color="C{0:d}".format(i), alpha=0.8) ax_h.plot(bins * convert_flux[i], gauss(bins, *coeff[i]), "--", color="C{0:d}".format(i), alpha=0.8) ax_h.set_xscale("log") ax_h.set_yscale("log") - ax_h.set_ylim([100.0, np.max([hist.max() for hist in histograms])]) - # ax_h.set_xlim([np.min(background * convert_flux) * 1e-2, np.max(background * convert_flux) * 1e2]) - ax_h.set_xlim([xmin, xmax]) + ax_h.set_ylim([5e1, np.max([hist.max() for hist in histograms])]) + ax_h.set_xlim([max(xmin, np.min(background * convert_flux) * 1e-2), min(xmax, np.max(background * convert_flux) * 1e2)]) ax_h.set_xlabel(r"Flux [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") ax_h.set_ylabel(r"Number of pixels in bin") ax_h.set_title("Histogram for each observation") @@ -135,7 +134,7 @@ def display_bkg(data, background, std_bkg, headers, histograms=None, binning=Non ax2.set(xlabel="pixel offset", ylabel="pixel offset", aspect="equal") fig2.subplots_adjust(hspace=0, wspace=0, right=1.0) - fig2.colorbar(im2, ax=ax2, location="right", aspect=50, pad=0.025, label=r"Flux [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") + fig2.colorbar(im2, ax=ax2, location="right", shrink=0.60, aspect=50, pad=0.025, label=r"Flux [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") if savename is not None: this_savename = deepcopy(savename) diff --git a/package/lib/plots.py b/package/lib/plots.py index 604d58f..a2278a9 100755 --- a/package/lib/plots.py +++ b/package/lib/plots.py @@ -428,7 +428,7 @@ def polarization_map( P_cut = 0.99 else: # MaskP on the confidence value - SNRp_cut = 5.0 + SNRp_cut = 3.0 maskP = confP > P_cut mask = np.logical_and(maskI, maskP) @@ -504,72 +504,32 @@ def polarization_map( if display.lower() in ["i", "intensity"]: # If no display selected, show intensity map display = "i" - if flux_lim is None: - if mask.sum() > 0.0: - vmin, vmax = ( - 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][mask]) * convert_flux), - np.max(stkI[stkI > 0.0] * convert_flux), - ) - else: - vmin, vmax = ( - 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][stkI > 0.0]) * convert_flux), - np.max(stkI[stkI > 0.0] * convert_flux), - ) - else: + if flux_lim is not None: vmin, vmax = flux_lim - im = ax.imshow( - stkI * convert_flux, - norm=LogNorm(vmin, vmax), - aspect="equal", - cmap=kwargs["cmap"], - alpha=1.0, - ) - fig.colorbar( - im, - ax=ax, - aspect=30, - shrink=0.60, - pad=0.025, - label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]", - ) + elif mask.sum() > 0.0: + vmin, vmax = 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][mask]) * convert_flux), np.max(stkI[stkI > 0.0] * convert_flux) + else: + vmin, vmax = 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][stkI > 0.0]) * convert_flux), np.max(stkI[stkI > 0.0] * convert_flux) + im = ax.imshow(stkI * convert_flux, norm=LogNorm(vmin, vmax), aspect="equal", cmap=kwargs["cmap"], alpha=1.0) + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") levelsI = np.array([0.8, 2.0, 5.0, 10.0, 20.0, 50.0]) / 100.0 * vmax print("Flux density contour levels : ", levelsI) ax.contour(stkI * convert_flux, levels=levelsI, colors="grey", linewidths=0.5) elif display.lower() in ["pf", "pol_flux"]: # Display polarization flux display = "pf" - if flux_lim is None: - if mask.sum() > 0.0: - vmin, vmax = ( - 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][mask]) * convert_flux), - np.max(stkI[stkI > 0.0] * convert_flux), - ) - pfmax = (stkI[mask] * pol[mask] * convert_flux).max() - else: - vmin, vmax = ( - 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][stkI > 0.0]) * convert_flux), - np.max(stkI[stkI > 0.0] * convert_flux), - ) - pfmax = (stkI[stkI > 0.0] * pol[stkI > 0.0] * convert_flux).max() - else: + if flux_lim is not None: vmin, vmax = flux_lim pfmax = (stkI[stkI > 0.0] * pol[stkI > 0.0] * convert_flux).max() - im = ax.imshow( - stkI * convert_flux * pol, - norm=LogNorm(vmin, vmax), - aspect="equal", - cmap=kwargs["cmap"], - alpha=1.0, - ) - fig.colorbar( - im, - ax=ax, - aspect=30, - shrink=0.60, - pad=0.025, - label=r"$F_{\lambda} \cdot P$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]", - ) - # levelsPf = np.linspace(0.0175, 0.50, 5) * pfmax + elif mask.sum() > 0.0: + vmin, vmax = 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][mask]) * convert_flux), np.max(stkI[stkI > 0.0] * convert_flux) + pfmax = (stkI[mask] * pol[mask] * convert_flux).max() + else: + vmin, vmax = 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][stkI > 0.0]) * convert_flux), np.max(stkI[stkI > 0.0] * convert_flux) + pfmax = (stkI[stkI > 0.0] * pol[stkI > 0.0] * convert_flux).max() + im = ax.imshow(stkI * convert_flux * pol, norm=LogNorm(vmin, vmax), aspect="equal", cmap=kwargs["cmap"], alpha=1.0) + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$F_{\lambda} \cdot P$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") + # levelsPf = np.linspace(0.0.60, 0.50, 5) * pfmax levelsPf = np.array([1.73, 13.0, 33.0, 66.0]) / 100.0 * pfmax print("Polarized flux density contour levels : ", levelsPf) ax.contour(stkI * convert_flux * pol, levels=levelsPf, colors="grey", linewidths=0.5) @@ -577,52 +537,24 @@ def polarization_map( # Display polarization degree map display = "p" vmin, vmax = 0.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=30, shrink=0.60, pad=0.025, label=r"$P$ [%]") + 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"]: # Display polarization degree map display = "pa" vmin, vmax = 0.0, 180.0 - im = ax.imshow( - princ_angle(pang), - vmin=vmin, - vmax=vmax, - aspect="equal", - cmap=kwargs["cmap"], - alpha=1.0, - ) - fig.colorbar(im, ax=ax, aspect=30, shrink=0.60, pad=0.025, label=r"$\theta_P$ [°]") + im = ax.imshow(princ_angle(pang), 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"$\theta_P$ [°]") elif display.lower() in ["s_p", "pol_err", "pol_deg_err"]: # Display polarization degree error map display = "s_p" if (SNRp > P_cut).any(): vmin, vmax = 0.0, np.max([pol_err[SNRp > P_cut].max(), 1.0]) * 100.0 - im = ax.imshow( - pol_err * 100.0, - vmin=vmin, - vmax=vmax, - aspect="equal", - cmap="inferno_r", - alpha=1.0, - ) + im = ax.imshow(pol_err * 100.0, vmin=vmin, vmax=vmax, aspect="equal", cmap="inferno_r", alpha=1.0) else: vmin, vmax = 0.0, 100.0 - im = ax.imshow( - pol_err * 100.0, - vmin=vmin, - vmax=vmax, - aspect="equal", - cmap="inferno_r", - alpha=1.0, - ) - fig.colorbar(im, ax=ax, aspect=30, shrink=0.60, pad=0.025, label=r"$\sigma_P$ [%]") + im = ax.imshow(pol_err * 100.0, vmin=vmin, vmax=vmax, aspect="equal", cmap="inferno_r", alpha=1.0) + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$\sigma_P$ [%]") elif display.lower() in ["s_i", "i_err"]: # Display intensity error map display = "s_i" @@ -631,73 +563,34 @@ def polarization_map( 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][stk_cov[0, 0] > 0.0]) * convert_flux), np.max(np.sqrt(stk_cov[0, 0][stk_cov[0, 0] > 0.0]) * convert_flux), ) - im = ax.imshow( - np.sqrt(stk_cov[0, 0]) * convert_flux, - norm=LogNorm(vmin, vmax), - aspect="equal", - cmap="inferno_r", - alpha=1.0, - ) + im = ax.imshow(np.sqrt(stk_cov[0, 0]) * convert_flux, norm=LogNorm(vmin, vmax), aspect="equal", cmap="inferno_r", alpha=1.0) else: - im = ax.imshow( - np.sqrt(stk_cov[0, 0]) * convert_flux, - aspect="equal", - cmap=kwargs["cmap"], - alpha=1.0, - ) - fig.colorbar( - im, - ax=ax, - aspect=30, - shrink=0.60, - pad=0.025, - label=r"$\sigma_I$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]", - ) + im = ax.imshow(np.sqrt(stk_cov[0, 0]) * convert_flux, aspect="equal", cmap=kwargs["cmap"], alpha=1.0) + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$\sigma_I$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") elif display.lower() in ["snri"]: # Display I_stokes signal-to-noise map display = "snri" vmin, vmax = 0.0, np.max(SNRi[np.isfinite(SNRi)]) if vmax * 0.99 > SNRi_cut: - im = ax.imshow( - SNRi, - vmin=vmin, - vmax=vmax, - aspect="equal", - cmap=kwargs["cmap"], - alpha=1.0, - ) + im = ax.imshow(SNRi, vmin=vmin, vmax=vmax, aspect="equal", cmap=kwargs["cmap"], alpha=1.0) levelsSNRi = np.linspace(SNRi_cut, vmax * 0.99, 5).astype(int) print("SNRi contour levels : ", levelsSNRi) ax.contour(SNRi, levels=levelsSNRi, colors="grey", linewidths=0.5) else: im = ax.imshow(SNRi, aspect="equal", cmap=kwargs["cmap"], alpha=1.0) - fig.colorbar( - im, - ax=ax, - aspect=30, - shrink=0.60, - pad=0.025, - label=r"$I_{Stokes}/\sigma_{I}$", - ) + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$I_{Stokes}/\sigma_{I}$") elif display.lower() in ["snr", "snrp"]: # Display polarization degree signal-to-noise map display = "snrp" vmin, vmax = 0.0, np.max(SNRp[np.isfinite(SNRp)]) if vmax * 0.99 > SNRp_cut: - im = ax.imshow( - SNRp, - vmin=vmin, - vmax=vmax, - aspect="equal", - cmap=kwargs["cmap"], - alpha=1.0, - ) + im = ax.imshow(SNRp, vmin=vmin, vmax=vmax, aspect="equal", cmap=kwargs["cmap"], alpha=1.0) levelsSNRp = np.linspace(P_cut, vmax * 0.99, 5).astype(int) print("SNRp contour levels : ", levelsSNRp) ax.contour(SNRp, levels=levelsSNRp, colors="grey", linewidths=0.5) else: im = ax.imshow(SNRp, aspect="equal", cmap=kwargs["cmap"], alpha=1.0) - fig.colorbar(im, ax=ax, aspect=30, shrink=0.60, pad=0.025, label=r"$P/\sigma_{P}$") + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$P/\sigma_{P}$") elif display.lower() in ["conf", "confp"]: # Display polarization degree signal-to-noise map display = "confp" @@ -706,34 +599,17 @@ def polarization_map( levelsconfp = np.array([0.500, 0.900, 0.990, 0.999]) print("confp contour levels : ", levelsconfp) ax.contour(confP, levels=levelsconfp, colors="grey", linewidths=0.5) - fig.colorbar(im, ax=ax, aspect=30, shrink=0.60, pad=0.025, label=r"$Conf_{P}$") + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$Conf_{P}$") else: # Defaults to intensity map - if mask.sum() > 0.0: - vmin, vmax = ( - 1.0 * np.mean(np.sqrt(stk_cov[0, 0][mask]) * convert_flux), - np.max(stkI[stkI > 0.0] * convert_flux), - ) + if flux_lim is not None: + vmin, vmax = flux_lim + elif mask.sum() > 0.0: + vmin, vmax = 1.0 * np.mean(np.sqrt(stk_cov[0, 0][mask]) * convert_flux), np.max(stkI[stkI > 0.0] * convert_flux) else: - vmin, vmax = ( - 1.0 * np.mean(np.sqrt(stk_cov[0, 0][stkI > 0.0]) * convert_flux), - np.max(stkI[stkI > 0.0] * convert_flux), - ) - im = ax.imshow( - stkI * convert_flux, - norm=LogNorm(vmin, vmax), - aspect="equal", - cmap=kwargs["cmap"], - alpha=1.0, - ) - fig.colorbar( - im, - ax=ax, - aspect=30, - shrink=0.60, - pad=0.025, - label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]", - ) + vmin, vmax = 1.0 * np.mean(np.sqrt(stk_cov[0, 0][stkI > 0.0]) * convert_flux), np.max(stkI[stkI > 0.0] * convert_flux) + im = ax.imshow(stkI * convert_flux, norm=LogNorm(vmin, vmax), aspect="equal", cmap=kwargs["cmap"], alpha=1.0) + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]") # Get integrated flux values from sum I_diluted = stkI[data_mask].sum() @@ -3148,12 +3024,7 @@ class pol_map(object): self.snr_conf = 0 b_snr_conf.label.set_text("Conf") self.s_P_cut = Slider( - self.ax_P_cut, - r"$SNR^{P}_{cut}$", - 1.0, - int(SNRp_max * 0.95), - valstep=1, - valinit=self.P_cut if P_cut >= 1.0 else 3, + self.ax_P_cut, r"$SNR^{P}_{cut}$", 1.0, max(int(SNRp_max * 0.95), 3), valstep=1, valinit=self.P_cut if P_cut >= 1.0 else 3 ) else: self.snr_conf = 1 diff --git a/package/lib/utils.py b/package/lib/utils.py index 214f0f2..86dcf8c 100755 --- a/package/lib/utils.py +++ b/package/lib/utils.py @@ -109,7 +109,7 @@ def CenterConf(mask, PA, sPA): conf[np.isfinite(PA)] = gammaincc(0.5, 0.5 * chi2[np.isfinite(PA)]) c0 = np.unravel_index(np.argmax(conf), conf.shape)[::-1] - result = minimize(chisq, c0, bounds=[(0, PA.shape[1]), (0.0, PA.shape[0])], method="TNC") + result = minimize(chisq, c0, bounds=[(0, PA.shape[1]), (0.0, PA.shape[0])], method="trust-constr") if result.success: print("Center of emission found") else: diff --git a/package/src/emission_center.py b/package/src/emission_center.py index c793206..c9fe314 100755 --- a/package/src/emission_center.py +++ b/package/src/emission_center.py @@ -44,7 +44,7 @@ def main(infile, P_cut=0.99, target=None, display="pf", output_dir=None): target = Stokes[0].header["TARGNAME"] fig = figure(figsize=(8, 8.5), layout="constrained") - fig, ax = polarization_map(Stokes, P_cut=P_cut, step_vec=1, scale_vec=3, display=display, fig=fig, width=0.33, linewidth=0.5) + fig, ax = polarization_map(Stokes, P_cut=P_cut, step_vec=1, scale_vec=5, display=display, fig=fig, width=0.33, linewidth=0.5) ax.plot(*Stokescenter, marker="+", color="k", lw=3) ax.plot(*Stokescenter, marker="+", color="red", lw=1.5, label="Best confidence for center: {0}".format(Stokespos.to_string("hmsdms"))) @@ -62,11 +62,11 @@ def main(infile, P_cut=0.99, target=None, display="pf", output_dir=None): # handles.append(Rectangle((0, 0), 1, 1, fill=False, ls="--", ec=snr3cont.get_edgecolor()[0])) # labels.append(r"$SNR_P \geq$ 4 contour") # handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=snr4cont.get_edgecolor()[0])) - ax.legend(handles=handles, labels=labels, bbox_to_anchor=(0.0, -0.02, 1.0, 0.01), loc="upper left", mode="expand", borderaxespad=0.0) + ax.legend(handles=handles, labels=labels, bbox_to_anchor=(-0.05, -0.02, 1.10, 0.01), loc="upper left", mode="expand", borderaxespad=0.0) if output_dir is not None: filename = pathjoin(output_dir, "%s_center.pdf" % target) - fig.savefig(filename, dpi=300, facecolor="None") + fig.savefig(filename, bbox_inches='tight', dpi=300, facecolor="None") output.append(filename) show() return output From efca472af168a4390a4a4e8ce73f8fdc2dfbcb69 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Tue, 11 Mar 2025 16:13:41 +0100 Subject: [PATCH 2/7] small improvments to overplot_pol --- package/lib/plots.py | 105 ++++++++++++++++++++++--------------------- 1 file changed, 53 insertions(+), 52 deletions(-) diff --git a/package/lib/plots.py b/package/lib/plots.py index a2278a9..40ec567 100755 --- a/package/lib/plots.py +++ b/package/lib/plots.py @@ -1599,7 +1599,7 @@ class overplot_pol(align_maps): def overplot( self, - levels=None, + levels="Default", P_cut=0.99, SNRi_cut=1.0, step_vec=1, @@ -1668,9 +1668,7 @@ class overplot_pol(align_maps): ["linewidth", [["linewidth", 0.3 * self.px_scale]]], ]: try: - test = kwargs[key] - if isinstance(test, LogNorm): - kwargs[key] = LogNorm(vmin, vmax) + _ = kwargs[key] except KeyError: for key_i, val_i in value: kwargs[key_i] = val_i @@ -1703,7 +1701,7 @@ class overplot_pol(align_maps): else: self.ax_overplot.set_facecolor("white") font_color = "black" - if hasattr(kwargs, "norm"): + if "norm" in kwargs.keys(): self.im = self.ax_overplot.imshow( other_data * self.other_convert, alpha=1.0, @@ -1766,50 +1764,51 @@ class overplot_pol(align_maps): # Display Stokes as contours disptypestr = "" - if disptype.lower() == "p": - disptypestr = "polarization degree" - if levels is None: - levels = np.array([2.0, 5.0, 10.0, 20.0, 90.0]) / 100.0 * np.max(pol[stkI > 0.0]) - cont_stk = self.ax_overplot.contour( - pol * 100.0, - levels=levels * 100.0, - colors="grey", - alpha=0.75, - transform=self.ax_overplot.get_transform(self.wcs_UV), - ) - if disptype.lower() == "pf": - disptypestr = "polarized flux" - if levels is None: - levels = np.array([2.0, 5.0, 10.0, 20.0, 90.0]) / 100.0 * np.max(stkI[stkI > 0.0] * pol[stkI > 0.0]) * self.map_convert - cont_stk = self.ax_overplot.contour( - stkI * pol * self.map_convert, - levels=levels, - colors="grey", - alpha=0.75, - transform=self.ax_overplot.get_transform(self.wcs_UV), - ) - if disptype.lower() == "snri": - disptypestr = "Stokes I signal-to-noise" - if levels is None: - levels = np.array([2.0, 5.0, 10.0, 20.0, 90.0]) / 100.0 * np.max(SNRi[stk_cov[0, 0] > 0.0]) - cont_stk = self.ax_overplot.contour( - SNRi, - levels=levels, - colors="grey", - alpha=0.75, - transform=self.ax_overplot.get_transform(self.wcs_UV), - ) - else: # default to intensity contours - disptypestr = "Stokes I" - if levels is None: - levels = np.array([2.0, 5.0, 10.0, 20.0, 90.0]) / 100.0 * np.max(stkI[stkI > 0.0]) * self.map_convert - cont_stk = self.ax_overplot.contour( - stkI * self.map_convert, - levels=levels, - colors="grey", - alpha=0.75, - transform=self.ax_overplot.get_transform(self.wcs_UV), - ) + if levels is not None: + if disptype.lower() == "p": + disptypestr = "polarization degree" + if levels == "Default": + levels = np.array([2.0, 5.0, 10.0, 20.0, 90.0]) / 100.0 * np.max(pol[stkI > 0.0]) + cont_stk = self.ax_overplot.contour( + pol * 100.0, + levels=levels * 100.0, + colors="grey", + alpha=0.75, + transform=self.ax_overplot.get_transform(self.wcs_UV), + ) + elif disptype.lower() == "pf": + disptypestr = "polarized flux" + if levels == "Default": + levels = np.array([2.0, 5.0, 10.0, 20.0, 90.0]) / 100.0 * np.max(stkI[stkI > 0.0] * pol[stkI > 0.0]) * self.map_convert + cont_stk = self.ax_overplot.contour( + stkI * pol * self.map_convert, + levels=levels, + colors="grey", + alpha=0.75, + transform=self.ax_overplot.get_transform(self.wcs_UV), + ) + elif disptype.lower() == "snri": + disptypestr = "Stokes I signal-to-noise" + if levels == "Default": + levels = np.array([2.0, 5.0, 10.0, 20.0, 90.0]) / 100.0 * np.max(SNRi[stk_cov[0, 0] > 0.0]) + cont_stk = self.ax_overplot.contour( + SNRi, + levels=levels, + colors="grey", + alpha=0.75, + transform=self.ax_overplot.get_transform(self.wcs_UV), + ) + else: # default to intensity contours + disptypestr = "Stokes I" + if levels == "Default": + levels = np.array([2.0, 5.0, 10.0, 20.0, 90.0]) / 100.0 * np.max(stkI[stkI > 0.0]) * self.map_convert + cont_stk = self.ax_overplot.contour( + stkI * self.map_convert, + levels=levels, + colors="grey", + alpha=0.75, + transform=self.ax_overplot.get_transform(self.wcs_UV), + ) # self.ax_overplot.clabel(cont_stk, inline=False, colors="k", fontsize=7) # Display pixel scale and North direction @@ -1881,8 +1880,10 @@ class overplot_pol(align_maps): handles[np.argmax([li == "{0:s} polarization map".format(self.map_observer) for li in labels])] = FancyArrowPatch( (0, 0), (0, 1), arrowstyle="-", fc="w", ec="k", lw=2 ) - labels.append("{0:s} {1:s} contour".format(self.map_observer, disptypestr)) - handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=cont_stk.get_edgecolor()[0])) + if disptypestr != "": + labels.append("{0:s} {1:s} contour".format(self.map_observer, disptypestr)) + handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=cont_stk.get_edgecolor()[0])) + disptypestr += " contours" self.legend = self.ax_overplot.legend( handles=handles, labels=labels, @@ -1893,7 +1894,7 @@ class overplot_pol(align_maps): ) self.fig_overplot.suptitle( - "{0:s} observation from {1:s} overplotted with {2:s} contours from {3:s}".format(obj, self.other_observer, vecstr + disptypestr, self.map_observer), + "{0:s} observation from {1:s} overplotted with {2:s} from {3:s}".format(obj, self.other_observer, vecstr + disptypestr, self.map_observer), wrap=True, ) From 1a765af73bd7d0d37a450f81f87b8608f478e5e6 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Tue, 11 Mar 2025 16:14:14 +0100 Subject: [PATCH 3/7] small improvments to ConfCenter et emission_center --- package/lib/utils.py | 2 +- package/src/emission_center.py | 95 ++++++++++++++++++++++++++++++---- 2 files changed, 85 insertions(+), 12 deletions(-) diff --git a/package/lib/utils.py b/package/lib/utils.py index 86dcf8c..214f0f2 100755 --- a/package/lib/utils.py +++ b/package/lib/utils.py @@ -109,7 +109,7 @@ def CenterConf(mask, PA, sPA): conf[np.isfinite(PA)] = gammaincc(0.5, 0.5 * chi2[np.isfinite(PA)]) c0 = np.unravel_index(np.argmax(conf), conf.shape)[::-1] - result = minimize(chisq, c0, bounds=[(0, PA.shape[1]), (0.0, PA.shape[0])], method="trust-constr") + result = minimize(chisq, c0, bounds=[(0, PA.shape[1]), (0.0, PA.shape[0])], method="TNC") if result.success: print("Center of emission found") else: diff --git a/package/src/emission_center.py b/package/src/emission_center.py index c9fe314..097362c 100755 --- a/package/src/emission_center.py +++ b/package/src/emission_center.py @@ -23,7 +23,12 @@ def main(infile, P_cut=0.99, target=None, display="pf", output_dir=None): stkI = Stokes["I_STOKES"].data QN, UN, QN_ERR, UN_ERR = np.full((4, stkI.shape[0], stkI.shape[1]), np.nan) for sflux, nflux in zip( - [Stokes["Q_STOKES"].data, Stokes["U_STOKES"].data, np.sqrt(Stokes["IQU_COV_MATRIX"].data[1, 1]), np.sqrt(Stokes["IQU_COV_MATRIX"].data[2, 2])], + [ + Stokes["Q_STOKES"].data, + Stokes["U_STOKES"].data, + np.sqrt(Stokes["IQU_COV_MATRIX"].data[1, 1]), + np.sqrt(Stokes["IQU_COV_MATRIX"].data[2, 2]), + ], [QN, UN, QN_ERR, UN_ERR], ): nflux[stkI > 0.0] = sflux[stkI > 0.0] / stkI[stkI > 0.0] @@ -44,10 +49,25 @@ def main(infile, P_cut=0.99, target=None, display="pf", output_dir=None): target = Stokes[0].header["TARGNAME"] fig = figure(figsize=(8, 8.5), layout="constrained") - fig, ax = polarization_map(Stokes, P_cut=P_cut, step_vec=1, scale_vec=5, display=display, fig=fig, width=0.33, linewidth=0.5) + fig, ax = polarization_map( + Stokes, + P_cut=P_cut, + step_vec=1, + scale_vec=5, + display=display, + fig=fig, + width=0.33, + linewidth=0.5, + ) ax.plot(*Stokescenter, marker="+", color="k", lw=3) - ax.plot(*Stokescenter, marker="+", color="red", lw=1.5, label="Best confidence for center: {0}".format(Stokespos.to_string("hmsdms"))) + ax.plot( + *Stokescenter, + marker="+", + color="red", + lw=1.5, + label="Best confidence for center: {0}".format(Stokespos.to_string("hmsdms")), + ) ax.contour(Stokescentconf, [0.01], colors="k", linewidths=3) confcentcont = ax.contour(Stokescentconf, [0.01], colors="red") # confcont = ax.contour(Stokesconf, [0.9905], colors="r") @@ -62,11 +82,18 @@ def main(infile, P_cut=0.99, target=None, display="pf", output_dir=None): # handles.append(Rectangle((0, 0), 1, 1, fill=False, ls="--", ec=snr3cont.get_edgecolor()[0])) # labels.append(r"$SNR_P \geq$ 4 contour") # handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=snr4cont.get_edgecolor()[0])) - ax.legend(handles=handles, labels=labels, bbox_to_anchor=(-0.05, -0.02, 1.10, 0.01), loc="upper left", mode="expand", borderaxespad=0.0) + ax.legend( + handles=handles, + labels=labels, + bbox_to_anchor=(-0.05, -0.02, 1.10, 0.01), + loc="upper left", + mode="expand", + borderaxespad=0.0, + ) if output_dir is not None: filename = pathjoin(output_dir, "%s_center.pdf" % target) - fig.savefig(filename, bbox_inches='tight', dpi=300, facecolor="None") + fig.savefig(filename, bbox_inches="tight", dpi=300, facecolor="None") output.append(filename) show() return output @@ -76,11 +103,57 @@ if __name__ == "__main__": import argparse parser = argparse.ArgumentParser(description="Look for the center of emission for a given reduced observation") - parser.add_argument("-t", "--target", metavar="targetname", required=False, help="the name of the target", type=str, default=None) - parser.add_argument("-f", "--file", metavar="path", required=False, help="The full or relative path to the data product", type=str, default=None) - parser.add_argument("-c", "--pcut", metavar="pcut", required=False, help="The polarization cut for the data mask", type=float, default=0.99) - parser.add_argument("-d", "--display", metavar="display", required=False, help="The map on which to display info", type=str, default="pf") - parser.add_argument("-o", "--output_dir", metavar="directory_path", required=False, help="output directory path for the plots", type=str, default="./data") + parser.add_argument( + "-t", + "--target", + metavar="targetname", + required=False, + help="the name of the target", + type=str, + default=None, + ) + parser.add_argument( + "-f", + "--file", + metavar="path", + required=False, + help="The full or relative path to the data product", + type=str, + default=None, + ) + parser.add_argument( + "-c", + "--pcut", + metavar="pcut", + required=False, + help="The polarization cut for the data mask", + type=float, + default=0.99, + ) + parser.add_argument( + "-d", + "--display", + metavar="display", + required=False, + help="The map on which to display info", + type=str, + default="pf", + ) + parser.add_argument( + "-o", + "--output_dir", + metavar="directory_path", + required=False, + help="output directory path for the plots", + type=str, + default="./data", + ) args = parser.parse_args() - exitcode = main(infile=args.file, P_cut=args.pcut, target=args.target, display=args.display, output_dir=args.output_dir) + exitcode = main( + infile=args.file, + P_cut=args.pcut, + target=args.target, + display=args.display, + output_dir=args.output_dir, + ) print("Written to: ", exitcode) From 8da199880b28ca1e696b87245acdb8cdc18af5ba Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Fri, 14 Mar 2025 14:20:38 +0100 Subject: [PATCH 4/7] rollback CenterConf to simple CDF of chi2 for 2 dof --- package/lib/utils.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/package/lib/utils.py b/package/lib/utils.py index 214f0f2..ee9fc5b 100755 --- a/package/lib/utils.py +++ b/package/lib/utils.py @@ -91,6 +91,7 @@ def CenterConf(mask, PA, sPA): chi2 = np.full(PA.shape, np.nan, dtype=np.float64) conf = np.full(PA.shape, -1.0, dtype=np.float64) yy, xx = np.indices(PA.shape) + Nobs = np.sum(mask) def ideal(c): itheta = np.full(PA.shape, np.nan) @@ -99,19 +100,18 @@ def CenterConf(mask, PA, sPA): return princ_angle(itheta) def chisq(c): - return np.sum((princ_angle(PA[mask]) - ideal((c[0], c[1]))[mask]) ** 2 / sPA[mask] ** 2) / np.sum(mask) + return np.sum((princ_angle(PA[mask]) - ideal((c[0], c[1]))[mask]) ** 2 / sPA[mask] ** 2) / (Nobs - 2) for x, y in zip(xx[np.isfinite(PA)], yy[np.isfinite(PA)]): chi2[y, x] = chisq((x, y)) from scipy.optimize import minimize - from scipy.special import gammaincc - conf[np.isfinite(PA)] = gammaincc(0.5, 0.5 * chi2[np.isfinite(PA)]) + conf[np.isfinite(PA)] = 0.5 * np.exp(-0.5 * chi2[np.isfinite(PA)]) c0 = np.unravel_index(np.argmax(conf), conf.shape)[::-1] - result = minimize(chisq, c0, bounds=[(0, PA.shape[1]), (0.0, PA.shape[0])], method="TNC") + result = minimize(chisq, c0, bounds=[(0, PA.shape[1]), (0.0, PA.shape[0])], method="trust-constr") if result.success: - print("Center of emission found") + print("Center of emission found: reduced chi_squared {0:.2f}/{1:d}={2:.2f}".format(chisq(result.x) * (Nobs - 2), Nobs - 2, chisq(result.x))) else: print("Center of emission not found", result) return conf, result.x From 749a08eae0fb3c3167e710d01614541e64dbd560 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Fri, 14 Mar 2025 14:22:18 +0100 Subject: [PATCH 5/7] default background estimation to Scott statistics --- package/lib/background.py | 22 +++++++++++----------- package/lib/reduction.py | 10 ++++------ 2 files changed, 15 insertions(+), 17 deletions(-) diff --git a/package/lib/background.py b/package/lib/background.py index d120ae6..99a09f2 100755 --- a/package/lib/background.py +++ b/package/lib/background.py @@ -294,7 +294,7 @@ def bkg_hist(data, error, mask, headers, sub_type=None, subtract_error=True, dis sub_type : str or int, optional If str, statistic rule to be used for the number of bins in counts/s. If int, number of bins for the counts/s histogram. - Defaults to "Freedman-Diaconis". + Defaults to "Scott". subtract_error : float or bool, optional If float, factor to which the estimated background should be multiplied If False the background is not subtracted. @@ -336,25 +336,25 @@ def bkg_hist(data, error, mask, headers, sub_type=None, subtract_error=True, dis if sub_type is not None: if isinstance(sub_type, int): n_bins = sub_type - elif sub_type.lower() in ["sqrt"]: + elif sub_type.lower() in ["square-root", "squareroot", "sqrt"]: n_bins = np.fix(np.sqrt(image[n_mask].size)).astype(int) # Square-root elif sub_type.lower() in ["sturges"]: n_bins = np.ceil(np.log2(image[n_mask].size)).astype(int) + 1 # Sturges elif sub_type.lower() in ["rice"]: n_bins = 2 * np.fix(np.power(image[n_mask].size, 1 / 3)).astype(int) # Rice - elif sub_type.lower() in ["scott"]: - n_bins = np.fix((image[n_mask].max() - image[n_mask].min()) / (3.5 * image[n_mask].std() / np.power(image[n_mask].size, 1 / 3))).astype( - int - ) # Scott - else: + elif sub_type.lower() in ["freedman-diaconis", "freedmandiaconis", "freedman", "diaconis"]: n_bins = np.fix( (image[n_mask].max() - image[n_mask].min()) / (2 * np.subtract(*np.percentile(image[n_mask], [75, 25])) / np.power(image[n_mask].size, 1 / 3)) ).astype(int) # Freedman-Diaconis - else: - n_bins = np.fix( - (image[n_mask].max() - image[n_mask].min()) / (2 * np.subtract(*np.percentile(image[n_mask], [75, 25])) / np.power(image[n_mask].size, 1 / 3)) - ).astype(int) # Freedman-Diaconis + else: # Fallback + n_bins = np.fix((image[n_mask].max() - image[n_mask].min()) / (3.5 * image[n_mask].std() / np.power(image[n_mask].size, 1 / 3))).astype( + int + ) # Scott + else: # Default statistic + n_bins = np.fix((image[n_mask].max() - image[n_mask].min()) / (3.5 * image[n_mask].std() / np.power(image[n_mask].size, 1 / 3))).astype( + int + ) # Scott hist, bin_edges = np.histogram(np.log(image[n_mask]), bins=n_bins) histograms.append(hist) diff --git a/package/lib/reduction.py b/package/lib/reduction.py index 7bd86f5..942a27f 100755 --- a/package/lib/reduction.py +++ b/package/lib/reduction.py @@ -446,9 +446,9 @@ def get_error( return_background=False, ): """ - Look for sub-image of shape sub_shape that have the smallest integrated - flux (no source assumption) and define the background on the image by the - standard deviation on this sub-image. + Estimate background intensity level from either fitting the intensity histogram + or by looking for the sub-image of smallest integrated intensity (no source assumption) + and define the background on the image by the standard deviation on this sub-image. ---------- Inputs: data_array : numpy.ndarray @@ -468,7 +468,7 @@ def get_error( If 'auto', look for optimal binning and fit intensity histogram with au gaussian. If str or None, statistic rule to be used for the number of bins in counts/s. If int, number of bins for the counts/s histogram. - If tuple, shape of the sub-image to look for. Must be odd. + If tuple, shape of the sub-image of lowest intensity to look for. Defaults to None. subtract_error : float or bool, optional If float, factor to which the estimated background should be multiplied @@ -1801,8 +1801,6 @@ def rotate_data(data_array, error_array, data_mask, headers): Updated list of headers corresponding to the reduced images accounting for the new orientation angle. """ - # Rotate I_stokes, Q_stokes, U_stokes using rotation matrix - old_center = np.array(data_array[0].shape) / 2 shape = np.fix(np.array(data_array[0].shape) * np.sqrt(2.5)).astype(int) new_center = np.array(shape) / 2 From fce787df39f435819d86bb00b1205440bc577d14 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Mon, 17 Mar 2025 14:33:47 +0100 Subject: [PATCH 6/7] fix rotation of WCS --- package/lib/reduction.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/package/lib/reduction.py b/package/lib/reduction.py index 942a27f..b663c3d 100755 --- a/package/lib/reduction.py +++ b/package/lib/reduction.py @@ -676,6 +676,7 @@ def rebin_array(data_array, error_array, headers, pxsize=2, scale="px", operatio nw.wcs.crpix /= Dxy nw.array_shape = new_shape new_header["NAXIS1"], new_header["NAXIS2"] = nw.array_shape + new_header["PXAREA"] *= Dxy[0] * Dxy[1] for key, val in nw.to_header().items(): new_header.set(key, val) new_header["SAMPLING"] = (str(pxsize) + scale, "Resampling performed during reduction") @@ -1725,10 +1726,6 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, header_st new_wcs.wcs.set() for key, val in new_wcs.to_header().items(): new_header_stokes.set(key, val) - if new_wcs.wcs.pc[0, 0] == 1.0: - new_header_stokes.set("PC1_1", 1.0) - if new_wcs.wcs.pc[1, 1] == 1.0: - new_header_stokes.set("PC2_2", 1.0) new_header_stokes["ORIENTAT"] += ang # Nan handling : From b577fc5afe5f6b8124e93f6fa3052bdad4835cbf Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Tue, 1 Apr 2025 17:02:31 +0200 Subject: [PATCH 7/7] fix WCS computation, cdelt should not be sorted --- package/FOC_reduction.py | 12 ++++++------ package/lib/fits.py | 6 +++--- package/lib/utils.py | 2 +- 3 files changed, 10 insertions(+), 10 deletions(-) diff --git a/package/FOC_reduction.py b/package/FOC_reduction.py index eab1da7..55bd806 100755 --- a/package/FOC_reduction.py +++ b/package/FOC_reduction.py @@ -45,8 +45,8 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= display_bkg = True # Data binning - pxsize = 0.05 - pxscale = "arcsec" # pixel, arcsec or full + pxsize = 4 + pxscale = "px" # pixel, arcsec or full rebin_operation = "sum" # sum or average # Alignement @@ -59,15 +59,15 @@ 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 = 0.075 # If None, no smoothing is done - smoothing_scale = "arcsec" # pixel or arcsec + smoothing_FWHM = 1.5 # If None, no smoothing is done + smoothing_scale = "px" # pixel or arcsec # Rotation rotate_North = True # Polarization map output - P_cut = 5 # if >=1.0 cut on the signal-to-noise else cut on the confidence level in Q, U - SNRi_cut = 5.0 # I measurments with SNR>30, which implies an uncertainty in P of 4.7%. + 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 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 diff --git a/package/lib/fits.py b/package/lib/fits.py index 89177c1..af8c10b 100755 --- a/package/lib/fits.py +++ b/package/lib/fits.py @@ -16,7 +16,7 @@ from astropy.io import fits from astropy.wcs import WCS from .convex_hull import clean_ROI -from .utils import wcs_PA +from .utils import wcs_PA, princ_angle def get_obs_data(infiles, data_folder="", compute_flux=False): @@ -72,13 +72,13 @@ def get_obs_data(infiles, data_folder="", compute_flux=False): for key in keys: header.remove(key, ignore_missing=True) new_cdelt = np.linalg.eigvals(wcs.wcs.cd) - new_cdelt.sort() + # 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"] + header["ORIENTAT"] = princ_angle(float(header["ORIENTAT"])) except KeyError: header["ORIENTAT"] = wcs_PA(new_wcs.wcs.pc[1, 0], np.diag(new_wcs.wcs.pc).mean()) diff --git a/package/lib/utils.py b/package/lib/utils.py index ee9fc5b..9e36dbc 100755 --- a/package/lib/utils.py +++ b/package/lib/utils.py @@ -177,4 +177,4 @@ def wcs_PA(PC21, PC22): orient = np.arccos(PC22) * 180.0 / np.pi elif (abs(PC21) < abs(PC22)) and (PC22 < 0): orient = -np.arccos(PC22) * 180.0 / np.pi - return orient + return princ_angle(orient)