Add EnergyOffsetArray

docs/background/energy_offset_array.rst

@@ -0,0 +1,24 @@
+.. _energyoffset_array:
+
+EnergyOffset Array
+==================
+
+The `~gammapy.background.EnergyOffsetArray` class represents a 2D array *(energy,offset)* that is filled with an eventlist. For a set of observations, by giving an energy binning and an offset binning, you fill the events in this histogram.
+
+
+Four Crab observations are located in the ``gammapy-extra`` repository as
+examples:
+`hess_events_simulated_023523.fits`_ , `hess_events_simulated_023526.fits`_ , `hess_events_simulated_023559.fits`_ and `hess_events_simulated_023592.fits`_
+
+An example script of how to fill this Array from these four observations and plots the result is given in the ``examples`` directory:
+:download:`example_energy_offset_array.py <../../examples/example_energy_offset_array.py>`
+
+.. plot:: ../examples/example_energy_offset_array.py
+   :include-source:
+
+
+
+.. _hess_events_simulated_023523.fits: https://github.com/gammapy/gammapy-extra/tree/master/datasets/hess-crab4/hess_events_simulated_023523.fits
+.. _hess_events_simulated_023526.fits: https://github.com/gammapy/gammapy-extra/tree/master/datasets/hess-crab4/hess_events_simulated_023526.fits
+.. _hess_events_simulated_023559.fits: https://github.com/gammapy/gammapy-extra/tree/master/datasets/hess-crab4/hess_events_simulated_023559.fits
+.. _hess_events_simulated_023592.fits: https://github.com/gammapy/gammapy-extra/tree/master/datasets/hess-crab4/hess_events_simulated_023592.fits

docs/background/index.rst

@@ -39,7 +39,8 @@ If you'd like to learn more about using `gammapy.background`, read the following
 
    models
    make_models
-
+   energy_offset_array
+
 Reference/API
 =============
 

examples/example_energy_offset_array.py

@@ -0,0 +1,26 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+from gammapy.data import DataStore
+from gammapy.datasets import gammapy_extra
+from gammapy.background import EnergyOffsetArray
+from gammapy.utils.energy import EnergyBounds
+
+
+def make_counts_array():
+    """Make an example counts array with energy and offset axes."""
+    data_store = DataStore.from_dir(gammapy_extra.dir / 'datasets/hess-crab4')
+
+    event_lists = data_store.load_all('events')
+    ebounds = EnergyBounds.equal_log_spacing(0.1, 100, 100, 'TeV')
+    offset = np.linspace(0, 2.5, 100)
+    array = EnergyOffsetArray(ebounds, offset)
+    array.fill_events(event_lists)
+
+    return array
+
+
+if __name__ == '__main__':
+    array = make_counts_array()
+    array.plot_image()
+    plt.show()

gammapy/background/__init__.py

@@ -2,6 +2,7 @@
 """
 Background estimation and modeling methods.
 """
+from .energy_offset_array import *
 from .fov import *
 from .maps import *
 from .cube import *

gammapy/background/energy_offset_array.py

@@ -0,0 +1,106 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+from __future__ import absolute_import, division, print_function, unicode_literals
+import numpy as np
+from astropy.coordinates import Angle
+from astropy.units import Quantity
+
+__all__ = [
+    'EnergyOffsetArray',
+]
+
+
+class EnergyOffsetArray(object):
+    """
+    Energy offset dependent array.
+
+    Parameters
+    ----------
+    energy : `~astropy.units.Quantity`
+         energy bin vector
+    offset : `~astropy.coordinates.Angle`
+        offset bin vector
+    data : `~numpy.ndarray`
+        data array (2D)
+
+    """
+
+    def __init__(self, energy, offset, data=None):
+        self.energy = Quantity(energy, 'TeV')
+        self.offset = Angle(offset, 'deg')
+        if data is None:
+            self.data = np.zeros((len(energy) - 1, len(offset) - 1))
+        else:
+            self.data = data
+
+    def fill_events(self, event_lists):
+        """Fill events histogram.
+
+        This add the counts to the existing value array.
+
+        Parameters
+        -------------
+        event_lists : list of `~gammapy.data.EventList`
+           Python list of event list objects.
+
+
+        """
+        # loop over the Lost of object EventList
+        for event_list in event_lists:
+            # Fill the events
+            counts = self._fill_one_event_list(event_list)
+            self.data += counts
+
+    def _fill_one_event_list(self, events):
+        """
+        histogram the counts of an EventList object in 2D (energy,offset)
+
+        Parameters
+        -------------
+        events :`~gammapy.data.EventList`
+           Event list objects.
+
+        """
+
+        offset = events.offset
+        ev_energy = events.energy
+
+        # stack the offset and energy array
+        ev_cube_array = np.vstack([ev_energy, offset]).T
+
+        # fill data cube into histogramdd
+        ev_cube_hist, ev_cube_edges = np.histogramdd(ev_cube_array,
+                                                     [self.energy.value,
+                                                      self.offset.value])
+        return ev_cube_hist
+
+    def plot_image(self, ax=None, offset=None, energy=None, **kwargs):
+        """
+        Plot Energy_offset Array image (x=offset, y=energy).
+        """
+        import matplotlib.pyplot as plt
+
+        kwargs.setdefault('cmap', 'afmhot')
+        kwargs.setdefault('origin', 'bottom')
+        kwargs.setdefault('interpolation', 'nearest')
+
+        ax = plt.gca() if ax is None else ax
+
+        if offset is None:
+            offset = self.offset
+
+        if energy is None:
+            energy = self.energy
+
+        extent = [
+            offset.value.min(), offset.value.max(),
+            energy.value.min(), energy.value.max(),
+        ]
+        ax.imshow(self.data, extent=extent, **kwargs)
+        ax.semilogy()
+        ax.set_xlabel('Offset ({0})'.format(offset.unit))
+        ax.set_ylabel('Energy ({0})'.format(energy.unit))
+        ax.set_title('Energy_offset Array')
+        ax.legend()
+        image = ax.imshow(self.data, extent=extent, **kwargs)
+        plt.colorbar(image)
+        return ax

