# This code is inspired by the code developed for ExoRad 2.
# Therefore, we attach here for ExoRad 2 license:
#
# BSD 3-Clause License
#
# Copyright (c) 2020, Lorenzo V. Mugnai, Enzo Pascale, "La Sapienza" Università di Roma
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
# 1. Redistributions of source code must retain the above copyright notice, this
# list of conditions and the following disclaimer.
#
# 2. Redistributions in binary form must reproduce the above copyright notice,
# this list of conditions and the following disclaimer in the documentation
# and/or other materials provided with the distribution.
#
# 3. Neither the names of the copyright holders nor the names of its
# contributors may be used to endorse or promote products derived from
# this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
import os
import astropy.units as u
import h5py
import matplotlib as mpl
import matplotlib.patches as mpatches
import matplotlib.pyplot as plt
import numpy as np
from astropy.io.misc.hdf5 import read_table_hdf5
from astropy.table import Table
from photutils.aperture import EllipticalAperture, RectangularAperture
import exosim.log as log
from exosim.output.hdf5.utils import load_signal
from exosim.utils.ascii_arts import observatory
from .utils import _create_ordered_cmap, prepare_channels_list
plt.rcParams.update({"font.size": 16})
[docs]
class RadiometricPlotter(log.Logger):
"""
Radiometric plotter.
This class handles the methods to plot the radiometric table produced by `exosim`.
Attributes
-----------
input: str or :class:`astropy.table.QTable`
input data
input_table: :class:`astropy.table.QTable`
input radiometric table
fig: :class:`matplotlib.figure.Figure`
produced figure
Examples
----------
The following example, given the `test_file.h5` preoduced by Exosim,
plots the radiometric table and stores the figure as `radiometric.png`.
>>> from exosim.plots import RadiometricPlotter
>>> radiometricPlotter = RadiometricPlotter(
... input="./test_file.h5"
... )
>>> radiometricPlotter.plot_table()
>>> radiometricPlotter.save_fig("radiometric.png")
"""
def __init__(self, input: str | Table, table_path: str = "radiometric") -> None:
"""
Parameters
----------
input: str or :class:`astropy.table.QTable`
input data
"""
self.set_log_name()
self.graphics(observatory)
self.announce("started")
if isinstance(input, str):
self.input_table = self.load_table(input, table_path)
else:
# Accept any QTable-like input; specific column checks will be deferred to plotting methods
if hasattr(input, "colnames") and len(input) == 0:
raise ValueError("Empty table")
self.input_table = input
[docs]
def load_table(self, input_file: str, table_path: str = "radiometric") -> Table:
"""
It loads the radiometric table from the input file:
Parameters
----------
input_file: str
input file name
Returns
-------
:class:`astropy.table.QTable`
loaded radiometric table
"""
with h5py.File(input_file, "r") as f:
tab = read_table_hdf5(f[table_path], path="table")
self.debug("radiometric table loaded")
return tab
[docs]
def plot_bands(
self,
ax: plt.Axes,
scale: str = "log",
channel_edges: bool = True,
add_legend: bool = True,
) -> plt.Axes:
"""
It plots the channels bands behind the indicated axes.
Parameters
-----------
ax: :class:`matplotlib.axes.Axes`
axes where to plot the bands
scale: str
x axes scale. Default is `log`.
channel_edges: bool
if ``True`` the x axes ticks are placed at the channel edges. Default is ``True``.
Returns
--------
:class:`matplotlib.axes.Axes`
axes with channel bands added
"""
channels, norm = prepare_channels_list(self.input_table)
cmap = _create_ordered_cmap("Pastel1", roll=-2, delete=-3)
tick_list, patches = [], []
for k, channel_name in enumerate(channels):
wl_min = min(
self.input_table["left_bin_edge"][
np.where(self.input_table["ch_name"] == channel_name)
]
)
if hasattr(wl_min, "unit"):
wl_min = wl_min.value
wl_max = max(
self.input_table["right_bin_edge"][
np.where(self.input_table["ch_name"] == channel_name)
]
)
if hasattr(wl_max, "unit"):
wl_max = wl_max.value
ax.axvspan(
wl_min,
wl_max,
alpha=0.3,
zorder=0,
color=cmap(
norm(k),
),
)
ax.axvspan(
wl_min,
wl_max,
alpha=0.3,
zorder=0,
color=cmap(
norm(k),
),
)
# wl_maxs += [wl_max]
tick_list.append(wl_min)
tick_list.append(wl_max)
patches += [
mpatches.Patch(color=cmap(norm(k)), alpha=0.3, label=channel_name)
]
# tick_list.append(max(wl_maxs))
ax.set_xscale(scale)
if channel_edges:
# Remove duplicate ticks and ticks that are too close
tick_list_sorted = sorted(set(tick_list))
# Filter out ticks that are too close (within 5% relative distance)
filtered_ticks = [tick_list_sorted[0]]
for tick in tick_list_sorted[1:]:
if (tick - filtered_ticks[-1]) / filtered_ticks[-1] > 0.05:
filtered_ticks.append(tick)
ax.set_xticks(filtered_ticks)
# Format with 3 significant figures
formatter = mpl.ticker.FuncFormatter(lambda x, p: f"{x:.3g}")
ax.get_xaxis().set_major_formatter(formatter)
if add_legend:
# Store patches for later use in combined legend
self._band_patches = patches
return ax
[docs]
def plot_noise(
self,
ax: plt.Axes,
scale: str = "log",
channel_edges: bool = True,
contribs: bool = False,
ch_lengend: bool = True,
ylim: tuple[float, float] | None = None,
) -> plt.Axes:
"""
It plots the noise components found in the input table in the indicated axes.
Parameters
-----------
ax: :class:`matplotlib.axes.Axes`
axes where to plot the noises
scale: str
x axes scale. Default is `log`.
channel_edges: bool
if ``True`` the x axes ticks are placed at the channel edges. Default is ``True``.
contribs: bool
if ``True`` all the contributions are plotted. Default is ``False``.
ch_lengend: bool
if ``True`` add a legend for the channels color. Default is ``True``.
ylim: tuple[float, float] | None
ylim for :class:`matplotlib.axes.Axes`.
Returns
--------
:class:`matplotlib.axes.Axes`
axes with noises plotted
"""
noise_keys = [
x for x in self.input_table.colnames if "noise" in x or "custom" in x
]
if not contribs:
noise_keys = [k for k in noise_keys if "photon_noise" not in k]
noise_keys += ["source_photon_noise", "foreground_photon_noise"]
self.debug(f"noise keys : {noise_keys}")
for _k, n in enumerate(noise_keys):
if n == "total_noise":
ax.plot(
self.input_table["wavelength"],
self.input_table[n],
zorder=9,
lw=1,
c="k",
marker=".",
markersize=5,
label="total_noise",
alpha=0.8,
) # , c='None')
# ignore variance columns
elif "variance" in n:
continue
else:
if self.input_table[n].unit == u.hr**0.5:
noise = self.input_table[n]
elif self.input_table[n].unit == u.ct / u.s:
self.debug(f"{n} rescaled by starSignal_inAperture")
noise = (
self.input_table[n]
/ self.input_table["source_signal_in_aperture"]
/ (u.hr.to(u.s)) ** 0.5
)
else:
self.error(f"{n} unit not valid: {self.input_table[n].unit}")
# ax.scatter(self.input_table['wavelength'], noise, label=n, zorder=10, s=5, color=palette[k])
ax.plot(
self.input_table["wavelength"],
noise,
zorder=9,
lw=1,
alpha=0.5,
marker=".",
label=n,
) # color=palette[k]) # c='None')
# ax.grid(zorder=0, which='both')
locmaj = mpl.ticker.LogLocator(base=10, numticks=12)
ax.yaxis.set_major_locator(locmaj)
locmin = mpl.ticker.LogLocator(
base=10.0, subs=(0.2, 0.4, 0.6, 0.8), numticks=12
)
ax.yaxis.set_minor_locator(locmin)
ax.yaxis.set_minor_formatter(mpl.ticker.NullFormatter())
ax.grid(axis="y", which="minor", alpha=0.3)
ax.grid(axis="y", which="major", alpha=0.5)
# ax.legend(bbox_to_anchor=(1, 1))
ax.set_title("Noise Budget")
ax.set_xlabel(r"wavelength [$\mu m$]")
ax.set_ylabel(r"relative noise [$\sqrt{{hr}}$]")
ax.set_yscale(scale)
if ylim:
ax.set_ylim(ylim)
# ax.set_xscale('log')
ax = self.plot_bands(ax, scale, channel_edges, add_legend=ch_lengend)
# Smart legend sizing based on number of items
handles, labels = ax.get_legend_handles_labels()
n_items = len(labels)
# More aggressive column scaling: 4-8 columns
ncol = min(max(4, int(np.ceil(n_items / 2))), 8)
# Larger font scaling: 8-12pt (was 6-10pt)
fontsize = max(8, 12 - n_items // 5)
# Calculate vertical offset based on expected legend height
legend_rows = int(np.ceil(n_items / ncol))
y_offset = -0.15 - (legend_rows - 1) * 0.10 # More space
# Create main legend for data at the bottom
data_legend = ax.legend(
handles=handles,
labels=labels,
prop={"size": fontsize},
loc="upper left",
ncol=ncol,
bbox_to_anchor=(0.0, y_offset),
labelspacing=0.6, # More spacing
handlelength=1.0, # Longer handles
columnspacing=1.0, # More column spacing
framealpha=0.9,
)
# Add channel band legend on the right side if requested
if ch_lengend and hasattr(self, "_band_patches"):
ax.add_artist(data_legend) # Keep the data legend
ax.legend(
handles=self._band_patches,
title="Channels",
loc="upper left",
bbox_to_anchor=(1.02, 1),
fontsize=11, # Larger font (was 9)
title_fontsize=12, # Larger title (was 10)
framealpha=0.9,
)
return ax
[docs]
def plot_signal(
self,
ax: plt.Axes,
ylim: tuple[float, float] | None = None,
scale: str = "log",
channel_edges: bool = True,
contribs: bool = False,
ch_lengend: bool = True,
) -> tuple[plt.Figure, tuple[plt.Axes, plt.Axes]]:
"""
It plots the signal components found in the input table in the indicated axes.
Parameters
-----------
ylim: float or (float, float)
ylim for :class:`matplotlib.axes.Axes`.
ax: :class:`matplotlib.axes.Axes`
axes where to plot the signals
scale: str
x axes scale. Default is `log`.
channel_edges: bool
if ``True`` the x axes ticks are placed at the channel edges. Default is ``True``.
contribs: bool
if ``True`` all the contributions are plotted. Default is ``False``.
ch_lengend: bool
if ``True`` add a legend for the channels color. Default is ``True``.
Returns
--------
:class:`matplotlib.axes.Axes`
axes with signals plotted
"""
keys = ["source_signal_in_aperture", "foreground_signal_in_aperture"]
if contribs:
keys = [
x
for x in self.input_table.colnames
if "signal_in_aperture" in x and "noise" not in x
]
self.debug(f"signal keys : {keys}")
for _k, s in enumerate(keys):
ax.plot(
self.input_table["wavelength"],
self.input_table[s],
zorder=9,
lw=1,
alpha=0.5,
marker=".",
label=s,
)
if ylim:
ax.set_ylim(ylim)
# elif ax.get_ylim()[0]<:
# ax.set_ylim(1e-3)
# ax.grid(zorder=0, which='both')
locmaj = mpl.ticker.LogLocator(base=10, numticks=12)
ax.yaxis.set_major_locator(locmaj)
locmin = mpl.ticker.LogLocator(
base=10.0, subs=(0.2, 0.4, 0.6, 0.8), numticks=12
)
ax.yaxis.set_minor_locator(locmin)
ax.yaxis.set_minor_formatter(mpl.ticker.NullFormatter())
ax.grid(axis="y", which="minor", alpha=0.3)
ax.grid(axis="y", which="major", alpha=0.5)
# ax.legend(bbox_to_anchor=(1, 1))
ax.set_title("Signals")
ax.set_xlabel(r"wavelength [$\mu m$]")
ax.set_ylabel("$ct/s$")
ax.set_yscale(scale)
ax = self.plot_bands(ax, scale, channel_edges, add_legend=ch_lengend)
# Smart legend sizing based on number of items
handles, labels = ax.get_legend_handles_labels()
n_items = len(labels)
# More aggressive column scaling: 4-8 columns
ncol = min(max(4, int(np.ceil(n_items / 2))), 8)
# Larger font scaling: 8-12pt (was 6-10pt)
fontsize = max(8, 12 - n_items // 5)
# Calculate vertical offset based on expected legend height
legend_rows = int(np.ceil(n_items / ncol))
y_offset = -0.15 - (legend_rows - 1) * 0.10 # More space
# Create main legend for data at the bottom
data_legend = ax.legend(
handles=handles,
labels=labels,
prop={"size": fontsize},
loc="upper left",
ncol=ncol,
bbox_to_anchor=(0.0, y_offset),
labelspacing=0.6, # More spacing
handlelength=1.0, # Longer handles
columnspacing=1.0, # More column spacing
framealpha=0.9,
)
# Add channel band legend on the right side if requested
if ch_lengend and hasattr(self, "_band_patches"):
ax.add_artist(data_legend) # Keep the data legend
ax.legend(
handles=self._band_patches,
title="Channels",
loc="upper left",
bbox_to_anchor=(1.02, 1),
fontsize=11, # Larger font (was 9)
title_fontsize=12, # Larger title (was 10)
framealpha=0.9,
)
return ax
[docs]
def plot_table(
self,
title: str = "Radiometric Table",
scale: str = "log",
channel_edges: bool = True,
contribs: bool = False,
signal_ylim: tuple[float, float] | None = None,
noise_ylim: tuple[float, float] | None = None,
) -> tuple[plt.Figure, tuple[plt.Axes, plt.Axes]]:
"""
It produces a figure with signal and noise for the input table.
Parameters
----------
scale: str
x axes scale. Default is `log`.
channel_edges: bool
if ``True`` the x axes ticks are placed at the channel edges. Default is ``True``.
contribs: bool
if ``True`` all the contributions are plotted. Default is ``False``.
signal_ylim: tuple[float, float] | None
ylim for signal :class:`matplotlib.axes.Axes`.
noise_ylim: tuple[float, float] | None
ylim for noise :class:`matplotlib.axes.Axes`.
Returns
--------
:class:`matplotlib.figure.Figure`
plotted figure
(:class:`matplotlib.axes.Axes`, :class:`matplotlib.axes.Axes`)
tuple of axis. First axes is for signal, second is for noise.
"""
self.info("plotting radiometric table")
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 12))
fig.suptitle(title)
ax1 = self.plot_signal(
ax1,
scale=scale,
channel_edges=channel_edges,
contribs=contribs,
ylim=signal_ylim,
ch_lengend=True, # Show channel legend on first plot
)
ax2 = self.plot_noise(
ax2,
scale=scale,
channel_edges=channel_edges,
contribs=contribs,
ch_lengend=False, # Hide channel legend on second plot (already shown above)
ylim=noise_ylim,
)
# Calculate needed space based on number of legend items
n_signal = len(
[
k
for k in self.input_table.colnames
if "signal_in_aperture" in k and "noise" not in k
]
)
n_noise = len(
[k for k in self.input_table.colnames if "noise" in k or "custom" in k]
)
max_items = max(n_signal if contribs else 2, n_noise)
# Estimate rows needed (assuming ~6 columns average)
legend_rows = int(np.ceil(max_items / 6))
bottom_margin = 0.15 + legend_rows * 0.05 # More space (was 0.12 + 0.04)
# Use tight_layout with padding, then fine-tune
plt.tight_layout(pad=2.0)
# Adjust margins: more space between plots (hspace) and for legends (bottom)
plt.subplots_adjust(
top=0.94,
bottom=bottom_margin,
hspace=0.65, # More space between subplots
right=0.85, # More space for channel legend on right
left=0.08,
)
self.fig = fig
return fig, (ax1, ax2)
[docs]
def plot_efficiency(
self,
scale: str = "log",
channel_edges: bool = False,
ch_lengend: bool = True,
) -> tuple[plt.Figure, tuple[plt.Axes, plt.Axes]]:
"""
It produces a figure with efficiencies for the input table.
Parameters
----------
scale: str
x axes scale. Default is `log`.
channel_edges: bool
if ``True`` the x axes ticks are placed at the channel edges. Default is ``True``.
ch_lengend: bool
if ``True`` add a legend for the channels color. Default is ``True``.
Returns
--------
:class:`matplotlib.figure.Figure`
plotted figure
(:class:`matplotlib.axes.Axes`, :class:`matplotlib.axes.Axes`)
tuple of axis. First axes is for signal, second is for noise.
"""
self.info("plotting efficiency table")
# fig.suptitle(self.input_table.meta['name'])
channels = set(self.input_table["ch_name"])
channels = list(channels)
channels.sort()
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 10))
with h5py.File(self.input, "r") as f:
for ch in channels:
eff_path = f"channels/{ch}/efficiency"
eff = load_signal(f[eff_path])
ax1.plot(eff.spectral, eff.data[0, 0], label=ch)
resp_path = f"channels/{ch}/responsivity"
resp = load_signal(f[resp_path])
ax2.plot(resp.spectral, resp.data[0, 0], label=ch)
# locmaj = matplotlib.ticker.LogLocator(base=10, numticks=12)
# ax.yaxis.set_major_locator(locmaj)
# locmin = matplotlib.ticker.LogLocator(base=10.0,
# subs=(0.2, 0.4, 0.6, 0.8),
# numticks=12)
# ax.yaxis.set_minor_locator(locmin)
# ax.yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())
# ax.grid(axis='y', which='minor', alpha=0.3)
# ax.grid(axis='y', which='major', alpha=0.5)
# # ax.legend(bbox_to_anchor=(1, 1))
ax1.set_title("Efficiency")
ax1.set_xlabel(r"wavelength [$\mu m$]")
# ax.set_ylabel('$ct/s$')
# ax1.set_yscale(scale)
ax2.set_title("Responsivity")
ax2.set_ylabel(rf"${resp.data_units}$")
ax2.set_xlabel(r"wavelength [$\mu m$]")
ax2.set_yscale("log")
ax1 = self.plot_bands(ax1, scale, channel_edges, add_legend=ch_lengend)
ax2 = self.plot_bands(ax2, scale, channel_edges, add_legend=ch_lengend)
ax2.legend(
prop={"size": 12},
loc="upper left",
ncol=7,
bbox_to_anchor=(0.05, -0.25),
labelspacing=1.2,
handlelength=1,
)
plt.tight_layout()
plt.subplots_adjust(top=0.92, bottom=0.12, hspace=0.7, right=0.8, left=0.1)
self.fig = fig
return fig, (ax1, ax2)
[docs]
def plot_apertures(self, channel_folder=None) -> plt.Figure:
"""
It produces a figure with apertures superimposed to the focal plane.
Returns
--------
:class:`matplotlib.figure.Figure`
plotted figure
"""
self.info("plotting apertures")
def _prepare_figure(
channel_folder=None,
) -> tuple[plt.Figure, plt.GridSpec, list[str]]:
if channel_folder is None:
self.debug("loading channel folder from default input file")
with h5py.File(self.input, "r") as f:
channel_folder = f["channels"]
ch_list = list(channel_folder.keys())
ch_list.sort()
# Pre-load the shape information since we'll lose access to the file
widths = []
for ch in ch_list:
size_x = channel_folder[ch]["focal_plane"]["spectral"].shape[0]
size_y = channel_folder[ch]["focal_plane"]["spatial"].shape[0]
widths += [int(np.ceil(size_x / size_y)), 0.1]
else:
ch_list = list(channel_folder.keys())
ch_list.sort()
widths = []
for ch in ch_list:
size_x = channel_folder[ch]["focal_plane"]["spectral"].shape[0]
size_y = channel_folder[ch]["focal_plane"]["spatial"].shape[0]
widths += [int(np.ceil(size_x / size_y)), 0.1]
heights = [1]
scale = np.ceil(len(widths) / len(heights))
size_y_fig = 10
size_x_fig = size_y_fig * scale
fig = plt.figure(
constrained_layout=True,
dpi=150,
figsize=(size_x_fig, size_y_fig),
)
spec = fig.add_gridspec(
ncols=len(widths),
nrows=len(heights),
width_ratios=widths,
height_ratios=heights,
wspace=0.1,
hspace=0.1,
)
return fig, spec, ch_list
def _load_apertures(ch):
center_spectral = self.input_table["spectral_center"][
self.input_table["ch_name"] == ch
]
spectral_size = self.input_table["spectral_size"][
self.input_table["ch_name"] == ch
]
center_spatial = self.input_table["spatial_center"][
self.input_table["ch_name"] == ch
]
spatial_size = self.input_table["spatial_size"][
self.input_table["ch_name"] == ch
]
shape = self.input_table["aperture_shape"][
self.input_table["ch_name"] == ch
]
aperture_shapes = {
"rectangular": RectangularAperture,
"elliptical": EllipticalAperture,
}
aper = []
for i in range(center_spectral.size):
aper += [
aperture_shapes[shape[i]](
(center_spectral[i], center_spatial[i]),
spectral_size[i],
spatial_size[i],
)
]
return aper
fig, spec, ch_list = _prepare_figure(channel_folder)
i = 0
for ch in ch_list:
if channel_folder is not None:
file_path = channel_folder[ch]
focal_plane = load_signal(file_path["focal_plane"])
foreground = load_signal(file_path["frg_focal_plane"])
osf = focal_plane.metadata["oversampling"]
final = focal_plane.data[0] + foreground.data[0]
ax0 = fig.add_subplot(spec[0, i])
im = ax0.imshow(
final[osf // 2 :: osf, osf // 2 :: osf],
interpolation="none",
)
else:
with h5py.File(self.input, "r") as f:
file_path = os.path.join("channels", ch)
focal_plane = load_signal(f[os.path.join(file_path, "focal_plane")])
foreground = load_signal(
f[os.path.join(file_path, "frg_focal_plane")]
)
osf = focal_plane.metadata["oversampling"]
final = focal_plane.data[0] + foreground.data[0]
ax0 = fig.add_subplot(spec[0, i])
im = ax0.imshow(
final[osf // 2 :: osf, osf // 2 :: osf],
interpolation="none",
)
ax0.set_title(ch)
apertures = _load_apertures(ch)
for aperture in apertures:
aperture.plot(color="r", lw=2)
i += 1
plt.colorbar(im, ax=ax0, cax=fig.add_subplot(spec[i]))
i += 1
self.fig = fig
return fig
[docs]
def save_fig(self, name: str) -> None:
"""
It saves the produced figure.
Parameters
--------
name: str
figure name
"""
dir_name = os.path.dirname(os.path.abspath(name))
if not os.path.exists(dir_name):
os.makedirs(dir_name)
try:
# Use bbox_inches='tight' to ensure legends are not cut off
self.fig.savefig(f"{name}", bbox_inches="tight")
self.info(f"plot saved in {name}")
except AttributeError:
self.error(
"the indicated figure is not available. Check if you have produced it."
)
