Add make_images method in MapMaker

gammapy/cube/exposure.py

@@ -1,6 +1,7 @@
 # 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 ..spectrum.models import PowerLaw
 from ..maps import WcsNDMap
 
 __all__ = [
@@ -42,3 +43,45 @@ def make_map_exposure_true_energy(pointing, livetime, aeff, geom):
     exposure = (exposure * livetime).to('m2 s')
 
     return WcsNDMap(geom, exposure.value, unit=exposure.unit)
+
+
+def _map_spectrum_weight(map, spectrum=None):
+    """Weight a map with a spectrum.
+
+    This requires map to have an "energy" axis.
+    The weights are normalised so that they sum to 1.
+    The mean and unit of the output image is the same as of the input cube.
+
+    At the moment this is used to get a weighted exposure image.
+
+    Parameters
+    ----------
+    map : `~gammapy.maps.Map`
+        Input map with an "energy" axis.
+    spectrum : `~gammapy.spectrum.models.SpectralModel`
+        Spectral model to compute the weights.
+        Default is power-law with spectral index of 2.
+
+    Returns
+    -------
+    map_weighted : `~gammapy.maps.Map`
+        Weighted image
+    """
+    if spectrum is None:
+        spectrum = PowerLaw(index=2.0)
+
+    # Compute weights vector
+    # Should we change to call spectrum.integrate ?
+    energy_axis = map.geom.get_axis_by_name("energy")
+    energy_center = energy_axis.center * energy_axis.unit
+    energy_edges = energy_axis.edges * energy_axis.unit
+    energy_width = np.diff(energy_edges)
+    weights = spectrum(energy_center) * energy_width
+    weights /= weights.sum()
+
+    # Make new map with weights applied
+    map_weighted = map.copy()
+    for img, idx in map_weighted.iter_by_image():
+        img *= weights[idx].value
+
+    return map_weighted

gammapy/cube/make.py

@@ -6,7 +6,7 @@
 from astropy.coordinates import Angle
 from ..maps import Map, WcsGeom
 from .counts import fill_map_counts
-from .exposure import make_map_exposure_true_energy
+from .exposure import make_map_exposure_true_energy, _map_spectrum_weight
 from .background import make_map_background_irf, _fov_background_norm
 
 __all__ = [
@@ -82,10 +82,10 @@ def run(self, obs_list, selection=None):
     def _process_obs(self, obs, selection):
         # Compute cutout geometry and slices to stack results back later
         cutout_map = Map.from_geom(self.geom).cutout(
-                position=obs.pointing_radec,
-                width=2 * self.offset_max,
-                mode='trim',
-            )
+            position=obs.pointing_radec,
+            width=2 * self.offset_max,
+            mode='trim',
+        )
 
         log.info('Processing observation {}'.format(obs.obs_id))
 
@@ -116,6 +116,32 @@ def _process_obs(self, obs, selection):
             data = maps_obs[name].quantity.to(self.maps[name].unit).value
             self.maps[name].fill_by_coord(coords, data)
 
+    def make_images(self, spectrum=None):
+        """Create 2D images by summing over the energy axis.
+
+        Exposure is weighted with an assumed spectrum,
+        resulting in a weighted mean exposure image.
+
+        Parameters
+        ----------
+        spectrum : `~gammapy.spectrum.models.SpectralModel`
+            Spectral model to compute the weights.
+            Default is power-law with spectral index of 2.
+
+        Returns
+        -------
+        images : dict of `~gammapy.maps.Map`
+            2D images
+        """
+        images = {}
+        for name, map in self.maps.items():
+            if name == 'exposure':
+                map = _map_spectrum_weight(map, spectrum)
+
+            images[name] = map.sum_over_axes()
+
+        return images
+
 
 class MapMakerObs(object):
     """Make maps for a single IACT observation.

gammapy/cube/tests/test_exposure.py

@@ -1,12 +1,14 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 from __future__ import absolute_import, division, print_function, unicode_literals
 import pytest
+import numpy as np
 from numpy.testing import assert_allclose
 from astropy.coordinates import SkyCoord
 from ...utils.testing import requires_data
-from ...maps import WcsGeom, HpxGeom, MapAxis
+from ...maps import WcsGeom, HpxGeom, MapAxis, WcsNDMap
 from ...irf import EffectiveAreaTable2D
-from ..exposure import make_map_exposure_true_energy
+from ..exposure import make_map_exposure_true_energy, _map_spectrum_weight
+from ...spectrum.models import ConstantModel
 
 pytest.importorskip('scipy')
 pytest.importorskip('healpy')
@@ -59,3 +61,16 @@ def test_make_map_exposure_true_energy(aeff, pars):
     assert m.data.shape == pars['shape']
     assert m.unit == 'm2 s'
     assert_allclose(m.data.sum(), pars['sum'], rtol=1e-5)
+
+
+def test_map_spectrum_weight():
+    axis = MapAxis.from_edges([0.1, 10, 1000], unit='TeV', name='energy')
+    expo_map = WcsNDMap.create(npix=10, binsz=1, axes=[axis], unit='m2 s')
+    expo_map.data += 1
+    spectrum = ConstantModel(42)
+
+    weighted_expo = _map_spectrum_weight(expo_map, spectrum)
+
+    assert weighted_expo.data.shape == (2, 10, 10)
+    assert weighted_expo.unit == 'm2 s'
+    assert_allclose(weighted_expo.data.sum(), 100)

gammapy/cube/tests/test_make.py

@@ -33,13 +33,15 @@ def geom(ebounds):
         'geom': geom(ebounds=[0.1, 1, 10]),
         'counts': 34366,
         'exposure': 3.99815e+11,
+        'exposure_image': 7.921993e+10,
         'background': 187528.89,
     },
     {
         # Test single energy bin
         'geom': geom(ebounds=[0.1, 10]),
         'counts': 34366,
         'exposure': 1.16866e+11,
+        'exposure_image': 1.16866e+11,
         'background': 1988492.8,
     },
     {
@@ -48,6 +50,7 @@ def geom(ebounds):
         'exclusion_mask': Map.from_geom(geom(ebounds=[0.1, 10])),
         'counts': 34366,
         'exposure': 1.16866e+11,
+        'exposure_image': 1.16866e+11,
         'background': 1988492.8,
     },
 
@@ -57,6 +60,7 @@ def test_map_maker(pars, obs_list):
         geom=pars['geom'],
         offset_max='2 deg',
     )
+
     maps = maker.run(obs_list)
 
     counts = maps['counts']
@@ -70,3 +74,17 @@ def test_map_maker(pars, obs_list):
     background = maps['background']
     assert background.unit == ""
     assert_allclose(background.data.sum(), pars['background'], rtol=1e-5)
+
+    images = maker.make_images()
+
+    counts = images['counts']
+    assert counts.unit == ""
+    assert_allclose(counts.data.sum(), pars['counts'], rtol=1e-5)
+
+    exposure = images['exposure']
+    assert exposure.unit == "m2 s"
+    assert_allclose(exposure.data.sum(), pars['exposure_image'], rtol=1e-5)
+
+    background = images['background']
+    assert background.unit == ""
+    assert_allclose(background.data.sum(), pars['background'], rtol=1e-5)