Remove Energy and EnergyBounds classes

docs/utils/index.rst

@@ -109,67 +109,10 @@ TODO: discuss when to use `~astropy.time.TimeDelta` or `~astropy.units.Quantity`
 or [MET]_ floats and where one needs to convert between those and what to watch
 out for.
 
-.. _energy_handling_gammapy:
-
-Energy handling in Gammapy
-==========================
-
-Basics
-------
-
-Most objects in Astronomy require an energy axis, e.g. counts spectra or
-effective area tables. In general, this axis can be defined in two ways.
-
-* As an array of energy values. E.g. the Fermi-LAT diffuse flux is given at
-  certain energies and those are stored in an ENERGY FITS table extension.
-  In Gammalib this is represented by GEnergy.
-* As an array of energy bin edges. This is usually stored in EBOUNDS tables,
-  e.g. for Fermi-LAT counts cubes. In Gammalib this is represented by GEbounds.
-
-In Gammapy both the use cases are handled by two seperate classes:
-`gammapy.utils.energy.Energy` for energy values and
-`gammapy.utils.energy.EnergyBounds` for energy bin edges
-
-Energy
-------
-
-The Energy class is a subclass of `~astropy.units.Quantity` and thus has the
-same functionality plus some convenience functions for fits I/O
-
-.. code-block:: python
-
-    >>> from gammapy.utils.energy import Energy
-    >>> energy = Energy([1,2,3], 'TeV')
-    >>> hdu = energy.to_fits()
-    >>> type(hdu)
-    <class 'astropy.io.fits.hdu.table.BinTableHDU'>
-
-EnergyBounds
-------------
-
-The EnergyBounds class is a subclass of Energy. Additional functions are
-available e.g. to compute the bin centers
-
-.. code-block:: python
-
-    >>> from gammapy.utils.energy import EnergyBounds
-    >>> ebounds = EnergyBounds.equal_log_spacing(1, 10, 8, 'GeV')
-    >>> ebounds.size
-    9
-    >>> ebounds.nbins
-    8
-    >>> center = ebounds.log_centers
-    >>> center
-    <Energy [ 1.15478198, 1.53992653, 2.05352503, 2.73841963, 3.65174127,
-              4.86967525, 6.49381632, 8.65964323] GeV>
 
 Reference/API
 =============
 
-.. automodapi:: gammapy.utils.energy
-    :no-inheritance-diagram:
-    :include-all-objects:
-
 .. automodapi:: gammapy.utils.units
     :no-inheritance-diagram:
     :include-all-objects:

gammapy/catalog/fermi.py

@@ -6,7 +6,6 @@
 from astropy.table import Table, Column
 from astropy.time import Time
 from ..utils.scripts import make_path
-from ..utils.energy import EnergyBounds
 from ..utils.table import table_standardise_units_inplace
 from ..maps import Map
 from ..spectrum import FluxPoints
@@ -48,7 +47,7 @@ class SourceCatalogObject3FGL(SourceCatalogObject):
     Catalog is represented by `~gammapy.catalog.SourceCatalog3FGL`.
     """
 
-    _ebounds = EnergyBounds([100, 300, 1000, 3000, 10000, 100000], "MeV")
+    _ebounds = u.Quantity([100, 300, 1000, 3000, 10000, 100000], "MeV")
     _ebounds_suffix = ["100_300", "300_1000", "1000_3000", "3000_10000", "10000_100000"]
     energy_range = u.Quantity([100, 100000], "MeV")
     """Energy range of the catalog.
@@ -374,10 +373,9 @@ def flux_points(self):
         """Flux points (`~gammapy.spectrum.FluxPoints`)."""
         table = Table()
         table.meta["SED_TYPE"] = "flux"
-        e_ref = self._ebounds.log_centers
-        table["e_ref"] = e_ref
-        table["e_min"] = self._ebounds.lower_bounds
-        table["e_max"] = self._ebounds.upper_bounds
+
+        table["e_min"] = self._ebounds[:-1]
+        table["e_max"] = self._ebounds[1:]
 
         flux = self._get_flux_values("Flux")
         flux_err = self._get_flux_values("Unc_Flux")
