# Licensed under a 3-clause BSD style license - see LICENSE.rst
import logging
import numpy as np
from astropy.units import Quantity
from astropy.io import fits
__all__ = ["Energy", "EnergyBounds"]
log = logging.getLogger(__name__)
[docs]class Energy(Quantity):
"""Energy quantity scalar or array.
This is a `~astropy.units.Quantity` sub-class that adds convenience methods
to handle common tasks for energy bin center arrays, like FITS I/O or generating
equal-log-spaced grids of energies.
See :ref:`energy_handling_gammapy` for further information.
Parameters
----------
energy : `~numpy.array`, scalar, `~astropy.units.Quantity`
Energy
unit : `~astropy.units.UnitBase`, str, optional
The unit of the value specified for the energy. This may be
any string that `~astropy.units.Unit` understands, but it is
better to give an actual unit object.
dtype : `~numpy.dtype`, optional
See `~astropy.units.Quantity`.
copy : bool, optional
See `~astropy.units.Quantity`.
"""
def __new__(cls, energy, unit=None, dtype=None, copy=True):
# Techniques to subclass Quantity taken from astropy.coordinates.Angle
# see: http://docs.scipy.org/doc/numpy/user/basics.subclassing.html
if isinstance(energy, str):
val, unit = energy.split()
energy = float(val)
# This is a pylint error false positive
# See https://github.com/PyCQA/pylint/issues/2335#issuecomment-415055075
self = super().__new__(
cls,
energy,
unit, # pylint:disable=redundant-keyword-arg
dtype=dtype,
copy=copy,
)
if not self.unit.is_equivalent("eV"):
raise ValueError(
"Given unit {} is not an" " energy".format(self.unit.to_string())
)
return self
def __array_finalize__(self, obj):
super().__array_finalize__(obj)
def __quantity_subclass__(self, unit):
if unit.is_equivalent("eV"):
return Energy, True
else:
return super().__quantity_subclass__(unit)[0], False
@property
def nbins(self):
"""The number of bins."""
return self.size
@property
def range(self):
"""The covered energy range (tuple)."""
return self[0 : self.size : self.size - 1]
[docs] @classmethod
def equal_log_spacing(cls, emin, emax, nbins, unit=None, per_decade=False):
"""Create Energy with equal log-spacing (`~gammapy.utils.energy.Energy`).
if no unit is given, it will be taken from emax
Parameters
----------
emin : `~astropy.units.Quantity`, float
Lowest energy bin
emax : `~astropy.units.Quantity`, float
Highest energy bin
nbins : int
Number of bins
unit : `~astropy.units.UnitBase`, str
Energy unit
per_decade : bool
Whether nbins is per decade.
"""
if unit is not None:
emin = Energy(emin, unit)
emax = Energy(emax, unit)
else:
emin = Energy(emin)
emax = Energy(emax)
unit = emax.unit
emin = emin.to(unit)
x_min, x_max = np.log10([emin.value, emax.value])
if per_decade:
nbins = (x_max - x_min) * nbins
energy = np.logspace(x_min, x_max, nbins)
return cls(energy, unit, copy=False)
[docs] @classmethod
def from_fits(cls, hdu, unit=None):
"""Read ENERGIES fits extension (`~gammapy.utils.energy.Energy`).
Parameters
----------
hdu: `~astropy.io.fits.BinTableHDU`
``ENERGIES`` extensions.
unit : `~astropy.units.UnitBase`, str, None
Energy unit
"""
header = hdu.header
fitsunit = header.get("TUNIT1")
if fitsunit is None:
if unit is not None:
log.warning(
"No unit found in the FITS header." " Setting it to {}".format(unit)
)
fitsunit = unit
else:
raise ValueError(
"No unit found in the FITS header." " Please specifiy a unit"
)
energy = cls(hdu.data["Energy"], fitsunit)
return energy.to(unit)
[docs] def to_fits(self):
"""Write ENERGIES fits extension
Returns
-------
hdu: `~astropy.io.fits.BinTableHDU`
ENERGIES fits extension
"""
col1 = fits.Column(name="Energy", format="D", array=self.value)
cols = fits.ColDefs([col1])
hdu = fits.BinTableHDU.from_columns(cols)
hdu.name = "ENERGIES"
hdu.header["TUNIT1"] = "{}".format(self.unit.to_string("fits"))
return hdu
[docs]class EnergyBounds(Energy):
"""EnergyBounds array.
This is a `~gammapy.utils.energy.Energy` sub-class that adds convenience
methods to handle common tasks for energy bin edges arrays, like FITS I/O or
generating arrays of bin centers.
See :ref:`energy_handling_gammapy` for further information.
Parameters
----------
energy : `~numpy.array`, scalar, `~astropy.units.Quantity`
EnergyBounds
unit : `~astropy.units.UnitBase`, str
The unit of the values specified for the energy. This may be any
string that `~astropy.units.Unit` understands, but it is better to
give an actual unit object.
"""
@property
def nbins(self):
"""The number of bins."""
return self.size - 1
@property
def log_centers(self):
"""Log centers of the energy bounds."""
center = np.sqrt(self[:-1] * self[1:])
return center.view(Energy)
@property
def upper_bounds(self):
"""Upper energy bin edges."""
return self[1:]
@property
def lower_bounds(self):
"""Lower energy bin edges."""
return self[:-1]
@property
def boundaries(self):
"""Energy range."""
return self[[0, -1]]
@property
def bands(self):
"""Width of the energy bins."""
upper = self.upper_bounds
lower = self.lower_bounds
return upper - lower
[docs] @classmethod
def from_lower_and_upper_bounds(cls, lower, upper, unit=None):
"""EnergyBounds from lower and upper bounds (`~gammapy.utils.energy.EnergyBounds`).
If no unit is given, it will be taken from upper.
Parameters
----------
lower,upper : `~astropy.units.Quantity`, float
Lowest and highest energy bin
unit : `~astropy.units.UnitBase`, str, None
Energy units
"""
# np.append renders Quantities dimensionless
# http://docs.astropy.org/en/latest/known_issues.html#quantity-issues
if unit is None:
unit = upper.unit
lower = cls(lower, unit)
upper = cls(upper, unit)
energy = np.append(lower.value, upper.value[-1])
return cls(energy, unit)
[docs] @classmethod
def equal_log_spacing(cls, emin, emax, nbins, unit=None):
"""EnergyBounds with equal log-spacing (`~gammapy.utils.energy.EnergyBounds`).
If no unit is given, it will be taken from emax.
Parameters
----------
emin : `~astropy.units.Quantity`, float
Lowest energy bin
emax : `~astropy.units.Quantity`, float
Highest energy bin
bins : int
Number of bins
unit : `~astropy.units.UnitBase`, str, None
Energy unit
"""
return super().equal_log_spacing(emin, emax, nbins + 1, unit)
[docs] @classmethod
def from_ebounds(cls, hdu):
"""Read EBOUNDS fits extension (`~gammapy.utils.energy.EnergyBounds`).
Parameters
----------
hdu: `~astropy.io.fits.BinTableHDU`
``EBOUNDS`` extensions.
"""
if hdu.name != "EBOUNDS":
log.warning(
"This does not seem like an EBOUNDS extension. " "Are you sure?"
)
header = hdu.header
unit = header.get("TUNIT2")
low = hdu.data["E_MIN"]
high = hdu.data["E_MAX"]
return cls.from_lower_and_upper_bounds(low, high, unit)
[docs] @classmethod
def from_rmf_matrix(cls, hdu):
"""Read MATRIX fits extension (`~gammapy.utils.energy.EnergyBounds`).
Parameters
----------
hdu: `~astropy.io.fits.BinTableHDU`
``MATRIX`` extensions.
"""
if hdu.name != "MATRIX":
log.warning("This does not seem like a MATRIX extension. " "Are you sure?")
header = hdu.header
unit = header.get("TUNIT1")
low = hdu.data["ENERG_LO"]
high = hdu.data["ENERG_HI"]
return cls.from_lower_and_upper_bounds(low, high, unit)
[docs] def find_energy_bin(self, energy):
"""Find the bins that contain the specified energy values.
Parameters
----------
energy : `~gammapy.utils.energy.Energy`
Array of energies to search for.
Returns
-------
bin_index : `~numpy.ndarray`
Indices of the energy bins containing the specified energies.
"""
# check that the specified energy is within the boundaries
if not self.contains(energy).all():
ss_error = "Specified energy {}".format(energy)
ss_error += " is outside the boundaries {}".format(self.boundaries)
raise ValueError(ss_error)
bin_index = np.searchsorted(self.upper_bounds, energy)
return bin_index
[docs] def contains(self, energy):
"""Check of energy is contained in boundaries.
Parameters
----------
energy : `~gammapy.utils.energy.Energy`
Array of energies to test
"""
return (energy > self[0]) & (energy < self[-1])
[docs] def to_dict(self):
"""Construct dict representing an energy range."""
if len(self) != 2:
raise ValueError(
"This is not an energy range. Nbins: {}".format(self.nbins)
)
d = dict(min=self[0].value, max=self[1].value, unit="{}".format(self.unit))
return d
[docs] @classmethod
def from_dict(cls, d):
"""Read dict representing an energy range."""
return cls((d["min"], d["max"]), d["unit"])