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

update scripts and plots to new fits data shape

Browse Source
This commit is contained in:
Thibault Barnouin
2025-09-04 11:17:14 +02:00
parent 8a8359fed0
commit f05f6b789c
5 changed files with 114 additions and 103 deletions
Download Patch File Download Diff File Expand all files Collapse all files

44
package/lib/plots.py
View File

@@ -181,9 +181,9 @@ def plot_Stokes(Stokes, savename=None, plots_folder=""):
Defaults to current folder.
"""
# Get data
stkI = Stokes["I_stokes"].data.copy()
stkQ = Stokes["Q_stokes"].data.copy()
stkU = Stokes["U_stokes"].data.copy()
stkI = Stokes["stokes"].data[0].copy()
stkQ = Stokes["stokes"].data[1].copy()
stkU = Stokes["stokes"].data[2].copy()
data_mask = Stokes["Data_mask"].data.astype(bool)
for dataset in [stkI, stkQ, stkU]:
@@ -288,10 +288,10 @@ def polarization_map(
The figure and ax created for interactive contour maps.
"""
# Get data
stkI = Stokes["I_stokes"].data.copy()
stkQ = Stokes["Q_stokes"].data.copy()
stkU = Stokes["U_stokes"].data.copy()
stk_cov = Stokes["IQU_COV_MATRIX"].data.copy()
stkI = Stokes["STOKES"].data[0].copy()
stkQ = Stokes["STOKES"].data[1].copy()
stkU = Stokes["STOKES"].data[2].copy()
stk_cov = Stokes["STOKES_COV"].data.copy()
pol = Stokes["Pol_deg_debiased"].data.copy()
pol_err = Stokes["Pol_deg_err"].data.copy()
pang = Stokes["Pol_ang"].data.copy()
@@ -314,7 +314,7 @@ def polarization_map(
for j in range(3):
stk_cov[i][j][np.logical_not(data_mask)] = np.nan
wcs = WCS(Stokes[0]).deepcopy()
wcs = WCS(Stokes[0]).celestial.deepcopy()
pivot_wav = Stokes[0].header["photplam"]
convert_flux = Stokes[0].header["photflam"]
@@ -968,10 +968,10 @@ class overplot_radio(align_maps):
self.wcs_UV = self.map_wcs
# Get Data
obj = self.Stokes_UV[0].header["targname"]
stkI = self.Stokes_UV["I_STOKES"].data
stkQ = self.Stokes_UV["Q_STOKES"].data
stkU = self.Stokes_UV["U_STOKES"].data
stk_cov = self.Stokes_UV["IQU_COV_MATRIX"].data
stkI = self.Stokes_UV["STOKES"].data[0]
stkQ = self.Stokes_UV["STOKES"].data[1]
stkU = self.Stokes_UV["STOKES"].data[2]
stk_cov = self.Stokes_UV["STOKES_COV"].data
pol = deepcopy(self.Stokes_UV["POL_DEG_DEBIASED"].data)
pol_err = self.Stokes_UV["POL_DEG_ERR"].data
pang = self.Stokes_UV["POL_ANG"].data
@@ -1187,10 +1187,10 @@ class overplot_chandra(align_maps):
self.wcs_UV = self.map_wcs
# Get Data
obj = self.Stokes_UV[0].header["targname"]
stkI = self.Stokes_UV["I_STOKES"].data
stkQ = self.Stokes_UV["Q_STOKES"].data
stkU = self.Stokes_UV["U_STOKES"].data
stk_cov = self.Stokes_UV["IQU_COV_MATRIX"].data
stkI = self.Stokes_UV["STOKES"].data[0]
stkQ = self.Stokes_UV["STOKES"].data[1]
stkU = self.Stokes_UV["STOKES"].data[2]
stk_cov = self.Stokes_UV["STOKES_COV"].data
pol = deepcopy(self.Stokes_UV["POL_DEG_DEBIASED"].data)
pol_err = self.Stokes_UV["POL_DEG_ERR"].data
pang = self.Stokes_UV["POL_ANG"].data
@@ -1398,10 +1398,10 @@ class overplot_pol(align_maps):
self.wcs_UV = self.map_wcs
# Get Data
obj = self.Stokes_UV[0].header["targname"]
stkI = self.Stokes_UV["I_STOKES"].data
stkQ = self.Stokes_UV["Q_STOKES"].data
stkU = self.Stokes_UV["U_STOKES"].data
stk_cov = self.Stokes_UV["IQU_COV_MATRIX"].data
stkI = self.Stokes_UV["STOKES"].data[0]
stkQ = self.Stokes_UV["STOKES"].data[1]
stkU = self.Stokes_UV["STOKES"].data[2]
stk_cov = self.Stokes_UV["STOKES_COV"].data
pol = deepcopy(self.Stokes_UV["POL_DEG_DEBIASED"].data)
pol_err = self.Stokes_UV["POL_DEG_ERR"].data
pang = self.Stokes_UV["POL_ANG"].data
@@ -1702,8 +1702,8 @@ class align_pol(object):
def single_plot(self, curr_map, wcs, v_lim=None, ax_lim=None, P_cut=3.0, SNRi_cut=3.0, savename=None, **kwargs):
# Get data
stkI = curr_map["I_STOKES"].data
stk_cov = curr_map["IQU_COV_MATRIX"].data
stkI = curr_map["STOKES"].data[0]
stk_cov = curr_map["STOKES_COV"].data
pol = deepcopy(curr_map["POL_DEG_DEBIASED"].data)
pol_err = curr_map["POL_DEG_ERR"].data
pang = curr_map["POL_ANG"].data

109
package/src/Combine.py
View File

@@ -1,5 +1,6 @@
#!/usr/bin/env python3
# -*- coding:utf-8 -*-
from copy import deepcopy
from pathlib import Path
from sys import path as syspath
@@ -26,7 +27,7 @@ def same_reduction(infiles):
except KeyError:
pass
test_IQU = True
for look in ["I_stokes", "Q_stokes", "U_stokes", "IQU_cov_matrix"]:
for look in ["STOKES", "STOKES_COV"]:
test_IQU *= look in datatype
params["IQU"].append(test_IQU)
# test for orientation and pixel size
@@ -88,73 +89,78 @@ def same_obs(infiles, data_folder):
def combine_Stokes(infiles):
"""
Combine I, Q, U from different observations of a same object.
Combine Stokes matrices from different observations of a same object.
"""
from astropy.io.fits import open as fits_open
from lib.reduction import align_data, zeropad
from lib.utils import remove_stokes_axis_from_header
from scipy.ndimage import shift as sc_shift
I_array, Q_array, U_array, IQU_cov_array, data_mask, headers = [], [], [], [], [], []
Stokes_array, Stokes_cov_array, Stokes_cov_stat_array, data_mask, headers = [], [], [], [], []
shape = np.array([0, 0])
for file in infiles:
with fits_open(file) as f:
headers.append(f[0].header)
I_array.append(f["I_stokes"].data)
Q_array.append(f["Q_stokes"].data)
U_array.append(f["U_stokes"].data)
IQU_cov_array.append(f["IQU_cov_matrix"].data)
Stokes_array.append(f["stokes"].data)
Stokes_cov_array.append(f["stokes_cov"].data)
Stokes_cov_stat_array.append(f["stokes_cov_stat"].data)
data_mask.append(f["data_mask"].data.astype(bool))
shape[0] = np.max([shape[0], f["I_stokes"].data.shape[0]])
shape[1] = np.max([shape[1], f["I_stokes"].data.shape[1]])
shape[0] = np.max([shape[0], f["stokes"].data[0].shape[0]])
shape[1] = np.max([shape[1], f["stokes"].data[0].shape[1]])
exposure_array = np.array([float(head["EXPTIME"]) for head in headers])
shape += np.array([5, 5])
data_mask = np.sum([zeropad(mask, shape) for mask in data_mask], axis=0).astype(bool)
I_array = np.array([zeropad(I, shape) for I in I_array])
Q_array = np.array([zeropad(Q, shape) for Q in Q_array])
U_array = np.array([zeropad(U, shape) for U in U_array])
IQU_cov_array = np.array([[[zeropad(cov[i, j], shape) for j in range(3)] for i in range(3)] for cov in IQU_cov_array])
Stokes_array = np.array([[zeropad(stk[i], shape) for i in range(4)] for stk in Stokes_array])
Stokes_cov_array = np.array([[[zeropad(cov[i, j], shape) for j in range(4)] for i in range(4)] for cov in Stokes_cov_array])
Stokes_cov_stat_array = np.array([[[zeropad(cov_stat[i, j], shape) for j in range(4)] for i in range(4)] for cov_stat in Stokes_cov_stat_array])
sI_array = np.sqrt(IQU_cov_array[:, 0, 0])
sQ_array = np.sqrt(IQU_cov_array[:, 1, 1])
sU_array = np.sqrt(IQU_cov_array[:, 2, 2])
I_array = deepcopy(Stokes_array[:, 0])
sI_array = deepcopy(np.sqrt(Stokes_cov_array[:, 0, 0]))
_, _, _, _, shifts, errors = align_data(I_array, headers, error_array=sI_array, data_mask=data_mask, ref_center="center", return_shifts=True)
heads = [remove_stokes_axis_from_header(head) for head in headers]
_, _, _, _, shifts, errors = align_data(
I_array, heads, error_array=sI_array, background=sI_array[:, 0, 0], data_mask=data_mask, ref_center="center", return_shifts=True
)
data_mask_aligned = np.sum([sc_shift(data_mask, s, order=1, cval=0.0) for s in shifts], axis=0).astype(bool)
I_aligned, sI_aligned = (
np.array([sc_shift(I, s, order=1, cval=0.0) for I, s in zip(I_array, shifts)]),
np.array([sc_shift(sI, s, order=1, cval=0.0) for sI, s in zip(sI_array, shifts)]),
Stokes_aligned = np.array([[sc_shift(stk[i], s, order=1, cval=0.0) for i in range(4)] for stk, s in zip(Stokes_array, shifts)])
Stokes_cov_aligned = np.array(
[[[sc_shift(cov[i, j], s, order=1, cval=0.0) for j in range(4)] for i in range(4)] for cov, s in zip(Stokes_cov_array, shifts)]
)
Q_aligned, sQ_aligned = (
np.array([sc_shift(Q, s, order=1, cval=0.0) for Q, s in zip(Q_array, shifts)]),
np.array([sc_shift(sQ, s, order=1, cval=0.0) for sQ, s in zip(sQ_array, shifts)]),
Stokes_cov_stat_aligned = np.array(
[[[sc_shift(cov_stat[i, j], s, order=1, cval=0.0) for j in range(4)] for i in range(4)] for cov_stat, s in zip(Stokes_cov_stat_array, shifts)]
)
U_aligned, sU_aligned = (
np.array([sc_shift(U, s, order=1, cval=0.0) for U, s in zip(U_array, shifts)]),
np.array([sc_shift(sU, s, order=1, cval=0.0) for sU, s in zip(sU_array, shifts)]),
)
IQU_cov_aligned = np.array([[[sc_shift(cov[i, j], s, order=1, cval=0.0) for j in range(3)] for i in range(3)] for cov, s in zip(IQU_cov_array, shifts)])
I_combined = np.sum([exp * I for exp, I in zip(exposure_array, I_aligned)], axis=0) / exposure_array.sum()
Q_combined = np.sum([exp * Q for exp, Q in zip(exposure_array, Q_aligned)], axis=0) / exposure_array.sum()
U_combined = np.sum([exp * U for exp, U in zip(exposure_array, U_aligned)], axis=0) / exposure_array.sum()
Stokes_combined = np.zeros((4, shape[0], shape[1]))
for i in range(4):
Stokes_combined[i] = np.sum([exp * stk for exp, stk in zip(exposure_array, Stokes_aligned[:, i])], axis=0) / exposure_array.sum()
IQU_cov_combined = np.zeros((3, 3, shape[0], shape[1]))
for i in range(3):
IQU_cov_combined[i, i] = np.sum([exp**2 * cov for exp, cov in zip(exposure_array, IQU_cov_aligned[:, i, i])], axis=0) / exposure_array.sum() ** 2
for j in [x for x in range(3) if x != i]:
IQU_cov_combined[i, j] = np.sqrt(
np.sum([exp**2 * cov**2 for exp, cov in zip(exposure_array, IQU_cov_aligned[:, i, j])], axis=0) / exposure_array.sum() ** 2
Stokes_cov_combined = np.zeros((4, 4, shape[0], shape[1]))
Stokes_cov_stat_combined = np.zeros((4, 4, shape[0], shape[1]))
for i in range(4):
Stokes_cov_combined[i, i] = np.sum([exp**2 * cov for exp, cov in zip(exposure_array, Stokes_cov_aligned[:, i, i])], axis=0) / exposure_array.sum() ** 2
Stokes_cov_stat_combined[i, i] = (
np.sum([exp**2 * cov_stat for exp, cov_stat in zip(exposure_array, Stokes_cov_stat_aligned[:, i, i])], axis=0) / exposure_array.sum() ** 2
)
for j in [x for x in range(4) if x != i]:
Stokes_cov_combined[i, j] = np.sqrt(
np.sum([exp**2 * cov**2 for exp, cov in zip(exposure_array, Stokes_cov_aligned[:, i, j])], axis=0) / exposure_array.sum() ** 2
)
IQU_cov_combined[j, i] = np.sqrt(
np.sum([exp**2 * cov**2 for exp, cov in zip(exposure_array, IQU_cov_aligned[:, j, i])], axis=0) / exposure_array.sum() ** 2
Stokes_cov_combined[j, i] = np.sqrt(
np.sum([exp**2 * cov**2 for exp, cov in zip(exposure_array, Stokes_cov_aligned[:, j, i])], axis=0) / exposure_array.sum() ** 2
)
Stokes_cov_stat_combined[i, j] = np.sqrt(
np.sum([exp**2 * cov_stat**2 for exp, cov_stat in zip(exposure_array, Stokes_cov_stat_aligned[:, i, j])], axis=0) / exposure_array.sum() ** 2
)
Stokes_cov_stat_combined[j, i] = np.sqrt(
np.sum([exp**2 * cov_stat**2 for exp, cov_stat in zip(exposure_array, Stokes_cov_stat_aligned[:, j, i])], axis=0) / exposure_array.sum() ** 2
)
header_combined = headers[0]
header_combined["EXPTIME"] = exposure_array.sum()
return I_combined, Q_combined, U_combined, IQU_cov_combined, data_mask_aligned, header_combined
return Stokes_combined, Stokes_cov_combined, Stokes_cov_stat_combined, data_mask_aligned, header_combined
def main(infiles, target=None, output_dir="./data/"):
@@ -190,21 +196,24 @@ def main(infiles, target=None, output_dir="./data/"):
infiles = new_infiles
I_combined, Q_combined, U_combined, IQU_cov_combined, data_mask_combined, header_combined = combine_Stokes(infiles=infiles)
I_combined, Q_combined, U_combined, IQU_cov_combined, data_mask_combined, header_combined = rotate_Stokes(
I_stokes=I_combined, Q_stokes=Q_combined, U_stokes=U_combined, Stokes_cov=IQU_cov_combined, data_mask=data_mask_combined, header_stokes=header_combined
Stokes_combined, Stokes_cov_combined, Stokes_cov_stat_combined, data_mask_combined, header_combined = combine_Stokes(infiles=infiles)
Stokes_combined, Stokes_cov_combined, data_mask_combined, header_combined, Stokes_cov_stat_combined = rotate_Stokes(
Stokes=Stokes_combined,
Stokes_cov=Stokes_cov_combined,
Stokes_cov_stat=Stokes_cov_stat_combined,
data_mask=data_mask_combined,
header_stokes=header_combined,
)
P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P = compute_pol(
I_stokes=I_combined, Q_stokes=Q_combined, U_stokes=U_combined, Stokes_cov=IQU_cov_combined, header_stokes=header_combined
Stokes=Stokes_combined, Stokes_cov=Stokes_cov_combined, Stokes_cov_stat=Stokes_cov_stat_combined, header_stokes=header_combined
)
filename = header_combined["FILENAME"]
figname = "_".join([target, filename[filename.find("FOC_") :], "combined"])
Stokes_combined = save_Stokes(
I_stokes=I_combined,
Q_stokes=Q_combined,
U_stokes=U_combined,
Stokes_cov=IQU_cov_combined,
Stokes_c = save_Stokes(
Stokes=Stokes_combined,
Stokes_cov=Stokes_cov_combined,
Stokes_cov_stat=Stokes_cov_stat_combined,
P=P,
debiased_P=debiased_P,
s_P=s_P,
@@ -219,7 +228,7 @@ def main(infiles, target=None, output_dir="./data/"):
return_hdul=True,
)
pol_map(Stokes_combined, **kwargs)
pol_map(Stokes_c, **kwargs)
return "/".join([data_folder, figname + ".fits"])

8
package/src/emission_center.py
View File

@@ -20,13 +20,13 @@ def main(infile, P_cut=0.99, target=None, display="pf", output_dir=None):
output = []
Stokes = fits_open(infile)
stkI = Stokes["I_STOKES"].data
s_I = np.sqrt(Stokes["IQU_COV_MATRIX"].data[0, 0])
stkI = Stokes["STOKES"].data[0]
s_I = np.sqrt(Stokes["STOKES_COV"].data[0, 0])
SNRi = np.zeros(stkI.shape)
SNRi[s_I > 0.0] = stkI[s_I > 0.0] / s_I[s_I > 0.0]
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["STOKES"].data[1], Stokes["STOKES"].data[2], np.sqrt(Stokes["STOKES_COV"].data[1, 1]), np.sqrt(Stokes["STOKES_COV"].data[2, 2])],
[QN, UN, QN_ERR, UN_ERR],
):
nflux[stkI > 0.0] = sflux[stkI > 0.0] / stkI[stkI > 0.0]
@@ -41,7 +41,7 @@ def main(infile, P_cut=0.99, target=None, display="pf", output_dir=None):
Stokescentconf, Stokescenter = CenterConf(Stokesconf > P_cut, Stokes["POL_ANG"].data, Stokes["POL_ANG_ERR"].data)
else:
Stokescentconf, Stokescenter = CenterConf(Stokessnr > P_cut, Stokes["POL_ANG"].data, Stokes["POL_ANG_ERR"].data)
Stokespos = WCS(Stokes[0].header).pixel_to_world(*Stokescenter)
Stokespos = WCS(Stokes[0].header).celestial.pixel_to_world(*Stokescenter)
if target is None:
target = Stokes[0].header["TARGNAME"]

10
package/src/overplot_IC5063.py
View File

@@ -10,7 +10,7 @@ from astropy.io import fits
from lib.plots import overplot_pol, overplot_radio
from matplotlib.colors import LogNorm
Stokes_UV = fits.open("./data/IC5063/5918/IC5063_FOC_b0.10arcsec_c0.20arcsec.fits")
Stokes_UV = fits.open("./data/IC5063/5918/IC5063_5918_F502M_FOC_b0.10arcsec_c0.15arcsec.fits")
# Stokes_18GHz = fits.open("./data/IC5063/radio/IC5063_18GHz.fits")
# Stokes_24GHz = fits.open("./data/IC5063/radio/IC5063_24GHz.fits")
# Stokes_103GHz = fits.open("./data/IC5063/radio/IC5063_103GHz.fits")
@@ -47,10 +47,12 @@ Stokes_IR = fits.open("./data/IC5063/IR/u2e65g01t_c0f_rot.fits")
G = overplot_pol(Stokes_UV, Stokes_IR, cmap="inferno")
G.plot(
P_cut=0.99,
SNRi_cut=1.0,
P_cut=3,
SNRi_cut=10,
savename="./plots/IC5063/IR_overplot.pdf",
scale_vec=None,
scale_vec=5,
norm=LogNorm(Stokes_IR[0].data.max() * Stokes_IR[0].header["photflam"] / 1e3, Stokes_IR[0].data.max() * Stokes_IR[0].header["photflam"]),
cmap="inferno",
disptype="pf",
levels="Default",
)

46
package/src/overplot_MRK463E.py
View File

@@ -11,9 +11,9 @@ from lib.plots import overplot_chandra, overplot_pol
from matplotlib.colors import LogNorm
Stokes_UV = fits.open("./data/MRK463E/5960/MRK463E_FOC_b0.05arcsec_c0.10arcsec.fits")
# Stokes_IR = fits.open("./data/MRK463E/WFPC2/IR_rot_crop.fits")
Stokes_IR = fits.open("./data/MRK463E/WFPC2/IR_rot_crop.fits")
# Stokes_Xr = fits.open("./data/MRK463E/Chandra/X_ray_crop.fits")
Radio = fits.open("./data/MRK463E/EMERLIN/Voorwerpjes_1356+1822_1356+1822_uniform-image.fits")
# Radio = fits.open("./data/MRK463E/EMERLIN/Voorwerpjes_1356+1822_1356+1822_uniform-image.fits")
# levels = np.geomspace(1.0, 99.0, 7)
@@ -21,26 +21,26 @@ Radio = fits.open("./data/MRK463E/EMERLIN/Voorwerpjes_1356+1822_1356+1822_unifor
# A.plot(levels=levels, P_cut=0.99, SNRi_cut=1.0, scale_vec=5, zoom=1, savename="./plots/MRK463E/Chandra_overplot.pdf")
# A.write_to(path1="./data/MRK463E/FOC_data_Chandra.fits", path2="./data/MRK463E/Chandra_data.fits", suffix="aligned")
# levels = np.array([0.8, 2, 5, 10, 20, 50]) / 100.0 * Stokes_UV[0].header["photflam"]
# B = overplot_pol(Stokes_UV, Stokes_IR, norm=LogNorm())
# B.plot(levels=levels, P_cut=0.99, SNRi_cut=1.0, scale_vec=5, norm=LogNorm(8.5e-18, 2.5e-15), savename="./plots/MRK463E/IR_overplot.pdf")
# B.write_to(path1="./data/MRK463E/FOC_data_WFPC.fits", path2="./data/MRK463E/WFPC_data.fits", suffix="aligned")
levels = np.array([0.8, 2, 5, 10, 20, 50]) / 100.0 * Stokes_UV[0].header["photflam"]
B = overplot_pol(Stokes_UV, Stokes_IR, norm=LogNorm())
B.plot(levels=levels, P_cut=0.99, SNRi_cut=1.0, scale_vec=5, norm=LogNorm(8.5e-18, 2.5e-15), savename="./plots/MRK463E/IR_overplot.pdf")
B.write_to(path1="./data/MRK463E/FOC_data_WFPC.fits", path2="./data/MRK463E/WFPC_data.fits", suffix="aligned")
# levels = np.array([0.8, 2, 5, 10, 20, 50]) / 100.0 * Stokes_UV[0].header["photflam"]
levels = np.array([5, 10, 20, 50])
C = overplot_pol(Stokes_UV, Radio, norm=LogNorm())
C.other_im.set(norm=LogNorm(1e-4, 2e-2))
C.plot(
levels=levels,
P_cut=0.99,
SNRi_cut=1.0,
step_vec=0,
scale_vec=3,
norm=LogNorm(1e-4, 2e-2),
cmap="inferno_r",
width=0.5,
linewidth=0.5,
disptype="snri",
savename="./plots/MRK463E/EMERLIN_snri_overplot.pdf",
)
C.write_to(path1="./data/MRK463E/FOC_data_EMERLIN.fits", path2="./data/MRK463E/EMERLIN_data.fits", suffix="aligned")
# levels = np.array([5, 10, 20, 50])
# C = overplot_pol(Stokes_UV, Radio, norm=LogNorm())
# C.other_im.set(norm=LogNorm(1e-4, 2e-2))
# C.plot(
# levels=levels,
# P_cut=0.99,
# SNRi_cut=1.0,
# step_vec=0,
# scale_vec=3,
# norm=LogNorm(1e-4, 2e-2),
# cmap="inferno_r",
# width=0.5,
# linewidth=0.5,
# disptype="snri",
# savename="./plots/MRK463E/EMERLIN_snri_overplot.pdf",
# )
# C.write_to(path1="./data/MRK463E/FOC_data_EMERLIN.fits", path2="./data/MRK463E/EMERLIN_data.fits", suffix="aligned")