Add .fill methods to SpectralCube and SkyMap

gammapy/cube/spectral_cube.py

@@ -19,8 +19,10 @@
 
 from ..spectrum import LogEnergyAxis, powerlaw
 from ..utils.energy import EnergyBounds
-from ..image import cube_to_image, solid_angle
+from ..image import cube_to_image, solid_angle, SkyMap
+from ..image.utils import _bin_events_in_cube 
 from ..utils.fits import table_to_fits_table
+from ..utils.wcs import get_wcs_ctype
 
 __all__ = ['SpectralCube']
 
@@ -49,6 +51,8 @@ class SpectralCube(object):
 
     Parameters
     ----------
+    name : str
+        Name of the spectral cube.
     data : `~astropy.units.Quantity`
         Data array (3-dim)
     wcs : `~astropy.wcs.WCS`
@@ -77,12 +81,12 @@ class SpectralCube(object):
     http://fermi.gsfc.nasa.gov/ssc/data/analysis/software/aux/gal_2yearp7v6_v0.fits
     """
 
-    def __init__(self, data, wcs, energy):
+    def __init__(self, name=None, data=None, wcs=None, energy=None):
         # TODO: check validity of inputs
+        self.name = name
         self.data = data
         self.wcs = wcs
-        self.header = wcs.to_header()
-
+
         # TODO: decide whether we want to use an EnergyAxis object or just use the array directly.
         self.energy = energy
         self.energy_axis = LogEnergyAxis(energy)
@@ -133,64 +137,92 @@ def read_hdu(cls, hdu_list):
         return cls(data, wcs, energy)
 
     @classmethod
-    def read(cls, filename):
+    def read(cls, filename, format='fermi'):
         """Read spectral cube from FITS file.
 
         Parameters
         ----------
         filename : str
             File name
+        format : {'fermi', 'fermi-counts'}
+            Fits file format.
 
         Returns
         -------
         spectral_cube : `SpectralCube`
             Spectral cube
         """
         data = fits.getdata(filename)
-        data = Quantity(data, '1 / (cm2 MeV s sr)')
         # 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 = Table.read(filename, 'ENERGIES')['Energy']
-        energy = Quantity(energy, 'MeV')
-
-        return cls(data, wcs, energy)
+        wcs = WCS(header).celestial
+        if format == 'fermi':
+            energy = Table.read(filename, 'ENERGIES')['Energy']
+            energy = Quantity(energy, 'MeV')
+            data = Quantity(data, '1 / (cm2 MeV s sr)')
+        elif format == 'fermi-counts':
+            energy = EnergyBounds.from_ebounds(fits.open(filename)['EBOUNDS'], unit='keV')
+            data = Quantity(data, 'counts')
+        else:
+            raise ValueError('Not a valid cube fits format')
+        return cls(data=data, wcs=wcs, energy=energy)
 
-    @classmethod
-    def read_counts(cls, filename):
-        """Read spectral cube of type "counts" from FITS file.
+    def fill(self, events, origin=0):
+        """
+        Fill spectral cube with events.
 
         Parameters
         ----------
-        filename : str
-            File name
+        events : `~astropy.table.Table`
+            Event list table
+        origin : {0, 1}
+            Pixel coordinate origin.
 
-        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)
-
+        self.data = _bin_events_in_cube(events, self.wcs, self.data.shape, self.energy)
+
     @classmethod
-    def empty_like(cls, skydata, energies=None):
+    def empty(cls, emin=0.5, emax=100, enbins=10, eunit='TeV', **kwargs):
         """
-        Create empty spectral cube from given skymap or spectral cube.
+        Create empty spectral cube with log equal energy binning from the scratch.
+        
+        Parameters
+        ----------
+        emin : float
+            Minimum energy.
+        emax : float
+            Maximum energy.
+        enbins : int
+            Number of energy bins.
+        eunit : str
+            Energy unit.
+        kwargs : dict
+            Keyword arguments passed to `~gammapy.image.SkyMap.empty` to create
+            the spatial part of the cube.
         """
+        refmap = SkyMap.empty(**kwargs)
+        energies = EnergyBounds.equal_log_spacing(emin, emax, enbins, eunit)
+        data = refmap.data * np.ones(len(energies)).reshape((-1, 1, 1))
+        return cls(data=data, wcs=refmap.wcs, energy=energies)
 
-
-        return cls(data, wcs, energy)
-
+    @classmethod
+    def empty_like(cls, refcube, fill=0):
+        """
+        Create an empty spectral cube with the same WCS and energy specification
+        as given spectral cube. 
+        
+        Parameters
+        ----------
+        refcube : `~gammapy.cube.SpectralCube`
+            Reference spectral cube.
+        fill : float, optional
+            Fill sky map with constant value. Default is 0.
+        """
+        wcs = refcube.wcs.copy()
+        data = fill * np.ones_like(refcube.data)
+        energies = refcube.energies.copy()
+        return cls(data=data, wcs=wcs, energy=energies)        
 
     def world2pix(self, lon, lat, energy, combine=False):
         """Convert world to pixel coordinates.
