# Licensed under a 3-clause BSD style license - see LICENSE.rst
import numpy as np
from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.io import fits
from astropy.nddata import NoOverlapError
from regions import PointSkyRegion
from gammapy.maps import HpxNDMap, Map, MapAxis, RegionNDMap
from gammapy.maps.hpx.io import HPX_FITS_CONVENTIONS, HpxConv
from gammapy.utils.scripts import make_path
from gammapy.utils.time import time_ref_from_dict
from . import LightCurveTemplateTemporalModel, Models, SkyModel, TemplateSpatialModel
__all__ = ["read_hermes_cube"]
def _template_model_from_cta_sdc(filename):
"""To create a `LightCurveTemplateTemporalModel` from the energy-dependent temporal model files of the cta-sdc1.
This format is subject to change.
"""
filename = str(make_path(filename))
with fits.open(filename) as hdul:
frame = hdul[0].header.get("frame", "icrs")
position = SkyCoord(
hdul[0].header["LONG"] * u.deg, hdul[0].header["LAT"] * u.deg, frame=frame
)
energy_hdu = hdul["ENERGIES"]
energy_axis = MapAxis.from_nodes(
nodes=energy_hdu.data,
unit=energy_hdu.header["TUNIT1"],
name="energy",
interp="log",
)
time_hdu = hdul["TIMES"]
time_header = time_hdu.header
time_header.setdefault("MJDREFF", 0.5)
time_header.setdefault("MJDREFI", 55555)
time_header.setdefault("scale", "tt")
time_min = time_hdu.data["Initial Time"]
time_max = time_hdu.data["Final Time"]
edges = np.append(time_min, time_max[-1]) * u.Unit(time_header["TUNIT1"])
data = hdul["SPECTRA"]
time_ref = time_ref_from_dict(time_header)
time_axis = MapAxis.from_edges(edges=edges, name="time", interp="log")
reg_map = RegionNDMap.create(
region=PointSkyRegion(center=position),
axes=[energy_axis, time_axis],
data=np.array(list(data.data) * u.Unit(data.header["UNITS"])),
)
return LightCurveTemplateTemporalModel(reg_map, t_ref=time_ref, filename=filename)
[docs]def read_hermes_cube(filename):
"""Read 3d templates produced with hermes."""
# add hermes conventions to the list used by gammapy
hermes_conv = HpxConv(
convname="hermes-template",
colstring="TFLOAT",
hduname="xtension",
frame="COORDTYPE",
quantity_type="differential",
)
HPX_FITS_CONVENTIONS["hermes-template"] = hermes_conv
maps = []
energy_nodes = []
with fits.open(filename) as hdulist:
# cannot read directly in 3d with Map.read because BANDS HDU is missing
# https://gamma-astro-data-formats.readthedocs.io/en/v0.2/skymaps/index.html#bands-hdu
# so we have to loop over hdus and create the energy axis
for hdu in hdulist[1:]:
template = HpxNDMap.from_hdu(hdu, format="hermes-template")
# fix missing/incompatible infos
template._unit = u.Unit(hdu.header["TUNIT1"]) # .from_hdu expect "BUNIT"
if template.geom.frame == "G":
template._geom._frame = "galactic"
maps.append(template)
energy_nodes.append(hdu.header["ENERGY"]) # SI unit (see header comment)
# create energy axis and set unit
energy_nodes *= u.Joule
energy_nodes = energy_nodes.to("GeV")
axis = MapAxis(
energy_nodes, interp="log", name="energy_true", node_type="center", unit="GeV"
)
return Map.from_stack(maps, axis=axis)
def cutout_template_models(models, cutout_kwargs, datasets_names=None):
"""apply cutout to template models"""
models_cut = Models()
for m in models:
if isinstance(m.spatial_model, TemplateSpatialModel):
try:
map_ = m.spatial_model.map.cutout(**cutout_kwargs)
except (NoOverlapError, ValueError):
continue
template_cut = TemplateSpatialModel(
map_,
normalize=m.spatial_model.normalize,
)
model_cut = SkyModel(
spatial_model=template_cut,
spectral_model=m.spectral_model,
datasets_names=datasets_names,
)
models_cut.append(model_cut)
else:
models_cut.append(m)
return models_cut