Add function to compute exposure cubes

gammapy/data/spectral_cube.py

@@ -157,6 +157,31 @@ def read(cls, filename):
 
         return cls(data, wcs, energy)
 
+    @classmethod
+    def read_counts(cls, filename):
+        """Read spectral cube of type "counts" from FITS file.
+
+        Parameters
+        ----------
+        filename : str
+            File name
+
+        Returns
+        -------
+        spectral_cube : `SpectralCube`
+            Spectral cube
+        """
+        data = fits.getdata(filename)
+        data = Quantity(data, '')
+        # Note: the energy axis of the FITS cube is unusable.
+        # We only use proj for LON, LAT and do ENERGY ourselves
+        header = fits.getheader(filename)
+        wcs = WCS(header)
+        energy=EnergyBounds.from_ebounds(fits.open(filename)['EBOUNDS'],unit='keV')
+        #keV to match convention in Fermi ST and CTOOLS, TODO read from header
+
+        return cls(data, wcs, energy)
+
     def world2pix(self, lon, lat, energy, combine=False):
         """Convert world to pixel coordinates.
 

gammapy/irf/__init__.py

@@ -6,3 +6,4 @@
 from .psf_table import *
 from .psf_analytical import *
 from .energy_dispersion import *
+from .exposure import *

gammapy/irf/effective_area_table.py

@@ -507,7 +507,7 @@ def evaluate(self, offset=None, energy=None):
         """Evaluate effective area for a given energy and offset.
 
         If a parameter is not given, the nodes from the FITS table are used.
-        2D input arrays are not supported yet.
+        Both offset and energy can be arbitrary ndarrays
 
         Parameters
         ----------
@@ -549,22 +549,25 @@ def _eval_linear(self, offset=None, energy=None):
 
         Parameters
         ----------
-        offset : float
+        offset : float or ndarray of floats
             offset in deg
-        energy : float
+        energy : float or ndarray of floats
             energy in TeV
 
         Returns
         -------
-        eff_area : float
+        eff_area : float or ndarray of floats
             Effective area
         """
         off = np.atleast_1d(offset)
         ener = np.atleast_1d(energy)
-        points = [(x,y) for x in off for y in ener]
-        shape = (off.size, ener.size)
+        shape_off = np.shape(off)
+        shape_ener = np.shape(ener)
+        off = off.flatten()
+        ener = ener.flatten()
+        points = [(x, y) for x in off for y in ener]
         val = self._linear(points)
-        val = np.reshape(val, shape).squeeze()
+        val = np.reshape(val, np.append(shape_off, shape_ener)).squeeze()
 
         return val
 
@@ -691,7 +694,6 @@ def info(self):
 
         return ss
 
-
     def _prepare_linear_interpolator(self):
         """Setup `~scipy.interpolate.RegularGridInterpolator`
         """
@@ -702,7 +704,7 @@ def _prepare_linear_interpolator(self):
         y = np.log10(self.energy.value)
         vals = self.eff_area.value
 
-        self._linear = RegularGridInterpolator((x,y),vals, bounds_error = True)
+        self._linear = RegularGridInterpolator((x, y), vals, bounds_error=True)
 
     def _prepare_spline_interpolator(self):
         """Only works for radial symmetric input files (N=2)

gammapy/irf/exposure.py