@@ -235,7 +267,7 @@ def pix2world(self, x, y, z):
         lon, lat, energy
         """
         origin = 0  # convention for gammapy
-        lon, lat, _ = self.wcs.wcs_pix2world(x, y, 0, origin)
+        lon, lat = self.wcs.wcs_pix2world(x, y, origin)
         energy = self.energy_axis.pix2world(z)
 
         lon = Quantity(lon, 'deg')
@@ -508,7 +540,7 @@ def writeto(self, filename, **kwargs):
 
     def __repr__(self):
         # Copied from `spectral-cube` package
-        s = "SpectralCube with shape={0}".format(self.data.shape)
+        s = "Spectral cube {} with shape={0}".format(self.name, self.data.shape)
         if self.data.unit is u.dimensionless_unscaled:
             s += ":\n"
         else:

gammapy/image/maps.py

@@ -13,7 +13,7 @@
 from astropy.extern import six
 
 from ..extern.bunch import Bunch
-from ..image.utils import make_header
+from ..image.utils import make_header, _bin_events_in_cube 
 from ..utils.wcs import get_wcs_ctype
 from ..utils.scripts import make_path
 
@@ -87,7 +87,7 @@ def read(cls, fobj, *args, **kwargs):
             name = header.get('EXTNAME')
         try:
             # Valitade unit string
-            unit = Unit(header['BUNIT']).to_string()
+            unit = Unit(header['BUNIT'], format='fits').to_string()
         except (KeyError, ValueError):
             unit = None
             log.warn('No valid units found for extension {}'.format(name))
@@ -178,6 +178,21 @@ def empty_like(cls, skymap, name=None, unit=None, fill=0, meta=None):
         return cls(name, data, wcs, unit, meta=wcs.to_header())
 
 
+    def fill(self, events, origin=0):
+        """
+        Fill sky map with events.
+
+        Parameters
+        ----------
+        events : `~astropy.table.Table`
+            Event list table
+        origin : {0, 1}
+            Pixel coordinate origin.
+
+        """
+        self.data = _bin_events_in_cube(events, self.wcs, self.data.shape, origin=origin).sum(axis=0)
+
+
     def write(self, filename, *args, **kwargs):
         """
         Write sky map to Fits file.
@@ -210,7 +225,9 @@ def coordinates(self, coord_type='world', origin=0, mode='center'):
         if mode == 'center':
             y, x = np.indices(self.data.shape)
         elif mode == 'edges':
-            raise NotImplementedError
+            shape = self.data.shape[0] + 1, self.data.shape[1] + 1 
+            y, x = np.indices(shape)
+            y, x = y - 0.5, x - 0.5
         else:
             raise ValueError('Invalid mode to compute coordinates.')
 
@@ -228,6 +245,15 @@ def coordinates(self, coord_type='world', origin=0, mode='center'):
                 raise ValueError("Not a valid coordinate type. Choose either"
                                  " 'world', 'pix' or 'skycoord'.")
 
+    def solid_angle(self):
+        """
+        Solid angle image
+        """
+        xsky, ysky = self.coordinates(mode='edges')
+        omega = -np.diff(xsky, axis=1)[1:, :] * np.diff(ysky, axis=0)[:, 1:]
+        return Quantity(omega, 'deg2').to('sr')
+
+
     def lookup(self, position, interpolation=None, origin=0):
         """
         Lookup value at given sky position.
@@ -371,11 +397,13 @@ def __repr__(self):
         String representation of the class.
         """
         info = "Name: {}\n".format(self.name)
