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

fix plots aspect ratio, fits saving of reduction parameters

Browse Source
This commit is contained in:
Thibault Barnouin
2024-07-08 14:48:18 +02:00
parent a6edf1c1b9
commit d8365e984d
2 changed files with 48 additions and 43 deletions
Download Patch File Download Diff File Expand all files Collapse all files

42
package/lib/fits.py
View File

@@ -16,6 +16,7 @@ from astropy.io import fits
from astropy.wcs import WCS from astropy.wcs import WCS
from .convex_hull import clean_ROI from .convex_hull import clean_ROI
from .utils import wcs_PA
def get_obs_data(infiles, data_folder="", compute_flux=False): def get_obs_data(infiles, data_folder="", compute_flux=False):
@@ -57,22 +58,20 @@ def get_obs_data(infiles, data_folder="", compute_flux=False):
if new_wcs.wcs.has_cd() or (new_wcs.wcs.cdelt[:2] == np.array([1.0, 1.0])).all(): if new_wcs.wcs.has_cd() or (new_wcs.wcs.cdelt[:2] == np.array([1.0, 1.0])).all():
# Update WCS with relevant information # Update WCS with relevant information
if new_wcs.wcs.has_cd(): if new_wcs.wcs.has_cd():
old_cd = new_wcs.wcs.cd
del new_wcs.wcs.cd 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"] 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: for key in keys:
header.remove(key, ignore_missing=True) header.remove(key, ignore_missing=True)
new_cdelt = np.linalg.eig(old_cd)[0] new_cdelt = np.linalg.eigvals(wcs.wcs.cd)
elif (new_wcs.wcs.cdelt == np.array([1.0, 1.0])).all() and (new_wcs.array_shape in [(512, 512), (1024, 512), (512, 1024), (1024, 1024)]): new_cdelt.sort()
old_cd = new_wcs.wcs.pc new_wcs.wcs.pc = wcs.wcs.cd.dot(np.diag(1.0 / new_cdelt))
new_wcs.wcs.pc = np.dot(old_cd, np.diag(1.0 / new_cdelt))
new_wcs.wcs.cdelt = new_cdelt new_wcs.wcs.cdelt = new_cdelt
for key, val in new_wcs.to_header().items(): for key, val in new_wcs.to_header().items():
header[key] = val header[key] = val
try: try:
_ = header["ORIENTAT"] _ = header["ORIENTAT"]
except KeyError: except KeyError:
header["ORIENTAT"] = -np.arccos(new_wcs.wcs.pc[0, 0]) * 180.0 / np.pi header["ORIENTAT"] = wcs_PA(new_wcs.wcs.pc[1, 0], np.diag(new_wcs.wcs.pc).mean())
# force WCS for POL60 to have same pixel size as POL0 and POL120 # force WCS for POL60 to have same pixel size as POL0 and POL120
is_pol60 = np.array([head["filtnam1"].lower() == "pol60" for head in headers], dtype=bool) is_pol60 = np.array([head["filtnam1"].lower() == "pol60" for head in headers], dtype=bool)
@@ -130,7 +129,6 @@ def save_Stokes(
Only returned if return_hdul is True. Only returned if return_hdul is True.
""" """
# Create new WCS object given the modified images # Create new WCS object given the modified images
exp_tot = header_stokes['exptime']
new_wcs = WCS(header_stokes).deepcopy() new_wcs = WCS(header_stokes).deepcopy()
if data_mask.shape != (1, 1): if data_mask.shape != (1, 1):
@@ -140,23 +138,23 @@ def save_Stokes(
new_wcs.wcs.crpix = np.array(new_wcs.wcs.crpix) - vertex[0::-2] new_wcs.wcs.crpix = np.array(new_wcs.wcs.crpix) - vertex[0::-2]
header = new_wcs.to_header() header = new_wcs.to_header()
header["TELESCOP"] = (header_stokes["telescop"] if "TELESCOP" in list(header_stokes.keys()) else "HST", "telescope used to acquire data") header["TELESCOP"] = (header_stokes["TELESCOP"] if "TELESCOP" in list(header_stokes.keys()) else "HST", "telescope used to acquire data")
header["INSTRUME"] = (header_stokes["instrume"] if "INSTRUME" in list(header_stokes.keys()) else "FOC", "identifier for instrument used to acuire data") header["INSTRUME"] = (header_stokes["INSTRUME"] if "INSTRUME" in list(header_stokes.keys()) else "FOC", "identifier for instrument used to acuire data")
header["PHOTPLAM"] = (header_stokes["photplam"], "Pivot Wavelength") header["PHOTPLAM"] = (header_stokes["PHOTPLAM"], "Pivot Wavelength")
header["PHOTFLAM"] = (header_stokes["photflam"], "Inverse Sensitivity in DN/sec/cm**2/Angst") header["PHOTFLAM"] = (header_stokes["PHOTFLAM"], "Inverse Sensitivity in DN/sec/cm**2/Angst")
header["EXPTIME"] = (exp_tot, "Total exposure time in sec") header["EXPTIME"] = (header_stokes["EXPTIME"], "Total exposure time in sec")
header["PROPOSID"] = (header_stokes["proposid"], "PEP proposal identifier for observation") header["PROPOSID"] = (header_stokes["PROPOSID"], "PEP proposal identifier for observation")
header["TARGNAME"] = (header_stokes["targname"], "Target name") header["TARGNAME"] = (header_stokes["TARGNAME"], "Target name")
header["ORIENTAT"] = (np.arccos(new_wcs.wcs.pc[0, 0]) * 180.0 / np.pi, "Angle between North and the y-axis of the image") header["ORIENTAT"] = (header_stokes["ORIENTAT"], "Angle between North and the y-axis of the image")
header["FILENAME"] = (filename, "Original filename") header["FILENAME"] = (filename, "ORIGINAL FILENAME")
header["BKG_TYPE"] = (header_stokes["BKG_TYPE"], "Bkg estimation method used during reduction") header["BKG_TYPE"] = (header_stokes["BKG_TYPE"], "Bkg estimation method used during reduction")
header["BKG_SUB"] = (header_stokes["BKG_SUB"], "Amount of bkg subtracted from images") header["BKG_SUB"] = (header_stokes["BKG_SUB"], "Amount of bkg subtracted from images")
header["SMOOTH"] = (header_stokes["SMOOTH"], "Smoothing method used during reduction") header["SMOOTH"] = (header_stokes["SMOOTH"] if "SMOOTH" in list(header_stokes.keys()) else "None", "Smoothing method used during reduction")
header["SAMPLING"] = (header_stokes["SAMPLING"], "Resampling performed during reduction") header["SAMPLING"] = (header_stokes["SAMPLING"] if "SAMPLING" in list(header_stokes.keys()) else "None", "Resampling performed during reduction")
header["P_INT"] = (header_stokes["P_int"], "Integrated polarization degree") header["P_INT"] = (header_stokes["P_INT"], "Integrated polarization degree")
header["sP_INT"] = (header_stokes["sP_int"], "Integrated polarization degree error") header["sP_INT"] = (header_stokes["sP_INT"], "Integrated polarization degree error")
header["PA_INT"] = (header_stokes["PA_int"], "Integrated polarization angle") header["PA_INT"] = (header_stokes["PA_INT"], "Integrated polarization angle")
header["sPA_INT"] = (header_stokes["sPA_int"], "Integrated polarization angle error") header["sPA_INT"] = (header_stokes["sPA_INT"], "Integrated polarization angle error")
# Crop Data to mask # Crop Data to mask
if data_mask.shape != (1, 1): if data_mask.shape != (1, 1):

49
package/lib/plots.py
View File

@@ -182,9 +182,11 @@ def plot_Stokes(Stokes, savename=None, plots_folder=""):
wcs = WCS(Stokes[0]).deepcopy() wcs = WCS(Stokes[0]).deepcopy()
# Plot figure # Plot figure
plt.rcParams.update({"font.size": 10}) plt.rcParams.update({"font.size": 12})
fig, (axI, axQ, axU) = plt.subplots(ncols=3, figsize=(20, 6), subplot_kw=dict(projection=wcs)) ratiox = max(int(stkI.shape[1]/stkI.shape[0]),1)
fig.subplots_adjust(hspace=0, wspace=0.75, bottom=0.01, top=0.99, left=0.08, right=0.95) ratioy = max(int(stkI.shape[0]/stkI.shape[1]),1)
fig, (axI, axQ, axU) = plt.subplots(ncols=3, figsize=(20*ratiox, 8*ratioy), subplot_kw=dict(projection=wcs))
fig.subplots_adjust(hspace=0, wspace=0.50, bottom=0.01, top=0.99, left=0.07, right=0.97)
fig.suptitle("I, Q, U Stokes parameters") fig.suptitle("I, Q, U Stokes parameters")
imI = axI.imshow(stkI, origin="lower", cmap="inferno") imI = axI.imshow(stkI, origin="lower", cmap="inferno")
@@ -320,9 +322,11 @@ def polarization_map(
print("No pixel with polarization information above requested SNR.") print("No pixel with polarization information above requested SNR.")
# Plot the map # Plot the map
plt.rcParams.update({"font.size": 10}) plt.rcParams.update({"font.size": 12})
plt.rcdefaults() plt.rcdefaults()
fig, ax = plt.subplots(figsize=(10, 10), layout="constrained", subplot_kw=dict(projection=wcs)) ratiox = max(int(stkI.shape[1]/stkI.shape[0]),1)
ratioy = max(int(stkI.shape[0]/stkI.shape[1]),1)
fig, ax = plt.subplots(figsize=(10*ratiox, 10*ratioy), layout="constrained", subplot_kw=dict(projection=wcs))
ax.set(aspect="equal", fc="k") ax.set(aspect="equal", fc="k")
# fig.subplots_adjust(hspace=0, wspace=0, left=0.102, right=1.02) # fig.subplots_adjust(hspace=0, wspace=0, left=0.102, right=1.02)
@@ -439,17 +443,17 @@ def polarization_map(
ax.transAxes, ax.transAxes,
"E", "E",
"N", "N",
length=-0.08, length=-0.05,
fontsize=0.025, fontsize=0.02,
loc=1, loc=1,
aspect_ratio=-1, aspect_ratio=-(stkI.shape[1]/stkI.shape[0]),
sep_y=0.01, sep_y=0.01,
sep_x=0.01, sep_x=0.01,
back_length=0.0, back_length=0.0,
head_length=10.0, head_length=10.0,
head_width=10.0, head_width=10.0,
angle=-Stokes[0].header["orientat"], angle=-Stokes[0].header["orientat"],
text_props={"ec": "k", "fc": "w", "alpha": 1, "lw": -0.2}, text_props={"ec": "k", "fc": "w", "alpha": 1, "lw": 0.4},
arrow_props={"ec": "k", "fc": "w", "alpha": 1, "lw": 1}, arrow_props={"ec": "k", "fc": "w", "alpha": 1, "lw": 1},
) )
@@ -666,7 +670,7 @@ class align_maps(object):
length=-0.08, length=-0.08,
fontsize=0.03, fontsize=0.03,
loc=1, loc=1,
aspect_ratio=-1, aspect_ratio=-(self.map_data.shape[1]/self.map_data.shape[0]),
sep_y=0.01, sep_y=0.01,
sep_x=0.01, sep_x=0.01,
angle=-self.map_header["orientat"], angle=-self.map_header["orientat"],
@@ -724,7 +728,7 @@ class align_maps(object):
length=-0.08, length=-0.08,
fontsize=0.03, fontsize=0.03,
loc=1, loc=1,
aspect_ratio=-1, aspect_ratio=-(self.other_data.shape[1]/self.other_data.shape[0]),
sep_y=0.01, sep_y=0.01,
sep_x=0.01, sep_x=0.01,
angle=-self.other_header["orientat"], angle=-self.other_header["orientat"],
@@ -988,7 +992,7 @@ class overplot_radio(align_maps):
length=-0.08, length=-0.08,
fontsize=0.03, fontsize=0.03,
loc=1, loc=1,
aspect_ratio=-1, aspect_ratio=-(stkI.shape[1]/stkI.shape[0]),
sep_y=0.01, sep_y=0.01,
sep_x=0.01, sep_x=0.01,
angle=-self.Stokes_UV[0].header["orientat"], angle=-self.Stokes_UV[0].header["orientat"],
@@ -1190,7 +1194,7 @@ class overplot_chandra(align_maps):
length=-0.08, length=-0.08,
fontsize=0.03, fontsize=0.03,
loc=1, loc=1,
aspect_ratio=-1, aspect_ratio=-(stkI.shape[1]/stkI.shape[0]),
sep_y=0.01, sep_y=0.01,
sep_x=0.01, sep_x=0.01,
angle=-self.Stokes_UV[0].header["orientat"], angle=-self.Stokes_UV[0].header["orientat"],
@@ -1329,7 +1333,6 @@ class overplot_pol(align_maps):
else: else:
self.scale_vec = scale_vec self.scale_vec = scale_vec
step_vec = 1 step_vec = 1
px_scale = np.abs(self.wcs_UV.wcs.get_cdelt()[0] / self.other_wcs.wcs.get_cdelt()[0])
self.X, self.Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0])) self.X, self.Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0]))
self.U, self.V = pol * np.cos(np.pi / 2.0 + pang * np.pi / 180.0), pol * np.sin(np.pi / 2.0 + pang * np.pi / 180.0) self.U, self.V = pol * np.cos(np.pi / 2.0 + pang * np.pi / 180.0), pol * np.sin(np.pi / 2.0 + pang * np.pi / 180.0)
self.Q = self.ax_overplot.quiver( self.Q = self.ax_overplot.quiver(
@@ -1339,7 +1342,7 @@ class overplot_pol(align_maps):
self.V[::step_vec, ::step_vec], self.V[::step_vec, ::step_vec],
units="xy", units="xy",
angles="uv", angles="uv",
scale=px_scale / self.scale_vec, scale=1. / self.scale_vec,
scale_units="xy", scale_units="xy",
pivot="mid", pivot="mid",
headwidth=0.0, headwidth=0.0,
@@ -1385,7 +1388,7 @@ class overplot_pol(align_maps):
length=-0.08, length=-0.08,
fontsize=0.03, fontsize=0.03,
loc=1, loc=1,
aspect_ratio=-1, aspect_ratio=-(stkI.shape[1]/stkI.shape[0]),
sep_y=0.01, sep_y=0.01,
sep_x=0.01, sep_x=0.01,
angle=-self.Stokes_UV[0].header["orientat"], angle=-self.Stokes_UV[0].header["orientat"],
@@ -1395,7 +1398,7 @@ class overplot_pol(align_maps):
self.ax_overplot.add_artist(north_dir) self.ax_overplot.add_artist(north_dir)
pol_sc = AnchoredSizeBar( pol_sc = AnchoredSizeBar(
self.ax_overplot.transData, self.ax_overplot.transData,
self.scale_vec / px_scale, self.scale_vec,
r"$P$= 100%", r"$P$= 100%",
4, 4,
pad=0.5, pad=0.5,
@@ -1550,7 +1553,7 @@ class align_pol(object):
length=-0.08, length=-0.08,
fontsize=0.025, fontsize=0.025,
loc=1, loc=1,
aspect_ratio=-1, aspect_ratio=-(stkI.shape[1]/stkI.shape[0]),
sep_y=0.01, sep_y=0.01,
sep_x=0.01, sep_x=0.01,
back_length=0.0, back_length=0.0,
@@ -1814,6 +1817,8 @@ class crop_map(object):
# Write cropped map to new HDUList # Write cropped map to new HDUList
self.header_crop = deepcopy(header) self.header_crop = deepcopy(header)
self.header_crop.update(self.wcs_crop.to_header()) self.header_crop.update(self.wcs_crop.to_header())
if self.header_crop["FILENAME"][-4:] != "crop":
self.header_crop["FILENAME"] += "_crop"
self.hdul_crop = fits.HDUList([fits.PrimaryHDU(self.data_crop, self.header_crop)]) self.hdul_crop = fits.HDUList([fits.PrimaryHDU(self.data_crop, self.header_crop)])
self.rect_selector.clear() self.rect_selector.clear()
@@ -1936,6 +1941,8 @@ class crop_Stokes(crop_map):
) )
for dataset in self.hdul_crop: for dataset in self.hdul_crop:
if dataset.header["FILENAME"][-4:] != "crop":
dataset.header["FILENAME"] += "_crop"
dataset.header["P_int"] = (P_diluted, "Integrated polarization degree") dataset.header["P_int"] = (P_diluted, "Integrated polarization degree")
dataset.header["sP_int"] = (np.ceil(P_diluted_err * 1000.0) / 1000.0, "Integrated polarization degree error") dataset.header["sP_int"] = (np.ceil(P_diluted_err * 1000.0) / 1000.0, "Integrated polarization degree error")
dataset.header["PA_int"] = (PA_diluted, "Integrated polarization angle") dataset.header["PA_int"] = (PA_diluted, "Integrated polarization angle")
@@ -2797,10 +2804,10 @@ class pol_map(object):
ax.transAxes, ax.transAxes,
"E", "E",
"N", "N",
length=-0.08, length=-0.05,
fontsize=0.025, fontsize=0.02,
loc=1, loc=1,
aspect_ratio=-1, aspect_ratio=-(self.I.shape[1]/self.I.shape[0]),
sep_y=0.01, sep_y=0.01,
sep_x=0.01, sep_x=0.01,
back_length=0.0, back_length=0.0,