@@ -0,0 +1,50 @@
+# 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 ..data import SpectralCube
+from astropy.coordinates import SkyCoord, Angle
+
+__all__ = [
+    'exposure_cube'
+]
+
+
+def exposure_cube(pointing,
+                  livetime,
+                  aeff2D,
+                  ref_cube):
+    """Calculate exposure cube
+
+    Parameters
+    ----------
+    pointing : `~astropy.coordinates.SkyCoord`
+        pointing direction
+    livetime : `~astropy.units.Quantity`
+        livetime
+    aeff2D : `~gammapy.irf.EffectiveAreaTable2D`
+        effective area table
+    ref_cube : `~gammapy.data.SpectralCube`
+        reference cube used to define geometry
+
+    Returns
+    -------
+    expcube : `~gammapy.data.SpectralCube`
+        3D exposure
+    """
+    ny, nx = ref_cube.data.shape[1:]
+    xx, yy = np.meshgrid(np.arange(nx), np.arange(ny))
+    lon, lat, en = ref_cube.pix2world(xx, yy, 0)
+    coord = SkyCoord(lon, lat, frame=ref_cube.wcs.wcs.radesys.lower())  # don't care about energy
+    offset = coord.separation(pointing)
+    ###this is a hack to make the current code work in spite of issues with offset binning from HAP export
+    ###TODO remove once the issue has been fixed
+    offset = np.clip(offset, Angle(0.2, 'deg'), Angle(2.2, 'deg'))
+
+    exposure = aeff2D.evaluate(offset, ref_cube.energy)
+    exposure = np.rollaxis(exposure, 2)
+    exposure *= livetime
+
+    expcube = SpectralCube(data=exposure,
+                           wcs=ref_cube.wcs,
+                           energy=ref_cube.energy)
+    return expcube

gammapy/irf/tests/test_effective_area.py

@@ -121,6 +121,17 @@ def test_EffectiveAreaTable2D(method):
     desired = np.zeros([nbinso, nbinse]).shape
     assert_equal(actual, desired)
 
+    # case 6: offset = 2Darray, energy = 1Darray
+    nbinse = 16
+    nx, ny = (12,3)
+    offset=np.linspace(1,0,nx*ny).reshape(nx,ny)
+    offset = Angle(offset, 'deg')
+    energy = Quantity(np.logspace(0, 1, nbinse), 'TeV')
+    actual = aeff.evaluate(offset=offset, energy=energy).shape
+    desired = np.zeros([nx, ny, nbinse]).shape
+    assert_equal(actual, desired)
+
+
     # Test ARF export
     offset = Angle(0.236, 'deg')
     e_axis = Quantity(np.logspace(0, 1, 20), 'TeV')

gammapy/irf/tests/test_exposure.py

@@ -0,0 +1,33 @@
+from __future__ import absolute_import, division, print_function, unicode_literals
+import numpy as np
+from numpy.testing import assert_allclose, assert_equal
+from astropy.units import Quantity
+from astropy.coordinates import SkyCoord
+from ...utils.testing import requires_dependency
+from ...utils.testing import requires_data
+from ...irf import exposure_cube
+from ...data import SpectralCube
+from ...irf import EffectiveAreaTable2D
+from ...datasets import gammapy_extra
+
+
+@requires_dependency('scipy')
+@requires_data('gammapy-extra')
+def test_exposure_cube():
+    exp_ref = Quantity(4.7e8, 'm^2 s')
+
+    aeff_filename = gammapy_extra.filename("datasets/hess-crab4/hess_aeff_023523.fits.gz")
+    ccube_filename = gammapy_extra.filename("datasets/hess-crab4/hess_events_simulated_023523_cntcube.fits")
+
+    pointing = SkyCoord(83.633, 21.514, frame='fk5', unit='deg')
+    livetime = Quantity(1581.17, 's')
+    aeff2D = EffectiveAreaTable2D.read(aeff_filename)
+    count_cube = SpectralCube.read_counts(ccube_filename)
+    exp_cube = exposure_cube(pointing, livetime, aeff2D, count_cube)
+
+    assert np.shape(exp_cube.data)[1:] == np.shape(count_cube.data)[1:]
+    assert np.shape(exp_cube.data)[0] == np.shape(count_cube.data)[0] + 1
+    assert exp_cube.wcs == count_cube.wcs
+    assert_equal(count_cube.energy, exp_cube.energy)
+    assert_allclose(exp_cube.data, exp_ref, rtol=100)
+    assert exp_cube.data.unit == exp_ref.unit