finish up example class

gammapy/irf/effective_area.py

@@ -1,7 +1,13 @@
-#Licensed under a 3-clause BSD style license - see LICENSE.rst
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
 
 from __future__ import absolute_import, division, print_function, unicode_literals
 from ..utils.nddata import NDDataArray, DataAxis, BinnedDataAxis
+import numpy as np
+import astropy.units as u
+
+__all__ = [
+    'EffectiveArea2D',
+]
 
 class EffectiveArea2D(NDDataArray):
     """2D effective area table
@@ -19,6 +25,9 @@ class EffectiveArea2D(NDDataArray):
         Nodes of Offset axis
     data : `~astropy.units.Quantity`
         Effective area
+    meta : dict
+        Optional meta information,
+        supported: ``low_threshold``, ``high_threshold``
 
     Examples
     --------
@@ -27,11 +36,9 @@ class EffectiveArea2D(NDDataArray):
     >>> from gammapy.irf import EffectiveArea2D
     >>> import astropy.units as u
     >>> import numpy as np
-
     >>> energy = np.logspace(0,1,11) * u.TeV
     >>> offset = np.linspace(0,1,4) * u.deg
     >>> data = np.ones(shape=(10,4)) * u.cm * u.cm
-
     >>> eff_area = EffectiveArea2D(energy=energy, offset=offset, data= data)
     >>> print(eff_area)
     Data array summary info
@@ -46,13 +53,170 @@ class EffectiveArea2D(NDDataArray):
     """Secondary axis: Offset from pointing position"""
     axis_names = ['energy', 'offset']
 
+    @property
+    def low_threshold(self):
+
+        """Low energy threshold"""
+        return self.meta.low_theshold
+
+    @property
+    def high_treshold(self):
+        """High energy threshold"""
+        return self.meta.high_threshold
+
     @classmethod
     def from_table(cls, t):
-        """This is a read for the format specified at
+        """This is a reader for the format specified at
         http://gamma-astro-data-formats.readthedocs.io/en/latest/irfs/effective_area/index.html#aeff-2d-format
         """
