#
# Copyright 2026 by University College Dublin. All rights reserved.
#
# Developed by: Derek O'Callaghan
# University College Dublin
# https://www.ucd.ie/
#
# Builds on:
# Gamma-ray Data Tools - Core Components (https://github.com/USRA-STI/gdt-core)
# Gamma-ray Data Tools - Fermi mission components (https://github.com/USRA-STI/gdt-fermi/)
#
# Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except
# in compliance with the License. You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software distributed under the License
# is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
# implied. See the License for the specific language governing permissions and limitations under the
# License.
#
from collections.abc import Callable
from dataclasses import dataclass, field
from matplotlib.colors import to_rgba
import numpy as np
import pandas as pd
import seaborn as sns
from typing import List, Tuple
from gdt.core.binning.unbinned import bin_by_time
from gdt.core.plot.plot import SAA, SkyHeatmap
from gdt.core.plot.earthplot import EarthPlot
from gdt.core.plot.sky import EquatorialPlot
from gdt.missions.gifts.region import HighParticleFluxRegion
from gdt.missions.gifts.trigdat import Trigdat
from gdt.missions.gifts.tte import GiftsTte
__all__ = ["GiftsEarthPlot", "GiftsEquatorialPlot", "GiftsLightcurve"]
[docs]
class GiftsEarthPlot(EarthPlot):
"""Class for plotting GIFTS orbit and high particle flux regions
Note:
This class requires installation of Cartopy.
Args:
regions: list of GIFTS high particle flux regions to be included in the plot
**kwargs: Options to pass to :class:`~gdt.core.plot.plot.GdtPlot`
"""
def __init__(self, regions: List[HighParticleFluxRegion], **kwargs):
# Default full latitude and longitude ranges will be used for GIFTS plots
super().__init__(**kwargs)
self._region_polygons = []
if regions:
for region in regions:
self.plot_high_particle_flux_region(region)
[docs]
def plot_high_particle_flux_region(self, region: HighParticleFluxRegion, color: str="darkred", alpha: float=0.3):
"""
Include a high particle flux region in the plot
Args:
region: GIFTS high particle flux region for a particular altitude
color: region polygon fill color
alpha: region polygon fill alpha
"""
self._region_polygons.append(SAA(region, self._m.projection, self._m, color=color, alpha=alpha))
[docs]
class GiftsEquatorialPlot(EquatorialPlot):
[docs]
def plot_heatmap(self, heatmap, lon_array, lat_array, **kwargs):
"""Plot a heatmap on the sky.
This is a workaround for the issue in the parent SkyPlot.plot_heatmap(),
where the meshgrid results in a lon_array starting/ending with 0 degrees.
This causes issues where the heatmap generated from a HEALPix location
is populated at both RA=0/360 degrees, resulting in an unexpected output
from the subsequent call to Matplotlib's pcolormesh().
The old GBM Data Tools didn't generate a meshgrid, and just passes the
specified coordinate arrays unchanged.
TODO: confirm whether this is still required with latest gdt-core
Args:
heatmap (np.array): A 2D array of values
ra_array (np.array): The array of longitudinal gridpoints
dec_array (np.array): The array of latitudinal gridpoints
**kwargs: Options to pass to :class:`~gdt.plot.plot.SkyHeatmap`
"""
heatmap = SkyHeatmap(
lon_array,
lat_array,
heatmap,
self.ax,
frame=self._frame,
flipped=self._flipped,
**kwargs
)
return heatmap
[docs]
@dataclass
class GiftsLightcurve:
"""
Temporally bin a set of TTEs, and generate a single lightcurve
plot for all detectors.
"""
trigdat: Trigdat
ttes: List[GiftsTte]
bin_method: Callable = bin_by_time
binning_times: List[float] = field(default_factory=lambda: [1.024 * 2**x for x in [0, -2, -4]])
time_range: Tuple[float, float] = None
energy_range: Tuple[float, float] = None
steps: bool = True
errorbars: bool = False
def __post_init__(self):
# Add some buffer time to end to ensure all steps are plotted
rate_time_range = [self.time_range[0], self.time_range[1] + 5] if self.time_range else None
df = self._tte_lightcurves(time_range=rate_time_range)
df["detector_label"] = df.detector.apply(lambda d: f"{d} {' Trig' if d in self.trigdat.triggered_detectors else ''}")
hue_order = df.binning_time.unique()
num_hue = hue_order.shape[0]
# https://stackoverflow.com/questions/66667334/python-seaborn-alpha-by-hue
hue_alpha = np.linspace(1, 0.15, 3) # May need to restrict number of time bins
palette = [to_rgba(p, a) for p, a in zip(sns.color_palette()[:num_hue], hue_alpha[:num_hue])]
if not self.time_range:
self.time_range = (df.time.min(), df.time.max())
rate_time_range = self.time_range
if self.errorbars:
# will be used for error bars x, matching "steps-post"
df["time_step"] = df.time + (df.binning_time / 2)
plot_df = df[(df.time >= rate_time_range[0]) & (df.time <= rate_time_range[1])]
with sns.axes_style("ticks"):
g = sns.relplot(
data=plot_df,
x="time",
y="rate",
hue="binning_time",
hue_order=hue_order,
size="binning_time",
col="detector_label",
kind="line",
sizes=(0.75, 1.75), # if shorter binning times should be larger, use: 1.75, 0.75),
palette=palette,
height=3,
aspect=1.5,
col_wrap=2,
drawstyle='steps-post' if self.steps else None,
facet_kws=dict(sharex=False),
)
if self.errorbars:
# https://stackoverflow.com/questions/75137582/how-can-i-plot-st-error-bars-with-seaborn-relplot
# add the errorbars
# Iterate through each facet of the facetgrid
for detector_label, ax in g.axes_dict.items():
detector_data = plot_df[plot_df.detector_label==detector_label]
for group, selected in detector_data.groupby("binning_time"):
# plot the errorbar with the correct color for each group
group_index = np.where(np.unique(df.binning_time)==group)[0][0]
ax.errorbar(data=selected,
x="time_step" if steps else "time",
y="rate",
yerr="rate_uncertainty",
marker="none",
color=palette[group_index],
ls="none")
# Plot x limit max will be the specified time range max, to ensure all steps are plotted
(g.set(xlim=self.time_range, ylim=(0, (plot_df.rate + plot_df.rate_uncertainty).max()))
.set_titles('{col_name}')
.set_xlabels("Time (s)")
.set_ylabels("Rate (count / bin s)")
)
g.legend.get_title().set_text("Time bin (s)")
g.fig.subplots_adjust(hspace=0.40, wspace=0.15, top=0.92)
energy_range_title = f"{self.energy_range[0]}-{self.energy_range[1]} keV" if self.energy_range else "all energy bins"
g.fig.suptitle(f"GIFTS {self.trigdat.headers[0]['OBJECT']} Light Curves ({energy_range_title})");
def _tte_lightcurves(self, time_range: Tuple[float, float] = None) -> pd.DataFrame:
"""
Generate a DataFrame containing PHAII rates for the specified TTEs
Args:
time_range: Requested TTE time range
Returns:
Lightcurves in a DataFrame
"""
if not time_range:
# As the underlying GIFTS TTEs can have different time ranges,
# use the narrowest range that's common to all detectors
tte_time_ranges = np.array([tte.data.time_range for tte in self.ttes])
time_range = (tte_time_ranges[:, 0].max(), tte_time_ranges[:, 1].min())
dfs = []
for tte in self.ttes:
for binning_time in self.binning_times:
phaii = tte.to_phaii(self.bin_method, dt=binning_time, time_range=time_range)
lc_data = phaii.to_lightcurve(energy_range=self.energy_range)
dfs.append(
pd.DataFrame(
{
"time": lc_data.lo_edges,
"rate": lc_data.rates,
"rate_uncertainty": lc_data.rate_uncertainty,
"detector": tte.detector,
"binning_time": binning_time,
}
)
)
return pd.concat(dfs)