@@ -462,7 +460,7 @@ class SourceCatalogObject1FHL(SourceCatalogObject):
     Catalog is represented by `~gammapy.catalog.SourceCatalog1FHL`.
     """
 
-    _ebounds = EnergyBounds([10, 30, 100, 500], "GeV")
+    _ebounds = u.Quantity([10, 30, 100, 500], "GeV")
     _ebounds_suffix = ["10_30", "30_100", "100_500"]
     energy_range = u.Quantity([0.01, 0.5], "TeV")
     """Energy range of the Fermi 1FHL source catalog"""
@@ -499,8 +497,8 @@ def flux_points(self):
         """Integral flux points (`~gammapy.spectrum.FluxPoints`)."""
         table = Table()
         table.meta["SED_TYPE"] = "flux"
-        table["e_min"] = self._ebounds.lower_bounds
-        table["e_max"] = self._ebounds.upper_bounds
+        table["e_min"] = self._ebounds[:-1]
+        table["e_max"] = self._ebounds[1:]
         table["flux"] = self._get_flux_values("Flux")
         flux_err = self._get_flux_values("Unc_Flux")
         table["flux_errn"] = np.abs(flux_err[:, 0])
@@ -534,7 +532,7 @@ class SourceCatalogObject2FHL(SourceCatalogObject):
     Catalog is represented by `~gammapy.catalog.SourceCatalog2FHL`.
     """
 
-    _ebounds = EnergyBounds([50, 171, 585, 2000], "GeV")
+    _ebounds = u.Quantity([50, 171, 585, 2000], "GeV")
     _ebounds_suffix = ["50_171", "171_585", "585_2000"]
     energy_range = u.Quantity([0.05, 2], "TeV")
     """Energy range of the Fermi 2FHL source catalog"""
@@ -571,8 +569,8 @@ def flux_points(self):
         """Integral flux points (`~gammapy.spectrum.FluxPoints`)."""
         table = Table()
         table.meta["SED_TYPE"] = "flux"
-        table["e_min"] = self._ebounds.lower_bounds
-        table["e_max"] = self._ebounds.upper_bounds
+        table["e_min"] = self._ebounds[:-1]
+        table["e_max"] = self._ebounds[1:]
         table["flux"] = self._get_flux_values("Flux")
         flux_err = self._get_flux_values("Unc_Flux")
         table["flux_errn"] = np.abs(flux_err[:, 0])
@@ -611,7 +609,7 @@ class SourceCatalogObject3FHL(SourceCatalogObject):
     energy_range = u.Quantity([0.01, 2], "TeV")
     """Energy range of the Fermi 1FHL source catalog"""
 
-    _ebounds = EnergyBounds([10, 20, 50, 150, 500, 2000], "GeV")
+    _ebounds = u.Quantity([10, 20, 50, 150, 500, 2000], "GeV")
 
     def __str__(self):
         return self.info()
@@ -832,10 +830,8 @@ def flux_points(self):
         """Flux points (`~gammapy.spectrum.FluxPoints`)."""
         table = Table()
         table.meta["SED_TYPE"] = "flux"
-        e_ref = self._ebounds.log_centers
-        table["e_ref"] = e_ref
-        table["e_min"] = self._ebounds.lower_bounds
-        table["e_max"] = self._ebounds.upper_bounds
+        table["e_min"] = self._ebounds[:-1]
+        table["e_max"] = self._ebounds[1:]
 
         flux = self.data["Flux_Band"]
         flux_err = self.data["Unc_Flux_Band"]

gammapy/data/event_list.py

@@ -7,7 +7,7 @@
 from astropy.coordinates.angle_utilities import angular_separation
 from astropy.table import Table
 from astropy.table import vstack as vstack_tables
-from ..utils.energy import EnergyBounds
+from ..utils.energy import energy_logspace
 from ..utils.fits import earth_location_from_dict
 from ..utils.scripts import make_path
 from ..utils.time import time_ref_from_dict
@@ -408,7 +408,7 @@ def select_parameter(self, parameter, band):
 
     def _default_plot_ebounds(self):
         energy = self.energy
-        return EnergyBounds.equal_log_spacing(energy.min(), energy.max(), 50)
+        return energy_logspace(energy.min(), energy.max(), 50)
 
     def _counts_spectrum(self, ebounds):
         from ..spectrum import CountsSpectrum

gammapy/irf/effective_area.py

@@ -7,7 +7,7 @@
 from ..utils.nddata import NDDataArray
 from ..maps import MapAxis
 from ..maps.utils import edges_from_lo_hi
