modifications to add raw image display

This commit is contained in:
2025-04-07 17:20:00 +02:00
parent 4b104d79d4
commit cebf4df8a4
4 changed files with 21 additions and 18 deletions

View File

@@ -46,8 +46,8 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
display_bkg = False display_bkg = False
# Data binning # Data binning
pxsize = 4 pxsize = 0.10
pxscale = "px" # pixel, arcsec or full pxscale = "arcsec" # pixel, arcsec or full
rebin_operation = "sum" # sum or average rebin_operation = "sum" # sum or average
# Alignement # Alignement
@@ -60,8 +60,8 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
# Smoothing # Smoothing
smoothing_function = "combine" # gaussian_after, weighted_gaussian_after, gaussian, weighted_gaussian or combine smoothing_function = "combine" # gaussian_after, weighted_gaussian_after, gaussian, weighted_gaussian or combine
smoothing_FWHM = 1.5 # If None, no smoothing is done smoothing_FWHM = 0.150 # If None, no smoothing is done
smoothing_scale = "px" # pixel or arcsec smoothing_scale = "arcsec" # pixel or arcsec
# Rotation # Rotation
rotate_data = False rotate_data = False
@@ -71,7 +71,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 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%. 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 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 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 # Pipeline start

View File

@@ -1,2 +1,2 @@
from .lib import * from .lib import *
from . import FOC_reduction from .FOC_reduction import main

View File

@@ -18,7 +18,6 @@ import matplotlib.dates as mdates
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import numpy as np import numpy as np
from astropy.time import Time from astropy.time import Time
from lib.plots import plot_obs
from matplotlib.colors import LogNorm from matplotlib.colors import LogNorm
from matplotlib.patches import Rectangle from matplotlib.patches import Rectangle
from scipy.optimize import curve_fit from scipy.optimize import curve_fit
@@ -136,6 +135,8 @@ def display_bkg(data, background, std_bkg, headers, histograms=None, binning=Non
fig2.subplots_adjust(hspace=0, wspace=0, right=1.0) fig2.subplots_adjust(hspace=0, wspace=0, right=1.0)
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}$]") 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}$]")
from .plots import plot_obs
if savename is not None: if savename is not None:
this_savename = deepcopy(savename) this_savename = deepcopy(savename)
if savename[-4:] not in [".png", ".jpg", ".pdf"]: if savename[-4:] not in [".png", ".jpg", ".pdf"]:

View File