-        info += "Data shape: {}\n".format(self.data.shape)
-        info += "Data type: {}\n".format(self.data.dtype)
-        info += "Data unit: {}\n".format(self.unit)
-        info += "Data mean: {:.3e}\n".format(np.nanmean(self.data))
-        info += "WCS type: {}\n".format(self.wcs.wcs.ctype)
+        if not self.data is None:
+            info += "Data shape: {}\n".format(self.data.shape)
+            info += "Data type: {}\n".format(self.data.dtype)
+            info += "Data unit: {}\n".format(self.unit)
+            info += "Data mean: {:.3e}\n".format(np.nanmean(self.data))
+        if not self.wcs is None:
+            info += "WCS type: {}\n".format(self.wcs.wcs.ctype)
         return info
 
     def __array__(self):

gammapy/image/measure.py

@@ -4,7 +4,7 @@
 from __future__ import absolute_import, division, print_function, unicode_literals
 import numpy as np
 from astropy.io import fits
-from ..image.utils import coordinates
+from .utils import coordinates
 
 __all__ = [
     'BoundingBox',

gammapy/image/tests/test_maps.py

@@ -6,6 +6,7 @@
 from astropy.io import fits
 
 from ..maps import SkyMap
+from ...data import DataStore
 from ...datasets import load_poisson_stats_image
 from ...utils.testing import requires_dependency, requires_data
 
@@ -67,3 +68,17 @@ def test_empty(self):
         empty = SkyMap.empty()
         assert empty.data.shape == (200, 200)
 
+
+    @requires_data('gammapy-extra')
+    def test_fill(self):
+        dirname = '$GAMMAPY_EXTRA/datasets/hess-crab4-hd-hap-prod2'
+        data_store = DataStore.from_dir(dirname)
+
+        events = data_store.load(obs_id=23523, filetype='events')
+
+        counts = SkyMap.empty(nxpix=200, nypix=200, xref=events.meta['RA_OBJ'],
+                              yref=events.meta['DEC_OBJ'], dtype='int', 
+                              coordsys='CEL')
+        counts.fill(events)
+        assert counts.data.sum() == 1233
+        assert counts.data.shape == (200, 200)

gammapy/image/tests/test_utils.py

@@ -23,8 +23,7 @@
     block_reduce_hdu,
     wcs_histogram2d,
     lon_lat_rectangle_mask,
-    SkyMap
-    bin_events_in_cube,
+    SkyMap,
 )
 
 
@@ -224,37 +223,12 @@ def test_bin_events_in_cube():
 
     events = data_store.load(obs_id=23523, filetype='events')
 
-    # Define WCS reference header
-    refheader = fits.Header()
-    refheader['WCSAXES'] = 3
-    refheader['NAXIS'] = 3
-    refheader['CRPIX1'] = 100.5
-    refheader['CRPIX2'] = 100.5
-    refheader['CRPIX3'] = 1.0
-    refheader['CDELT1'] = -0.02
-    refheader['CDELT2'] = 0.02
-
-    # shouldn't matter, but must contain sufficient number of digits,
-    # so that CDELT1 and CDELT2 are not truncated, when wcs.to_header() is called
-    # seems to be a bug...
-    refheader['CDELT3'] = 2.02  
-
-    refheader['CTYPE1'] = 'RA---CAR'
-    refheader['CTYPE2'] = 'DEC--CAR'
-    refheader['CTYPE3'] = 'log_Energy'  # shouldn't matter
-    refheader['CUNIT1'] = 'deg'
-    refheader['CUNIT2'] = 'deg'
-    refheader['CRVAL1'] = events.meta['RA_OBJ']
-    refheader['CRVAL2'] = events.meta['DEC_OBJ']
-    refheader['CRVAL3'] = 10.0  # shouldn't matter
-
-    energies = EnergyBounds.equal_log_spacing(0.5, 80, 8, 'TeV')
-    data = Quantity(np.zeros((len(energies), 200, 200)))
-    wcs = WCS(refheader)
-    refcube = SpectralCube(data, wcs, energy=energies)
-
-    # Counts cube
-    counts_hdu = bin_events_in_cube(events, refcube, energies)
-    counts = SpectralCube(Quantity(counts_hdu.data, 'count'), wcs, energies)
-
-    assert counts.data.sum().value == 1233
+    counts = SpectralCube.empty(emin=0.5, emax=80, enbins=8, eunit='TeV',
+                                nxpix=200, nypix=200, xref=events.meta['RA_OBJ'],
+                                yref=events.meta['DEC_OBJ'], dtype='int', 
+                                coordsys='CEL')
+    counts.fill(events)
+    assert counts.data.sum() == 1233
+
+
+