Tibeuleu/FOC_Reduction
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

save raw flux in fits file and display

Browse Source 
fix rebase display on main

rebase display on main
This commit is contained in:
Thibault Barnouin
2025-03-19 16:38:36 +01:00
parent c94d646b02
commit 3da6cbcfaf
4 changed files with 201 additions and 112 deletions
Download Patch File Download Diff File Expand all files Collapse all files

36
package/FOC_reduction.py
View File

@@ -20,6 +20,7 @@ import lib.reduction as proj_red # Functions used in reduction pipeline
import numpy as np
from lib.utils import princ_angle, sci_not
from matplotlib.colors import LogNorm
from astropy.wcs import WCS
def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=False, interactive=False):
@@ -63,6 +64,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
smoothing_scale = "arcsec" # pixel or arcsec
# Rotation
rotate_data = False
rotate_North = True
# Polarization map output
@@ -142,7 +144,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
)
# Rotate data to have same orientation
rotate_data = np.unique([np.round(float(head["ORIENTAT"]), 3) for head in headers]).size != 1
rotate_data = rotate_data or np.unique([np.round(float(head["ORIENTAT"]), 3) for head in headers]).size != 1
if rotate_data:
ang = np.mean([head["ORIENTAT"] for head in headers])
for head in headers:
@@ -158,7 +160,6 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
vmin=data_array[data_array > 0.0].min() * headers[0]["photflam"], vmax=data_array[data_array > 0.0].max() * headers[0]["photflam"]
),
)
# Align and rescale images with oversampling.
data_array, error_array, headers, data_mask, shifts, error_shifts = proj_red.align_data(
data_array,
@@ -243,6 +244,8 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
)
flux_data, flux_error, flux_mask, flux_head = proj_red.rotate_data(np.array([flux_data]), np.array([flux_error]), flux_mask, [flux_head])
flux_data, flux_error, flux_head = flux_data[0], flux_error[0], flux_head[0]
elif not rotate_data:
figname += "_noROT"
# Compute polarimetric parameters (polarization degree and angle).
P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P = proj_red.compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, header_stokes)
@@ -250,7 +253,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
# Step 4:
# Save image to FITS.
figname = "_".join([figname, figtype]) if figtype != "" else figname
savename = "_".join([figname, figtype]) if figtype != "" else figname
Stokes_hdul = proj_fits.save_Stokes(
I_stokes,
Q_stokes,
@@ -265,18 +268,18 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
s_PA_P,
header_stokes,
data_mask,
figname,
savename,
data_folder=data_folder,
return_hdul=True,
flux_data=flux_data,
flux_data=np.array([flux_data, flux_error, flux_mask]),
flux_head=flux_head,
)
outfiles.append("/".join([data_folder, Stokes_hdul["I_STOKES"].header["FILENAME"] + ".fits"]))
outfiles.append("/".join([data_folder, Stokes_hdul[0].header["FILENAME"] + ".fits"]))
# Step 5:
# crop to desired region of interest (roi)
if crop:
figname += "_crop"
savename += "_crop"
stokescrop = proj_plots.crop_Stokes(deepcopy(Stokes_hdul), norm=LogNorm())
stokescrop.crop()
stokescrop.write_to("/".join([data_folder, figname + ".fits"]))
@@ -286,13 +289,8 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
data_mask = Stokes_hdul["data_mask"].data.astype(bool)
print(
"F_int({0:.0f} Angs) = ({1} ± {2})e{3} ergs.cm^-2.s^-1.Angs^-1".format(
header_stokes["PHOTPLAM"],
*sci_not(
Stokes_hdul["I_STOKES"].data[data_mask].sum() * header_stokes["PHOTFLAM"],
np.sqrt(Stokes_hdul["IQU_COV_MATRIX"].data[0, 0][data_mask].sum()) * header_stokes["PHOTFLAM"],
2,
out=int,
),
flux_head["PHOTPLAM"],
*sci_not(flux_data[flux_mask].sum() * flux_head["PHOTFLAM"], np.sqrt(np.sum(flux_error[flux_mask] ** 2)) * flux_head["PHOTFLAM"], 2, out=int),
)
)
print("P_int = {0:.1f} ± {1:.1f} %".format(header_stokes["p_int"] * 100.0, np.ceil(header_stokes["sP_int"] * 1000.0) / 10.0))
@@ -315,12 +313,12 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
flux_lim=flux_lim,
step_vec=step_vec,
scale_vec=scale_vec,
savename="_".join([figname]),
savename=figname,
plots_folder=plots_folder,
)
for figtype, figsuffix in zip(
["Intensity", "Pol_flux", "Pol_deg", "Pol_ang", "I_err", "P_err", "SNRi", "SNRp", "confp"],
["I", "P_flux", "P", "PA", "I_err", "P_err", "SNRi", "SNRp", "confP"],
["FluxDensity", "Intensity", "Pol_flux", "Pol_deg", "Pol_ang", "I_err", "P_err", "SNRi", "SNRp", "confp"],
["F", "I", "P_flux", "P", "PA", "I_err", "P_err", "SNRi", "SNRp", "confP"],
):
try:
proj_plots.polarization_map(
@@ -331,7 +329,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
flux_lim=flux_lim,
step_vec=step_vec,
scale_vec=scale_vec,
savename="_".join([figname, figsuffix]),
savename="_".join([savename, figsuffix]),
plots_folder=plots_folder,
display=figtype,
)
@@ -339,7 +337,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
pass
elif not interactive:
proj_plots.polarization_map(
deepcopy(Stokes_hdul), data_mask, P_cut=P_cut, SNRi_cut=SNRi_cut, savename=figname, plots_folder=plots_folder, display="integrate"
deepcopy(Stokes_hdul), data_mask, P_cut=P_cut, SNRi_cut=SNRi_cut, savename=savename, plots_folder=plots_folder, display="integrate"
)
elif pxscale.lower() not in ["full", "integrate"]:
proj_plots.pol_map(Stokes_hdul, P_cut=P_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, scale_vec=scale_vec, flux_lim=flux_lim)