@@ -440,7 +440,7 @@ def polarization_map(
im = ax.imshow( im = ax.imshow(
imflux, transform=ax.get_transform(WCS(Stokes[0].header).celestial), norm=LogNorm(vmin, vmax), aspect="equal", cmap=kwargs["cmap"], alpha=1.0 imflux, transform=ax.get_transform(WCS(Stokes[0].header).celestial), 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$]") fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]")
levelsF = np.array([0.8, 2.0, 5.0, 10.0, 20.0, 50.0]) / 100.0 * vmax levelsF = np.array([0.8, 2.0, 5.0, 10.0, 20.0, 50.0]) / 100.0 * vmax
print("Flux density contour levels : ", levelsF) print("Flux density contour levels : ", levelsF)
ax.contour(flux, levels=levelsF, transform=ax.get_transform(WCS(Stokes[0].header).celestial), colors="grey", linewidths=0.5) ax.contour(flux, levels=levelsF, transform=ax.get_transform(WCS(Stokes[0].header).celestial), colors="grey", linewidths=0.5)
@@ -457,7 +457,7 @@ def polarization_map(
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) 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) im = ax.imshow(stkI * convert_flux, norm=LogNorm(vmin, vmax), aspect="equal", cmap=kwargs["cmap"], alpha=1.0)
ax.plot(*WCS(Stokes[1]).wcs.crpix, "g+") ax.plot(*WCS(Stokes[1]).wcs.crpix, "g+")
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}$]") fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, 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 # levelsI = np.array([0.8, 2.0, 5.0, 10.0, 20.0, 50.0]) / 100.0 * vmax
# print("Stokes I contour levels : ", levelsI) # print("Stokes I contour levels : ", levelsI)
# ax.contour(stkI * convert_flux, levels=levelsI, colors="grey", linewidths=0.5) # ax.contour(stkI * convert_flux, levels=levelsI, colors="grey", linewidths=0.5)
@@ -476,7 +476,9 @@ def polarization_map(
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) 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() 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) 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}$]") fig.colorbar(
im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, 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.linspace(0.0.60, 0.50, 5) * pfmax
levelsPf = np.array([1.73, 13.0, 33.0, 66.0]) / 100.0 * pfmax levelsPf = np.array([1.73, 13.0, 33.0, 66.0]) / 100.0 * pfmax
print("Polarized flux density contour levels : ", levelsPf) print("Polarized flux density contour levels : ", levelsPf)
@@ -486,13 +488,13 @@ def polarization_map(
display = "p" display = "p"
vmin, vmax = 0.0, 100.0 vmin, vmax = 0.0, 100.0
im = ax.imshow(pol * 100.0, vmin=vmin, vmax=vmax, aspect="equal", cmap=kwargs["cmap"], alpha=1.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$ [%]") fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$P$ [%]")
elif display.lower() in ["pa", "pang", "pol_ang"]: elif display.lower() in ["pa", "pang", "pol_ang"]:
# Display polarization degree map # Display polarization degree map
display = "pa" display = "pa"
vmin, vmax = 0.0, 180.0 vmin, vmax = 0.0, 180.0
im = ax.imshow(princ_angle(pang), vmin=vmin, vmax=vmax, aspect="equal", cmap=kwargs["cmap"], alpha=1.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=50, shrink=0.60, pad=0.025, label=r"$\theta_P$ [°]") fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$\theta_P$ [°]")
elif display.lower() in ["s_p", "pol_err", "pol_deg_err"]: elif display.lower() in ["s_p", "pol_err", "pol_deg_err"]:
# Display polarization degree error map # Display polarization degree error map
display = "s_p" display = "s_p"
@@ -502,7 +504,7 @@ def polarization_map(
else: else:
vmin, vmax = 0.0, 100.0 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) 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$ [%]") fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$\sigma_P$ [%]")
elif display.lower() in ["s_i", "i_err"]: elif display.lower() in ["s_i", "i_err"]:
# Display intensity error map # Display intensity error map
display = "s_i" display = "s_i"
@@ -514,7 +516,7 @@ def polarization_map(
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: else:
im = ax.imshow(np.sqrt(stk_cov[0, 0]) * convert_flux, aspect="equal", cmap=kwargs["cmap"], alpha=1.0) 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}$]") fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$\sigma_I$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]")
elif display.lower() in ["snri"]: elif display.lower() in ["snri"]:
# Display I_stokes signal-to-noise map # Display I_stokes signal-to-noise map
display = "snri" display = "snri"
@@ -526,7 +528,7 @@ def polarization_map(
ax.contour(SNRi, levels=levelsSNRi, colors="grey", linewidths=0.5) ax.contour(SNRi, levels=levelsSNRi, colors="grey", linewidths=0.5)
else: else:
im = ax.imshow(SNRi, aspect="equal", cmap=kwargs["cmap"], alpha=1.0) im = ax.imshow(SNRi, aspect="equal", cmap=kwargs["cmap"], alpha=1.0)
fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$I_{Stokes}/\sigma_{I}$") fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$I_{Stokes}/\sigma_{I}$")
elif display.lower() in ["snr", "snrp"]: elif display.lower() in ["snr", "snrp"]:
# Display polarization degree signal-to-noise map # Display polarization degree signal-to-noise map
display = "snrp" display = "snrp"
@@ -538,7 +540,7 @@ def polarization_map(
ax.contour(SNRp, levels=levelsSNRp, colors="grey", linewidths=0.5) ax.contour(SNRp, levels=levelsSNRp, colors="grey", linewidths=0.5)
else: else:
im = ax.imshow(SNRp, aspect="equal", cmap=kwargs["cmap"], alpha=1.0) im = ax.imshow(SNRp, aspect="equal", cmap=kwargs["cmap"], alpha=1.0)
fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$P/\sigma_{P}$") fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$P/\sigma_{P}$")
elif display.lower() in ["conf", "confp"]: elif display.lower() in ["conf", "confp"]:
# Display polarization degree signal-to-noise map # Display polarization degree signal-to-noise map
display = "confp" display = "confp"
@@ -547,7 +549,7 @@ def polarization_map(
levelsconfp = np.array([0.500, 0.900, 0.990, 0.999]) levelsconfp = np.array([0.500, 0.900, 0.990, 0.999])
print("confp contour levels : ", levelsconfp) print("confp contour levels : ", levelsconfp)
ax.contour(confP, levels=levelsconfp, colors="grey", linewidths=0.5) ax.contour(confP, levels=levelsconfp, colors="grey", linewidths=0.5)
fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$Conf_{P}$") fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$Conf_{P}$")
else: else:
# Defaults to intensity map # Defaults to intensity map
if flux_lim is not None: if flux_lim is not None:
@@ -562,7 +564,7 @@ def polarization_map(
cmap=kwargs["cmap"], cmap=kwargs["cmap"],
alpha=1.0, 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$]") fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.015, fraction=0.03, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]")
I_diluted, I_diluted_err = np.sum(flux[flux_mask]), np.sqrt(np.sum(flux_error[flux_mask] ** 2)) I_diluted, I_diluted_err = np.sum(flux[flux_mask]), np.sqrt(np.sum(flux_error[flux_mask] ** 2))
plt.rcParams.update({"font.size": 11}) plt.rcParams.update({"font.size": 11})