-from ..utils.energy import EnergyBounds
+from ..utils.energy import energy_logcenter
 from ..utils.scripts import make_path
 
 __all__ = ["EffectiveAreaTable", "EffectiveAreaTable2D"]
@@ -139,7 +139,7 @@ def from_parametrization(cls, energy, instrument="HESS"):
         instrument : {'HESS', 'HESS2', 'CTA'}
             Instrument name
         """
-        energy = EnergyBounds(energy)
+        energy = u.Quantity(energy)
         # Put the parameters g in a dictionary.
         # Units: g1 (cm^2), g2 (), g3 (MeV)
         # Note that whereas in the paper the parameter index is 1-based,
@@ -155,7 +155,7 @@ def from_parametrization(cls, energy, instrument="HESS"):
             ss += "Valid instruments: HESS, HESS2, CTA"
             raise ValueError(ss)
 
-        xx = energy.log_centers.to_value("MeV")
+        xx = energy_logcenter(energy).to_value("MeV")
 
         g1 = pars[instrument][0]
         g2 = pars[instrument][1]
@@ -165,7 +165,7 @@ def from_parametrization(cls, energy, instrument="HESS"):
         data = u.Quantity(value, "cm2", copy=False)
 
         return cls(
-            energy_lo=energy.lower_bounds, energy_hi=energy.upper_bounds, data=data
+            energy_lo=energy[:-1], energy_hi=energy[1:], data=data
         )
 
     @classmethod
@@ -441,11 +441,10 @@ def to_effective_area_table(self, offset, energy=None):
         if energy is None:
             energy = self.data.axis("energy").edges
 
-        energy = EnergyBounds(energy)
-        area = self.data.evaluate(offset=offset, energy=energy.log_centers)
+        area = self.data.evaluate(offset=offset, energy=energy_logcenter(energy))
 
         return EffectiveAreaTable(
-            energy_lo=energy.lower_bounds, energy_hi=energy.upper_bounds, data=area
+            energy_lo=energy[:-1], energy_hi=energy[1:], data=area
         )
 
     def plot_energy_dependence(self, ax=None, offset=None, energy=None, **kwargs):

gammapy/irf/energy_dispersion.py

@@ -8,10 +8,10 @@
 from astropy.table import Table
 from ..maps import MapAxis
 from ..maps.utils import edges_from_lo_hi
-from ..utils.energy import EnergyBounds, Energy
 from ..utils.scripts import make_path
 from ..utils.nddata import NDDataArray
 from ..utils.fits import energy_axis_to_ebounds
+from ..utils.energy import energy_logcenter
 
 __all__ = ["EnergyDispersion", "EnergyDispersion2D"]
 
@@ -260,14 +260,19 @@ def from_hdulist(cls, hdulist, hdu1="MATRIX", hdu2="EBOUNDS"):
                     ] = l.field("MATRIX")[m_start : m_start + l.field("N_CHAN")[k]]
                     m_start += l.field("N_CHAN")[k]
 
-        e_reco = EnergyBounds.from_ebounds(ebounds_hdu)
-        e_true = EnergyBounds.from_rmf_matrix(matrix_hdu)
+        unit = ebounds_hdu.header.get("TUNIT2")
+        e_reco_lo = Quantity(ebounds_hdu.data["E_MIN"], unit=unit)
+        e_reco_hi = Quantity(ebounds_hdu.data["E_MAX"], unit=unit)
+
+        unit = matrix_hdu.header.get("TUNIT1")
+        e_true_lo = Quantity(matrix_hdu.data["ENERG_LO"], unit=unit)
+        e_true_hi = Quantity(matrix_hdu.data["ENERG_HI"], unit=unit)
 
         return cls(
-            e_true_lo=e_true.lower_bounds,
-            e_true_hi=e_true.upper_bounds,
-            e_reco_lo=e_reco.lower_bounds,
-            e_reco_hi=e_reco.upper_bounds,
+            e_true_lo=e_true_lo,
+            e_true_hi=e_true_hi,
+            e_reco_lo=e_reco_lo,
+            e_reco_hi=e_reco_hi,
             data=pdf_matrix,
         )
 
@@ -684,7 +689,7 @@ class EnergyDispersion2D:
     Read energy dispersion IRF from disk:
 
     >>> from gammapy.irf import EnergyDispersion2D
-    >>> from gammapy.utils.energy import EnergyBounds
+    >>> from gammapy.utils.energy import energy_logspace
     >>> filename = '$GAMMAPY_DATA/tests/irf/hess/pa/hess_edisp_2d_023523.fits.gz'
     >>> edisp2d = EnergyDispersion2D.read(filename, hdu='ENERGY DISPERSION')
     >>> print(edisp2d)
@@ -698,7 +703,7 @@ class EnergyDispersion2D:
     Create energy dispersion matrix (`~gammapy.irf.EnergyDispersion`)
     for a given field of view offset and energy binning:
 
-    >>> energy = EnergyBounds.equal_log_spacing(0.1,20,60, 'TeV')
+    >>> energy = energy_logspace(0.1,20,60, 'TeV')
     >>> edisp = edisp2d.to_energy_dispersion(offset='1.2 deg', e_reco=energy, e_true=energy)
     >>> print(edisp)
     EnergyDispersion
@@ -782,9 +787,9 @@ def from_gauss(cls, e_true, migra, bias, sigma, offset, pdf_threshold=1e-6):
         pdf_threshold : float, optional
             Zero suppression threshold
         """