-        raise NotImplementedError()
+        # TODO: only define columns here, move actual reader to the base class
+        energy_col = 'ENERG'
+        offset_col = 'THETA'
+        data_col = 'EFFAREA'
+
+        energy_lo = t['{}_LO'.format(energy_col)].quantity[0]
+        energy_hi = t['{}_HI'.format(energy_col)].quantity[0]
+        energy = np.append(energy_lo.value, energy_hi[-1].value) * energy_lo.unit
+        offset = t['{}_HI'.format(offset_col)].quantity[0]
+        # see https://github.com/gammasky/hess-host-analyses/issues/32
+        data = t['{}'.format(data_col)].quantity[0].transpose()
+
+        return cls(offset=offset, energy=energy, data=data)
+
+    def to_effective_area(self):
+        """Create effective area table"""
+        raise NotImplementedError
+
+    def plot_energy_dependence(self, ax=None, offset=None, energy=None, **kwargs):
+        """Plot effective area versus energy for a given offset.
+
+        Parameters
+        ----------
+        ax : `~matplolib.axes`, optional
+            Axis
+        offset : `~astropy.coordinates.Angle`
+            Offset
+        energy : `~astropy.units.Quantity`
+            Energy axis
+        kwargs : dict
+            Forwarded tp plt.plot()
+
+        Returns
+        -------
+        ax : `~matplolib.axes`
+            Axis
+        """
+        import matplotlib.pyplot as plt
+
+        ax = plt.gca() if ax is None else ax
+
+        if offset is None:
+            off_min, off_max = self.offset.nodes[[0, -1]].value
+            offset = np.linspace(off_min, off_max, 4) * self.offset.unit
+
+        if energy is None:
+            energy = self.energy.nodes
+
+        for off in offset:
+            area = self.evaluate(offset=off, energy=energy)
+            label = 'offset = {:.1f}'.format(off)
+            ax.plot(energy, area.value, label=label, **kwargs)
+
+        ax.set_xscale('log')
+        ax.set_xlabel('Energy [{0}]'.format(self.energy.unit))
+        ax.set_ylabel('Effective Area [{0}]'.format(self.data.unit))
+        ax.set_xlim(min(energy.value), max(energy.value))
+        ax.legend(loc='upper left')
+
+        return ax
+
+    def plot_offset_dependence(self, ax=None, offset=None, energy=None, **kwargs):
+        """Plot effective area versus offset for a given energy
+
+        Parameters
+        ----------
+        ax : `~matplolib.axes`, optional
+            Axis
+        offset : `~astropy.coordinates.Angle`
+            Offset axis
+        energy : `~gammapy.utils.energy.Energy`
+            Energy
+
+        Returns
+        -------
+        ax : `~matplolib.axes`
+            Axis
+        """
+        import matplotlib.pyplot as plt
+
+        ax = plt.gca() if ax is None else ax
+
+        if energy is None:
+            e_min, e_max = np.log10(self.energy.nodes[[0, -1]].value)
+            energy = np.logspace(e_min, e_max, 4) * self.energy.unit
+
+        if offset is None:
+            off_lo, off_hi = self.offset.nodes[[0, -1]].to('deg').value
+            offset = np.linspace(off_lo, off_hi, 100) * u.deg
+
+        for ee in energy:
+            area = self.evaluate(offset=offset, energy=ee)
+            area /= np.nanmax(area)
+            if np.isnan(area).all():
+                continue
+            label = 'energy = {:.1f}'.format(ee)
+            ax.plot(offset, area, label=label, **kwargs)
+
+        ax.set_ylim(0, 1.1)
+        ax.set_xlabel('Offset ({0})'.format(self.offset.unit))
+        ax.set_ylabel('Relative Effective Area')
+        ax.legend(loc='best')
+
+        return ax
+
+    def plot_image(self, ax=None, offset=None, energy=None, **kwargs):
+        """Plot effective area 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:
+            vals = self.offset.nodes.value
+            offset = np.linspace(vals.min(), vals.max(), 100)
+            offset = offset * self.offset.unit
+
+        if energy is None:
+            vals = np.log10(self.energy.nodes.value)
+            energy = np.logspace(vals.min(), vals.max(), 100) * self.energy.unit
+
+        aeff = self.evaluate(offset=offset, energy=energy)
+        extent = [
+            offset.value.min(), offset.value.max(),
+            energy.value.min(), energy.value.max(),
+        ]
+        ax.imshow(aeff.value, extent=extent, **kwargs)
+
+        ax.set_yscale('log')
+        ax.set_xlabel('Offset ({0})'.format(offset.unit))
+        ax.set_ylabel('Energy ({0})'.format(energy.unit))
+
+        ax.set_title('Effective Area ({0})'.format(aeff.unit))
+
+        ax.legend()
+
+        return ax
 
-    def plot():
-        """Plot image"""
-        raise NotImplementedError()
+    def peek(self, figsize=(15, 5)):
+        """Quick-look summary plots."""
+        import matplotlib.pyplot as plt
+        fig, axes = plt.subplots(nrows=1, ncols=3, figsize=figsize)
+        self.plot_image(ax=axes[0])
+        self.plot_energy_dependence(ax=axes[1])
+        self.plot_offset_dependence(ax=axes[2])
+        plt.tight_layout()
+        plt.show()

gammapy/irf/tests/test_effective_area2.py

@@ -17,11 +17,19 @@ def test_EffectiveArea2D_generic():
     energy = np.logspace(0, 1, 4) * u.TeV
     offset = [0.2, 0.3] * u.deg
     effective_area = np.arange(6).reshape(3, 2) * u.cm * u.cm
-    aeff = EffectiveArea2D(offset=offset, energy=energy, data=effective_area)
-
+    meta = dict(name = 'example')
+    aeff = EffectiveArea2D(offset=offset, energy=energy, data=effective_area,
+                           meta=meta)
     assert (aeff.axes[0].data == energy).all()
     assert (aeff.axes[1].data == offset).all()
     assert (aeff.data == effective_area).all()
+    assert aeff.meta.name == 'example'
+
+    wrong_data = np.arange(8).reshape(4,2) * u.cm * u.cm
+
+    with pytest.raises(ValueError):
+        aeff.data = wrong_data
+        aeff = EffectiveArea(offset=offset, energy=energy, data=wrong_data)
 
     # Test evaluate function 
     # Check that nodes are evaluated correctly
@@ -55,6 +63,7 @@ def test_EffectiveArea2D_generic():
     assert_equal(actual, desired)
 
     # Case 2: offset = scalar, energy = 1Darray
+
     offset = 0.564 * u.deg
     nbins = 10 
     energy = np.logspace(3, 4, nbins) * u.GeV
@@ -98,30 +107,40 @@ def test_EffectiveArea2D_generic():
     print(aeff)
 
 
-@pytest.mark.xfail()
+@requires_dependency('scipy')
+@requires_dependency('matplotlib')
 @requires_data('gammapy-extra')
 def test_EffectiveArea2D(tmpdir):
 
     filename = gammapy_extra.filename('datasets/hess-crab4-hd-hap-prod2/run023400-023599/run023523/hess_aeff_2d_023523.fits.gz')
     aeff = EffectiveArea2D.read(filename)
 
+    assert aeff.energy.nbins == 73
+    assert aeff.offset.nbins == 6
+    assert aeff.data.shape == (73, 6)
 
+    assert aeff.energy.unit == 'TeV'
+    assert aeff.offset.unit == 'deg'
+    assert aeff.data.unit == 'm2'
+
+    aeff.plot_image()
+    aeff.plot_energy_dependence()
+    aeff.plot_offset_dependence()
 
     # Test ARF export
-    offset = Angle(0.236, 'deg')
-    e_axis = Quantity(np.logspace(0, 1, 20), 'TeV')
-    effareafrom2d = aeff.to_effective_area_table(offset, e_axis)
-    energy = EnergyBounds(e_axis).log_centers
-    area = aeff.evaluate(offset, energy)
-    effarea1d = EffectiveAreaTable(e_axis, area)
-    test_energy = Quantity(2.34, 'TeV')
-    actual = effareafrom2d.evaluate(test_energy)
-    desired = effarea1d.evaluate(test_energy)
-    assert_equal(actual, desired)
+    # offset = Angle(0.236, 'deg')
+    # e_axis = Quantity(np.logspace(0, 1, 20), 'TeV')
+    # effareafrom2d = aeff.to_effective_area_table(offset, e_axis)
+    # energy = EnergyBounds(e_axis).log_centers
+    # area = aeff.evaluate(offset, energy)
+    # effarea1d = EffectiveAreaTable(e_axis, area)
+    # test_energy = Quantity(2.34, 'TeV')
+    # actual = effareafrom2d.evaluate(test_energy)
+    # desired = effarea1d.evaluate(test_energy)
+    # assert_equal(actual, desired)
 
     # Test ARF export #2
-    effareafrom2dv2 = aeff.to_effective_area_table('1.2 deg')
-    actual = effareafrom2dv2.effective_area
-    desired = aeff.evaluate(offset='1.2 deg')
-    assert_equal(actual, desired)
-
+    # effareafrom2dv2 = aeff.to_effective_area_table('1.2 deg')
+    # actual = effareafrom2dv2.effective_area
+    # desired = aeff.evaluate(offset='1.2 deg')
+    # assert_equal(actual, desired)

gammapy/utils/nddata.py

@@ -4,6 +4,7 @@
 import itertools
 import numpy as np
 import abc
+from ..extern.bunch import Bunch
 from astropy.units import Quantity
 from astropy.table import Table, Column
 from astropy.extern import six
@@ -37,7 +38,11 @@ def __init__(self, **kwargs):
 
         # TODO : Dynamically generate function signature
         # https://github.com/astropy/astropy/blob/ffc0a89b2c42fd440eb19bcb2f93db90cab3c98b/astropy/utils/codegen.py#L30
-        self._data = kwargs.pop('data', None)
+        data = kwargs.pop('data', None)
+
+        meta = kwargs.pop('meta', None)
+        if meta is not None:
+            self.meta = Bunch(meta)
 
         self._axes = list()
         for axis_name in self.axis_names:
@@ -46,6 +51,7 @@ def __init__(self, **kwargs):
             axis.data = Quantity(value)
             self._axes.append(axis)
 
+        self.data = data
         self._regular_grid_interp = None
 
     @property
@@ -69,28 +75,27 @@ def data(self, data):
         data : `~astropy.units.Quantity`, array-like
             Data array
         """
-        data = np.array(data)
         dimension = len(data.shape)
         if dimension != self.dim:
             raise ValueError('Overall dimensions to not match. '
                              'Data: {}, Hist: {}'.format(dimension, self.dim))
 
         for dim in np.arange(self.dim):
-            if self.axes[dim].nbins != data.shape[dim]:
-                axis = self.axes[dim]
+            axis = self.axes[dim]
+            if axis.nbins != data.shape[dim]:
                 msg = 'Data shape does not match in dimension {d}\n'
-                msg += 'Axis "{n}":{sa}, Data {sd}'
-                raise ValueError(msg.format(d=dim, n=axis.name, sa=axis.nbins,
-                                            sd=data.shape[dim]))
-        self._data = Quantity(data)
+                msg += 'Axis {n} : {sa}, Data {sd}'
+                raise ValueError(msg.format(d=dim, n=self.axis_names[dim],
+                                            sa=axis.nbins, sd=data.shape[dim]))
+        self._data = data
 
     @property
     def dim(self):
         """Dimension (number of axes)"""
         return len(self.axes)
 
     def to_table(self):
-        """Convert to astropy.Table"""
+        """Convert to `~astropy.table.Table`"""
 
         raise NotImplementedError("Broken")