Merge pull request #538 from OlgaVorokh/solid_angle

Delete solid_angle function and update all usages of it

gammapy/cube/core.py

@@ -19,10 +19,9 @@
 
 from ..spectrum import LogEnergyAxis, powerlaw
 from ..utils.energy import EnergyBounds
-from ..image import cube_to_image, solid_angle, SkyMap
+from ..image import 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__ = ['SkyCube']
 
@@ -318,10 +317,13 @@ def to_sherpa_data3d(self):
     def solid_angle_image(self):
         """Solid angle image in steradian (`~astropy.units.Quantity`)"""
         cube_hdu = fits.ImageHDU(self.data, self.wcs.to_header())
+
+        from .utils import cube_to_image
         image_hdu = cube_to_image(cube_hdu)
         image_hdu.header['WCSAXES'] = 2
 
-        return solid_angle(image_hdu).to('sr')
+        sky_map = SkyMap.read(image_hdu)
+        return sky_map.solid_angle()
 
     def flux(self, lon, lat, energy):
         """Differential flux.

gammapy/cube/tests/test_utils.py

@@ -0,0 +1,17 @@
+from __future__ import absolute_import, division, print_function, unicode_literals
+from numpy.testing import assert_allclose
+
+from ..utils import cube_to_image
+from ...image import images_to_cube, SkyMap
+
+
+def test_cube_to_image():
+    layer = SkyMap.empty(nxpix=101, nypix=101, fill=1.).to_image_hdu()
+    hdu_list = [layer, layer, layer, layer]
+    cube = images_to_cube(hdu_list)
+    case1 = cube_to_image(cube)
+    case2 = cube_to_image(cube, slicepos=1)
+    # Check that layers are summed if no layer is specified (case1),
+    # or only a specified layer is extracted (case2)
+    assert_allclose(case1.data, 4 * layer.data)
+    assert_allclose(case2.data, layer.data)

gammapy/cube/utils.py

@@ -6,12 +6,19 @@
 from __future__ import absolute_import, division, print_function, unicode_literals
 import numpy as np
 
+from astropy.io.fits import ImageHDU
 from astropy.units import Quantity
 from astropy.coordinates import Angle
 
 from .core import SkyCube
+from ..image.maps import SkyMap
 
-__all__ = ['compute_npred_cube', 'convolve_cube']
+
+__all__ = ['compute_npred_cube',
+           'convolve_cube',
+           'cube_to_image',
+           'cube_to_spec',
+           ]
 
 
 def compute_npred_cube(flux_cube, exposure_cube, energy_bins,
@@ -96,3 +103,62 @@ def convolve_cube(cube, psf, offset_max):
                              energy=cube.energy)
     return convolved_cube
 
+
+def cube_to_image(cube, slicepos=None):
+    """Slice or project 3-dim cube into a 2-dim image.
+
+    Parameters
+    ----------
+    cube : `~astropy.io.fits.ImageHDU`
+        3-dim FITS cube
+    slicepos : int or None, optional
+        Slice position (None means to sum along the spectral axis)
+
+    Returns
+    -------
+    image : `~astropy.io.fits.ImageHDU`
+        2-dim FITS image
+    """
+    header = cube.header.copy()
+    header['NAXIS'] = 2
+
+    for key in ['NAXIS3', 'CRVAL3', 'CDELT3', 'CTYPE3', 'CRPIX3', 'CUNIT3']:
+        if key in header:
+            del header[key]
+
+    if slicepos is None:
+        data = cube.data.sum(0)
+    else:
+        data = cube.data[slicepos]
+    return ImageHDU(data, header)
+
+
+def cube_to_spec(cube, mask, weighting='none'):
+    """Integrate spatial dimensions of a FITS cube to give a spectrum.
+
+    TODO: give formulas.
+
+    Parameters
+    ----------
+    cube : `~astropy.io.fits.ImageHDU`
+        3-dim FITS cube
+    mask : numpy.array
+        2-dim mask array.
+    weighting : {'none', 'solid_angle'}, optional
+        Weighting factor to use.
+
+    Returns
+    -------
+    spectrum : numpy.array
+        Summed spectrum of pixels in the mask.
+    """
+
+    # TODO: clean up API and implementation and add test
+    value = cube.dat
+    sky_map = SkyMap.read(cube)
+    A = sky_map.solid_angle()
+    # Note that this is the correct way to get an average flux:
+
+    spec = (value * A).sum(-1).sum(-1)
+    return spec
+

gammapy/image/tests/test_utils.py

@@ -10,15 +10,13 @@
 from ...datasets import FermiGalacticCenter
 from ...data import DataStore
 from ...utils.energy import EnergyBounds
-from ...cube import SkyCube
+from ...cube import SkyCube, cube_to_image
 from ...image import (
     binary_disk,
     binary_ring,
     make_header,
     contains,
-    solid_angle,
     images_to_cube,
-    cube_to_image,
     block_reduce_hdu,
     wcs_histogram2d,
     lon_lat_rectangle_mask,
@@ -73,15 +71,6 @@ def test_contains(self):
         inside = contains(self.image, x, y)
         assert_equal(inside, np.ones((3, 2), dtype=bool))
 
-    # TODO: this works on my machine, but fails for unknown reasons
-    # with an IndexError with the `numpy` used here:
-    # https://travis-ci.org/gammapy/gammapy/jobs/26836201#L1123
-    @pytest.mark.xfail
-    def test_image_area(self):
-        actual = solid_angle(self.image)
-        expected = 99.61946869
-        assert_allclose(actual, expected)
-
 
 @pytest.mark.xfail
 def test_process_image_pixels():
@@ -161,18 +150,6 @@ def test_ref_pixel():
     assert_allclose(footprint[0], footprint_1[0])
 
 
-def test_cube_to_image():
-    layer = SkyMap.empty(nxpix=101, nypix=101, fill=1.).to_image_hdu()
-    hdu_list = [layer, layer, layer, layer]
-    cube = images_to_cube(hdu_list)
-    case1 = cube_to_image(cube)
-    case2 = cube_to_image(cube, slicepos=1)
-    # Check that layers are summed if no layer is specified (case1),
-    # or only a specified layer is extracted (case2)
-    assert_allclose(case1.data, 4 * layer.data)
-    assert_allclose(case2.data, layer.data)
-
-
 def test_wcs_histogram2d():
     # A simple test case that can by checked by hand:
     header = make_header(nxpix=2, nypix=1, binsz=10, xref=0, yref=0, proj='CAR')

gammapy/image/utils.py

@@ -24,8 +24,6 @@
     'binary_ring',
     'block_reduce_hdu',
     'contains',
-    'cube_to_image',
-    'cube_to_spec',
     'crop_image',
     'dict_to_hdulist',
     'disk_correlate',
@@ -40,7 +38,6 @@
     'process_image_pixels',
     'ring_correlate',
     'shape_2N',
-    'solid_angle',
     'threshold',
     'upsample_2N',
     'wcs_histogram2d',
@@ -765,35 +762,6 @@ def binary_opening_circle(input, radius):
     return binary_opening(input, structure)
 
 
-def solid_angle(image):
-    """Compute the solid angle of each pixel.
-
-    This will only give correct results for CAR maps!
-
-    Parameters
-    ----------
-    image : `~astropy.io.fits.ImageHDU`
-        Input image
-
-    Returns
-    -------
-    area_image : `~astropy.units.Quantity`
-        Solid angle image (matching the input image) in steradians.
-    """
-    # Area of one pixel at the equator
-    cdelt0 = image.header['CDELT1']
-    cdelt1 = image.header['CDELT2']
-    equator_area = Quantity(abs(cdelt0 * cdelt1), 'deg2')
-
-    # Compute image with fraction of pixel area at equator
-    glat = coordinates(image)[1]
-    area_fraction = np.cos(np.radians(glat))
-
-    result = area_fraction * equator_area.to('sr')
-
-    return result
-
-
 def make_header(nxpix=100, nypix=100, binsz=0.1, xref=0, yref=0,
                 proj='CAR', coordsys='GAL',
                 xrefpix=None, yrefpix=None):
@@ -884,63 +852,6 @@ def crop_image(image, bounding_box):
     return fits.ImageHDU(data=data, header=header)
 
 
-def cube_to_image(cube, slicepos=None):
-    """Slice or project 3-dim cube into a 2-dim image.
-
-    Parameters
-    ----------
-    cube : `~astropy.io.fits.ImageHDU`
-        3-dim FITS cube
-    slicepos : int or None, optional
-        Slice position (None means to sum along the spectral axis)
-
-    Returns
-    -------
-    image : `~astropy.io.fits.ImageHDU`
-        2-dim FITS image
-    """
-    from astropy.io.fits import ImageHDU
-    header = cube.header.copy()
-    header['NAXIS'] = 2
-
-    for key in ['NAXIS3', 'CRVAL3', 'CDELT3', 'CTYPE3', 'CRPIX3', 'CUNIT3']:
-        if key in header:
-            del header[key]
-
-    if slicepos is None:
-        data = cube.data.sum(0)
-    else:
-        data = cube.data[slicepos]
-    return ImageHDU(data, header)
-
-
-def cube_to_spec(cube, mask, weighting='none'):
-    """Integrate spatial dimensions of a FITS cube to give a spectrum.
-
-    TODO: give formulas.
-
-    Parameters
-    ----------
-    cube : `~astropy.io.fits.ImageHDU`
-        3-dim FITS cube
-    mask : numpy.array
-        2-dim mask array.
-    weighting : {'none', 'solid_angle'}, optional
-        Weighting factor to use.
-
-    Returns
-    -------
-    spectrum : numpy.array
-        Summed spectrum of pixels in the mask.
-    """
-    value = cube.dat
-    A = solid_angle(cube)
-    # Note that this is the correct way to get an average flux:
-
-    spec = (value * A).sum(-1).sum(-1)
-    return spec
-
-
 def contains(image, x, y, world=True):
     """Check if given pixel or world positions are in an image.