-        e_true = EnergyBounds(e_true)
+        e_true = Quantity(e_true)
         # erf does not work with Quantities
-        true = e_true.log_centers.to_value("TeV")
+        true = energy_logcenter(e_true).to_value("TeV")
 
         true2d, migra2d = np.meshgrid(true, migra)
 
@@ -872,9 +877,9 @@ def to_energy_dispersion(self, offset, e_true=None, e_reco=None):
         ----------
         offset : `~astropy.coordinates.Angle`
             Offset
-        e_true : `~gammapy.utils.energy.EnergyBounds`, None
+        e_true : `~astropy.units.Quantity`, None
             True energy axis
-        e_reco : `~gammapy.utils.energy.EnergyBounds`
+        e_reco : `~astropy.units.Quantity`
             Reconstructed energy axis
 
         Returns
@@ -885,11 +890,9 @@ def to_energy_dispersion(self, offset, e_true=None, e_reco=None):
         offset = Angle(offset)
         e_true = self.data.axis("e_true").edges if e_true is None else e_true
         e_reco = self.data.axis("e_true").edges if e_reco is None else e_reco
-        e_true = EnergyBounds(e_true)
-        e_reco = EnergyBounds(e_reco)
 
         data = []
-        for energy in e_true.log_centers:
+        for energy in energy_logcenter(e_true):
             vec = self.get_response(offset=offset, e_true=energy, e_reco=e_reco)
             data.append(vec)
 
@@ -914,9 +917,9 @@ def get_response(self, offset, e_true, e_reco=None, migra_step=5e-3):
 
         Parameters
         ----------
-        e_true : `~gammapy.utils.energy.Energy`
+        e_true : `~astropy.units.Quantity`
             True energy
-        e_reco : `~gammapy.utils.energy.EnergyBounds`, None
+        e_reco : `~astropy.units.Quantity`, None
             Reconstructed energy axis
         offset : `~astropy.coordinates.Angle`
             Offset
@@ -928,14 +931,14 @@ def get_response(self, offset, e_true, e_reco=None, migra_step=5e-3):
         rv : `~numpy.ndarray`
             Redistribution vector
         """
-        e_true = Energy(e_true)
+        e_true = Quantity(e_true)
 
         if e_reco is None:
             # Default: e_reco nodes = migra nodes * e_true nodes
-            e_reco = EnergyBounds(self.data.axis("migra").edges * e_true)
+            e_reco = self.data.axis("migra").edges * e_true
         else:
             # Translate given e_reco binning to migra at bin center
-            e_reco = EnergyBounds(e_reco)
+            e_reco = Quantity(e_reco)
 
         # migration value of e_reco bounds
         migra_e_reco = e_reco / e_true
@@ -991,9 +994,9 @@ def plot_migration(self, ax=None, offset=None, e_true=None, migra=None, **kwargs
         else:
             offset = np.atleast_1d(Angle(offset))
         if e_true is None:
-            e_true = Energy([0.1, 1, 10], "TeV")
+            e_true = Quantity([0.1, 1, 10], "TeV")
         else:
-            e_true = np.atleast_1d(Energy(e_true))
+            e_true = np.atleast_1d(Quantity(e_true))
         migra = self.data.axis("migra").center if migra is None else migra
 
         for ener in e_true: