Merge 6e1b0af into f30420f

gammapy/image/profile.py

@@ -2,8 +2,12 @@
 """Tools to create profiles (i.e. 1D "slices" from 2D images)"""
 from __future__ import print_function, division
 import numpy as np
+from astropy.table import Table
+from astropy.units import Quantity
+from .utils import coordinates
 
-__all__ = ['compute_binning', 'FluxProfile']
+__all__ = ['compute_binning', 'FluxProfile',
+           'image_lat_profile', 'image_lon_profile']
 
 
 def compute_binning(data, n_bins, method='equal width', eps=1e-10):
@@ -186,3 +190,199 @@ def save(self, filename):
         ----------
         """
         raise NotImplementedError
+
+
+def image_lat_profile(image, lats, lons, binsz, counts=None,
+                      mask_array=None, invert=False, errors=False,
+                      standard_error=0.1):
+    """Creates a latitude profile from input flux image HDU.
+
+    Parameters
+    ----------
+    image : `~astropy.io.fits.ImageHDU`
+        Image HDU object to produce GLAT profile.
+    lats : array_like
+        Specified as [GLAT_min, GLAT_max], with GLAT_min and GLAT_max
+        as floats. A 1x2 array specifying the maximum and minimum latitudes to
+        include in the region of the image for which the profile is formed, 
+        which should be within the spatial bounds of the image.
+    lons : array_like
+        Specified as [GLON_min, GLON_max], with GLON_min and GLON_max
+        as floats. A 1x2 array specifying the maximum and minimum longitudes to
+        include in the region of the image for which the profile is formed, 
+        which should be within the spatial bounds of the image.
+    binsz : float
+        Latitude bin size of the resulting latitude profile. This should be
+        no less than 5 times the pixel resolution.
+    counts : `~astropy.io.fits.ImageHDU`
+        Counts image to allow Poisson errors to be calculated. If not provided,
+        a standard_error should be provided, or zero errors will be returned.
+        (Optional).
+    mask_array : array_like
+        2D mask array, matching spatial dimensions of input image. (Optional).
+    invert : bool
+        If True, inverts the mask array. Default False.
+    errors : bool
+        If True, computes errors, if possible, according to provided inputs.
+        If False (default), returns all errors as zero.
+    standard_error : float
+        If counts image is not provided, but error values required, this
+        specifies a standard fractional error to be applied to values.
+        Default = 0.1.
+
+    Returns
+    -------
+    table : `~astropy.table.Table`
+        Galactic latitude profile as table, with latitude bin boundaries,
+        profile values and errors.
+    """
+    res_bound = 5 * np.abs(image.header['CDELT2'])
+    # Assert that the called binsz is at least 5 times
+    # GLAT image resolution to avoid binning bug
+    if binsz < res_bound:
+        raise ValueError
+
+    lon, lat = coordinates(image)
+    mask_init = (lats[0] <= lat) & (lat < lats[1])
+    mask = mask_init & (lons[0] <= lon) & (lon < lons[1])
+    if mask_array != None:
+        if invert == True:
+            mask_array = np.invert(mask_array)
+        mask = mask & mask_array
+
+    # Need to preserve shape here so use multiply
+    cut_image = image.data * mask
+    if counts != None:
+        cut_counts = counts.data * mask
+    values = []
+    count_vals = []
+    bins = np.arange((lats[1] - lats[0]) / binsz)
+
+    glats_min = lats[0] + bins[:-1] * binsz
+    glats_max = lats[0] + bins[1:] * binsz
+
+    # Table is required here to avoid rounding problems with floats
+    bounds = Table([glats_min, glats_max], names=('GLAT_MIN', 'GLAT_MAX'))
+    for bin in bins[:-1]:
+        lat_mask = (bounds['GLAT_MIN'][bin] <= lat) & (lat < bounds['GLAT_MAX'][bin])
+        lat_band = cut_image[lat_mask]
+        values.append(lat_band.sum())
+        if counts != None:
+            count_band = cut_counts[lat_mask]
+            count_vals.append(count_band.sum())
+        else:
+            count_vals.append(0)
+
+    if errors == True:
+        if counts != None:
+            rel_errors = 1. / np.sqrt(count_vals)
+            error_vals = values * rel_errors
+        else:
+            error_vals = values * (np.ones(len(values)) * standard_error)
+    else:
+        error_vals = np.zeros_like(values)
+
+    table = Table([Quantity(glats_min, 'deg'),
+                   Quantity(glats_max, 'deg'),
+                   values,
+                   error_vals],
+                  names=('GLAT_MIN', 'GLAT_MAX', 'BIN_VALUE', 'BIN_ERR'))
+    return table
+
+
+def image_lon_profile(image, lats, lons, binsz, counts=None,
+                      mask_array=None, invert=False, errors=False,
+                      standard_error=0.1):
+    """Creates a longitude profile from input flux image HDU.
+
+    Parameters
+    ----------
+    image : `~astropy.io.fits.ImageHDU`
+        Image HDU object to produce GLON profile.
+    lats : array_like
+        Specified as [GLAT_min, GLAT_max], with GLAT_min and GLAT_max
+        as floats. A 1x2 array specifying the maximum and minimum latitudes to
+        include in the region of the image for which the profile is formed, 
+        which should be within the spatial bounds of the image.
+    lons : array_like
+        Specified as [GLON_min, GLON_max], with GLON_min and GLON_max
+        as floats. A 1x2 array specifying the maximum and minimum longitudes to
+        include in the region of the image for which the profile is formed, 
+        which should be within the spatial bounds of the image.
+    binsz : float
+        Latitude bin size of the resulting longitude profile. This should be
+        no less than 5 times the pixel resolution.
+    counts : `~astropy.io.fits.ImageHDU`
+        Counts image to allow Poisson errors to be calculated. If not provided,
+        a standard_error should be provided, or zero errors will be returned.
+        (Optional).
+    mask_array : array_like
+        2D mask array, matching spatial dimensions of input image. (Optional).
+    invert : bool
+        If True, inverts the mask array. Default False.
+    errors : bool
+        If True, computes errors, if possible, according to provided inputs.
+        If False (default), returns all errors as zero.
+    standard_error : float
+        If counts image is not provided, but error values required, this
+        specifies a standard fractional error to be applied to values.
+        Default = 0.1.
+
+    Returns
+    -------
+    table : `~astropy.table.Table`
+        Galactic longitude profile as table, with longitude bin boundaries,
+        profile values and errors.
+    """
+    res_bound = 5 * np.abs(image.header['CDELT1'])
+    # Assert that the called binsz is at least 5 times
+    # GLON image resolution to avoid binning bug
+    if binsz < res_bound:
+        raise ValueError
+
+    lon, lat = coordinates(image)
+    mask_init = (lats[0] <= lat) & (lat < lats[1])
+    mask = mask_init & (lons[0] <= lon) & (lon < lons[1])
+    if mask_array != None:
+        if invert == True:
+            mask_array = np.invert(mask_array)
+        mask = mask & mask_array
+
+    # Need to preserve shape here so use multiply
+    cut_image = image.data * mask
+    if counts != None:
+        cut_counts = counts.data * mask
+    values = []
+    count_vals = []
+    bins = np.arange((lons[1] - lons[0]) / binsz)
+
+    glons_min = lons[0] + bins[:-1] * binsz
+    glons_max = lons[0] + bins[1:] * binsz
+
+    # Table is required here to avoid rounding problems with floats
+    bounds = Table([glons_min, glons_max], names=('GLON_MIN', 'GLON_MAX'))
+    for bin in bins[:-1]:
+        lon_mask = (bounds['GLON_MIN'][bin] <= lon) & (lon < bounds['GLON_MAX'][bin])
+        lon_band = cut_image[lon_mask]
+        values.append(lon_band.sum())
+        if counts != None:
+            count_band = cut_counts[lon_mask]
+            count_vals.append(count_band.sum())
+        else:
+            count_vals.append(0)
+
+    if errors == True:
+        if counts != None:
+            rel_errors = 1. / np.sqrt(count_vals)
+            error_vals = values * rel_errors
+        else:
+            error_vals = values * (np.ones(len(values)) * standard_error)
+    else:
+        error_vals = np.zeros_like(values)
+
+    table = Table([Quantity(glons_min, 'deg'),
+                   Quantity(glons_max, 'deg'),
+                   values,
+                   error_vals],
+                  names=('GLON_MIN', 'GLON_MAX', 'BIN_VALUE', 'BIN_ERR'))
+    return table

gammapy/image/tests/test_profile.py

@@ -1,9 +1,14 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 from __future__ import print_function, division
+import numpy as np
 from numpy.testing import assert_allclose
 from astropy.tests.helper import pytest
+from astropy.units import Quantity
+from ...datasets import FermiGalacticCenter
+from ..utils import coordinates
 from .. import profile
 
+
 try:
     import pandas
     HAS_PANDAS = True
@@ -20,3 +25,91 @@ def test_compute_binning():
     bin_edges = profile.compute_binning(data, n_bins=3, method='equal entries')
     # TODO: create test-cases that have been verified by hand here!
     assert_allclose(bin_edges, [1,  2,  2.66666667,  4])
+
+
+def test_image_lat_profile():
+    """Tests GLAT profile with image of 1s of known size and shape."""
+    image = FermiGalacticCenter.counts()
+
+    lons, lats = coordinates(image)
+    image.data = np.ones_like(image.data)
+
+    counts = FermiGalacticCenter.counts()
+    counts.data = np.ones_like(counts.data)
+
+    mask = np.zeros_like(image.data)
+    # Select Full Image
+    lat = [lats.min(), lats.max()]
+    lon = [lons.min(), lons.max()]
+    # Pick minimum valid binning
+    binsz = 0.5
+    mask_array = np.zeros_like(image.data)
+    # Test output
+    lat_profile1 = profile.image_lat_profile(image, lat, lon, binsz,
+                                             errors=True)
+    # atol 0.1 is sufficient to check if correct number of pixels are included
+    assert_allclose(lat_profile1['BIN_VALUE'].data.astype(float),
+                    2000 * np.ones(39), rtol=1, atol=0.1)
+    assert_allclose(lat_profile1['BIN_ERR'].data,
+                    0.1 * lat_profile1['BIN_VALUE'].data)
+
+    lat_profile2 = profile.image_lat_profile(image, lat, lon, binsz,
+                                             counts, errors=True)
+    # atol 0.1 is sufficient to check if correct number of pixels are included
+    assert_allclose(lat_profile2['BIN_ERR'].data,
+                    44.721359549995796 * np.ones(39), rtol=1, atol=0.1)
+
+    lat_profile3 = profile.image_lat_profile(image, lat, lon, binsz, counts,
+                                    mask_array, invert=False, errors=True)
+
+    assert_allclose(lat_profile3['BIN_VALUE'].data, np.zeros(39))
+
+    lat_profile4 = profile.image_lat_profile(image, lat, lon, binsz, counts,
+                                    mask_array, invert=True, errors=True)
+    # Result should be the same as in case 2 as zero mask_array is inverted
+    assert_allclose(lat_profile4['BIN_VALUE'].data,
+                    lat_profile2['BIN_VALUE'].data)
+
+
+def test_image_lon_profile():
+    """Tests GLON profile with image of 1s of known size and shape."""
+    image = FermiGalacticCenter.counts()
+
+    lons, lats = coordinates(image)
+    image.data = np.ones_like(image.data)
+
+    counts = FermiGalacticCenter.counts()
+    counts.data = np.ones_like(counts.data)
+
+    mask = np.zeros_like(image.data)
+    # Select Full Image
+    lat = [lats.min(), lats.max()]
+    lon = [lons.min(), lons.max()]
+    # Pick minimum valid binning
+    binsz = 0.5
+    mask_array = np.zeros_like(image.data)
+    # Test output
+    lon_profile1 = profile.image_lon_profile(image, lat, lon, binsz,
+                                             errors=True)
+    # atol 0.1 is sufficient to check if correct number of pixels are included
+    assert_allclose(lon_profile1['BIN_VALUE'].data.astype(float),
+                    1000 * np.ones(79), rtol=1, atol=0.1)
+    assert_allclose(lon_profile1['BIN_ERR'].data,
+                    0.1 * lon_profile1['BIN_VALUE'].data)
+
+    lon_profile2 = profile.image_lon_profile(image, lat, lon, binsz,
+                                             counts, errors=True)
+    # atol 0.1 is sufficient to check if correct number of pixels are included
+    assert_allclose(lon_profile2['BIN_ERR'].data,
+                    31.622776601683793 * np.ones(79), rtol=1, atol=0.1)
+
+    lon_profile3 = profile.image_lon_profile(image, lat, lon, binsz, counts,
+                                    mask_array, invert=False, errors=True)
+
+    assert_allclose(lon_profile3['BIN_VALUE'].data, np.zeros(79))
+
+    lon_profile4 = profile.image_lon_profile(image, lat, lon, binsz, counts,
+                                    mask_array, invert=True, errors=True)
+    # Result should be the same as in case 2 as zero mask_array is inverted
+    assert_allclose(lon_profile4['BIN_VALUE'].data,
+                    lon_profile2['BIN_VALUE'].data)