gammapy/background/tests/test_energy_offset_array.py

@@ -0,0 +1,21 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+from __future__ import absolute_import, division, print_function, unicode_literals
+import numpy as np
+from ..energy_offset_array import EnergyOffsetArray
+from ...data import DataStore
+from ...datasets import gammapy_extra
+from ...utils.testing import requires_dependency, requires_data
+from ...utils.energy import EnergyBounds
+
+
+@requires_dependency('matplotlib')
+@requires_data('gammapy-extra')
+def test_energy_offset_array():
+    dir = str(gammapy_extra.dir) + '/datasets/hess-crab4'
+    data_store = DataStore.from_dir(dir)
+    ev_list = data_store.load_all('events')
+    ebounds = EnergyBounds.equal_log_spacing(0.1, 100, 100, 'TeV')
+    offset = np.linspace(0, 2.5, 100)
+    array = EnergyOffsetArray(ebounds, offset)
+    array.fill_events(ev_list)
+    array.plot_image()

gammapy/data/data_store.py

@@ -209,6 +209,25 @@ def load(self, obs_id, filetype):
         else:
             raise ValueError('Invalid filetype.')
 
+    def load_all(self, filetype):
+        """Load a given file type for all observations
+
+        Parameters
+        ----------
+        filetype : {'events', 'aeff', 'edisp', 'psf', 'background'}
+            Type of file.
+
+        Returns
+        -------
+        list : python list of object
+            Object depends on type, e.g. for `events` it is a list of `~gammapy.data.EventList`.
+        """
+        data_lists = []
+        for obs_id in self.obs_table['OBS_ID']:
+            data_list = self.load(obs_id, filetype)
+            data_lists.append(data_list)
+        return data_lists
+
     def check_integrity(self, logger):
         """Check integrity, i.e. whether index table and files match.
         """

gammapy/data/event_list.py

@@ -133,6 +133,14 @@ def altaz(self):
         lon, lat = self['AZ'], self['ALT']
         return SkyCoord(lon, lat, unit='deg', frame=altaz_frame)
 
+    @property
+    def offset(self):
+        """Event offset (`~astropy.coordinates.Angle`)"""
+        position = self.radec
+        center = self.pointing_radec
+        offset = center.separation(position)
+        return Angle(offset, unit='deg')
+
     @property
     def energy(self):
         """Event energies (`~astropy.units.Quantity`)"""
@@ -620,14 +628,14 @@ class EventListDatasetChecker(object):
         Logger to use (use module-level Gammapy logger by default)
     """
     _AVAILABLE_CHECKS = OrderedDict(
-            misc='check_misc',
-            times='check_times',
-            coordinates='check_coordinates',
+        misc='check_misc',
+        times='check_times',
+        coordinates='check_coordinates',
     )
 
     accuracy = OrderedDict(
-            angle=Angle('1 arcsec'),
-            time=Quantity(1, 'microsecond'),
+        angle=Angle('1 arcsec'),
+        time=Quantity(1, 'microsecond'),
 
     )
 
@@ -712,9 +720,9 @@ def check_times(self):
 
         # http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_Data/Time_in_ScienceTools.html
         telescope_met_refs = OrderedDict(
-                FERMI=Time('2001-01-01T00:00:00'),
-                HESS=Time('2000-01-01T12:00:00.000'),
-                # TODO: Once CTA has specified their MET reference add check here
+            FERMI=Time('2001-01-01T00:00:00'),
+            HESS=Time('2000-01-01T12:00:00.000'),
+            # TODO: Once CTA has specified their MET reference add check here
         )
 
         telescope = self.dset.event_list.meta['TELESCOP']

gammapy/data/tests/test_data_store.py

@@ -4,7 +4,8 @@
 from ...data import EventList
 from ...data import DataStore
 from ...utils.testing import requires_data, data_manager
-
+from ...datasets import gammapy_extra
+from numpy.testing import assert_allclose
 
 @pytest.mark.xfail
 @requires_data('hess')
@@ -59,6 +60,15 @@ def test_DataStore_load(data_manager):
     assert isinstance(events, EventList)
 
 
+@requires_data('hess')
+def test_DataStore_load_all():
+    """Test loading data and IRF files via the DataStore"""
+    dir = str(gammapy_extra.dir) + '/datasets/hess-crab4'
+    data_store = DataStore.from_dir(dir)
+    event_lists = data_store.load_all(filetype='events')
+    assert_allclose(event_lists[0]['ENERGY'][0], 1.1156039)
+    assert_allclose(event_lists[-1]['ENERGY'][0], 1.0204216)
+
 @pytest.mark.xfail
 @requires_data('hess')
 def test_DataStore_other(data_manager):

gammapy/data/tests/test_event_list.py

@@ -19,6 +19,7 @@ def test_EventList():
     assert event_list.radec[0].to_string() == '82.7068 19.8186'
     assert event_list.galactic[0].to_string() == '185.956 -7.69277'
     assert event_list.altaz[0].to_string() == '46.2059 31.2001'
+    assert_allclose(event_list.offset[0].value, 1.904497742652893, rtol=1e-5)
     assert '{:1.5f}'.format(event_list.energy[0]) == '11.64355 TeV'
 
     lon, lat, height = event_list.observatory_earth_location.to_geodetic()