Add gammapy.region and reflected region computation

docs/index.rst

@@ -75,6 +75,7 @@ The Gammapy toolbox
   irf/index
   morphology/index
   obs/index
+  region/index
   spectrum/index
   stats/index
   time/index

docs/region/index.rst

@@ -0,0 +1,58 @@
+.. include:: ../references.txt
+
+.. _region:
+
+**************************
+Regions (`gammapy.region`)
+**************************
+
+.. currentmodule:: gammapy.region
+
+Introduction
+============
+
+`gammapy.region` contains classes and functions for region handling.
+
+
+Getting Started
+===============
+
+Creating regions
+----------------
+
+.. code-block:: python
+
+    >>> from gammapy.region import SkyCircleRegion
+    >>> from astropy.coordinates import Angle, SkyCoord
+    >>> pos = SkyCoord(80.2, 23.5, unit='deg', frame='icrs')
+    >>> radius = Angle(0.4, 'deg')
+    >>> region = SkyCircleRegion(pos=pos, radius=radius)
+
+Containment
+-----------
+
+Masks
+-----
+
+Rotate
+------
+
+Read / write
+------------
+
+Reflected Regions
+-----------------
+
+Details on the reflected regions method can be found in [Berge2007]_
+
+The following example illustrates how to create reflected regions
+for a given ImageHDU
+
+.. plot:: region/make_reflected_regions.py
+   :include-source:
+
+Reference/API
+=============
+
+.. automodapi:: gammapy.region
+    :no-inheritance-diagram:

docs/region/make_reflected_regions.py

@@ -0,0 +1,28 @@
+from astropy.coordinates import SkyCoord, Angle
+from astropy.wcs import WCS
+
+from gammapy.region import find_reflected_regions, SkyCircleRegion
+from gammapy.image import ExclusionMask, make_empty_image
+from gammapy.utils.testing import requires_data
+
+hdu = make_empty_image(nxpix=801, nypix=701, binsz=0.01,
+                       coordsys='CEL', xref=83.2, yref=22.7)
+mask = ExclusionMask.create_random(hdu, n=8, min_rad=30, max_rad=80)
+pos = SkyCoord(80.2, 23.5, unit='deg', frame='icrs')
+radius = Angle(0.4, 'deg')
+test_region = SkyCircleRegion(pos=pos, radius=radius)
+center = SkyCoord(82.8, 22.5, unit='deg', frame='icrs')
+regions = find_reflected_regions(test_region, center, mask)
+
+import matplotlib.pyplot as plt
+wcs = WCS(hdu.header)
+fig = plt.figure(figsize=(8, 5), dpi=80)
+ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection=wcs)
+mask.plot(ax)
+for reg in regions:
+    patch = reg.plot(ax)
+    patch.set_facecolor('red')
+patch2 =test_region.plot(ax)
+patch2.set_facecolor('blue')
+
+plt.show()

examples/test_regions.py

@@ -0,0 +1,10 @@
+from astropy.coordinates import Angle, SkyCoord
+from gammapy.region import SkyCircleRegion
+
+pos1 = SkyCoord(2, 1, unit='deg', frame='galactic')
+radius = Angle(1, 'deg')
+reg = SkyCircleRegion(pos=pos1, radius=radius)
+print(reg)
+
+
+import IPython; IPython.embed()

gammapy/data/counts_spectrum.py

@@ -10,6 +10,7 @@
 
 __all__ = ['CountsSpectrum']
 
+
 class CountsSpectrum(object):
     """Counts spectrum dataset
 
@@ -27,7 +28,7 @@ class CountsSpectrum(object):
     --------
 
     .. code-block:: python
-    
+
         from gammapy.spectrum import Energy, EnergyBounds
         from gammapy.data import CountsSpectrum
         ebounds = EnergyBounds.equal_log_spacing(1,10,10,'TeV') 
@@ -37,24 +38,24 @@ class CountsSpectrum(object):
     """
 
     def __init__(self, counts, energy, livetime=None, backscal=1):
-        
+
         if np.asarray(counts).size != energy.nbins:
             raise ValueError("Dimension of {0} and {1} do not"
                              " match".format(counts, energy))
-        
+
         if not isinstance(energy, Energy):
             raise ValueError("energy must be an Energy object")
-        
+
         self.counts = np.asarray(counts)
 
-        if isinstance (energy, EnergyBounds):
+        if isinstance(energy, EnergyBounds):
             self.energy = energy.log_centers
         else:
             self.energy = energy
 
         self._livetime = livetime
         self._backscal = backscal
-        self.channels = np.arange(1,self.energy.nbins+1,1)
+        self.channels = np.arange(1, self.energy.nbins + 1, 1)
 
     @property
     def entries(self):
@@ -93,7 +94,7 @@ def from_fits(cls, hdu):
     def from_eventlist(cls, event_list, bins):
         """Create CountsSpectrum from fits 'EVENTS' extension (`CountsSpectrum`).
 
-        Subsets of the event list should be choosen via the approriate methods
+        Subsets of the event list should be chosen via the appropriate methods
         in `~gammapy.data.EventList`.
 
         Parameters
@@ -104,15 +105,15 @@ def from_eventlist(cls, event_list, bins):
             Energy bin edges
         """
 
-        if isinstance (event_list,fits.BinTableHDU):
+        if isinstance(event_list, fits.BinTableHDU):
             event_list = EventList.read(event_list)
         elif isinstance(event_list, EventListDataset):
             event_list = event_list.event_list
         elif isinstance(event_list, str):
-            event_list = EventList.read(event_list, hdu='EVENTS') 
+            event_list = EventList.read(event_list, hdu='EVENTS')
 
         energy = Energy(event_list.energy).to(bins.unit)
-        val, dummy = np.histogram(energy,bins.value)
+        val, dummy = np.histogram(energy, bins.value)
         livetime = event_list.observation_live_time_duration
 
         return cls(val, bins, livetime)
@@ -122,13 +123,26 @@ def write(self, filename, bkg=None, corr=None, rmf=None, arf=None,
         """Write PHA to FITS file.
 
         Calls `~astropy.io.fits.HDUList.writeto`, forwarding all arguments.
+
+        Parameters
+        ----------
+        arf : str
+            ARF file to be written into FITS header
+        rmf : str
+            RMF file to be written into FITS header
+        corr : str
+            CORR file to be written into FITS header
+        bkg : str
+            BKG file to be written into FITS header
+        filename : str
+            File to be written
         """
         self.to_fits(bkg=bkg, corr=corr, rmf=rmf, arf=arf).writeto(
             filename, *args, **kwargs)
 
-    def to_fits(self, bkg=None, corr=None, rmf=None, arf=None ):
+    def to_fits(self, bkg=None, corr=None, rmf=None, arf=None):
         """Output OGIP pha file
-        
+
         Returns
         -------
         pha : `~astropy.io.fits.BinTableHDU`
@@ -141,16 +155,16 @@ def to_fits(self, bkg=None, corr=None, rmf=None, arf=None ):
         ogip_92_007_summary.html
         """
         col1 = fits.Column(name="CHANNEL", array=self.channels,
-                           format = 'I')
+                           format='I')
         col2 = fits.Column(name="COUNTS", array=self.counts,
                            format='J', unit='count')
-        
+
         cols = fits.ColDefs([col1, col2])
 
         hdu = fits.BinTableHDU.from_columns(cols)
         header = hdu.header
 
-        #TODO: option to store meta info in the class
+        # TODO: option to store meta info in the class
         header['EXTNAME'] = 'SPECTRUM', 'name of this binary table extension'
         header['TELESCOP'] = 'DUMMY', 'Telescope (mission) name'
         header['INSTRUME'] = 'DUMMY', 'Instrument name'
@@ -161,7 +175,7 @@ def to_fits(self, bkg=None, corr=None, rmf=None, arf=None ):
         header['CORRFILE'] = corr, 'Correlation FITS file'
         header['RESPFILE'] = rmf, 'Redistribution matrix file (RMF)'
         header['ANCRFILE'] = arf, 'Ancillary response file (ARF)'
-        
+
         header['HDUCLASS'] = 'OGIP', 'Format conforms to OGIP/GSFC spectral standards'
         header['HDUCLAS1'] = 'SPECTRUM', 'Extension contains a spectrum'
         header['HDUVERS '] = '1.2.1', 'Version number of the format'
@@ -182,18 +196,18 @@ def to_fits(self, bkg=None, corr=None, rmf=None, arf=None ):
         header['STAT_ERR'] = 0, 'No statisitcal error was specified'
         header['SYS_ERR'] = 0, 'No systematic error was specified'
         header['GROUPING'] = 0, 'No grouping data has been specified'
-        
+
         header['QUALITY '] = 0, 'No data quality information specified'
-        
+
         header['AREASCAL'] = 1., 'Nominal effective area'
         header['BACKSCAL'] = self.backscal, 'Background scale factor'
         header['CORRSCAL'] = 0., 'Correlation scale factor'
-        
+
         header['FILENAME'] = 'several', 'Spectrum was produced from more than one file'
         header['ORIGIN'] = 'dummy', 'origin of fits file'
         header['DATE'] = datetime.datetime.today().strftime('%Y-%m-%d'), 'FITS file creation date (yyyy-mm-dd)'
         header['PHAVERSN'] = '1992a', 'OGIP memo number for file format'
-        
+
         return hdu
 
     def plot(self, ax=None, filename=None, **kwargs):
@@ -204,12 +218,11 @@ def plot(self, ax=None, filename=None, **kwargs):
 
         ax = plt.gca() if ax is None else ax
 
-        plt.plot(self.energy, self.counts, 'o' , **kwargs)
+        plt.plot(self.energy, self.counts, 'o', **kwargs)
         plt.xlabel('Energy [{0}]'.format(self.energy.unit))
         plt.ylabel('Counts')
         if filename is not None:
             plt.savefig(filename)
             log.info('Wrote {0}'.format(filename))
 
         return ax
-

gammapy/data/event_list.py

@@ -8,6 +8,7 @@
 from astropy.time import Time
 from astropy.coordinates import SkyCoord, Angle, AltAz
 from astropy.table import Table
+
 from ..extern.pathlib import Path
 from ..image import wcs_histogram2d
 from ..data import GoodTimeIntervals, TelescopeArray
@@ -283,6 +284,51 @@ def select_sky_box(self, lon_lim, lat_lim, frame='icrs'):
         from ..catalog import select_sky_box
         return select_sky_box(self, lon_lim, lat_lim, frame)
 
+    def select_circular_region(self, region):
+        """Select events in circular regions
+
+        TODO: Extend to support generic regions
+
+        Parameters
+        ----------
+        region : `~gammapy.region.SkyRegionList`, `~gammapy.region.SkyCircleRegion`
+            (List of) sky region(s)
+
+        Returns
+        -------
+        event_list : `EventList`
+            Copy of event list with selection applied.
+        """
+
+        if not isinstance(region, list):
+            region = list([region])
+        mask = self.filter_circular_region(region)
+        return self[mask]
+
+    def filter_circular_region(self, region):
+        """Create selection mask for event in given circular regions
+
+        TODO: Extend to support generic regions
+
+        Parameters
+        ----------
+        region : `~gammapy.region.SkyRegionList`
+            List of sky regions
+
+        Returns
+        -------
+        index_array : `np.array`
+            Index array of seleced events
+        """
+
+        position = self.radec
+        mask = np.array([], dtype=int)
+        for reg in region:
+            separation = reg.pos.separation(position)
+            temp = np.where(separation < reg.radius)[0]
+            mask = np.union1d(mask, temp)
+        return mask
+
     def fill_counts_image(self, image):
         """Fill events in counts image.
 

gammapy/data/tests/test_event_list.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
 from numpy.testing import assert_allclose
+from astropy.coordinates import Angle, SkyCoord
 from astropy.tests.helper import pytest
 from ...utils.testing import requires_data
 from ...data import EventList, EventListDataset, EventListDatasetChecker
@@ -30,6 +31,22 @@ def test_EventList():
     assert_allclose(event_list.observation_dead_time_fraction, 0.03576320037245795)
 
 
+@requires_data('gammapy-extra')
+def test_EventList_region():
+    from gammapy.region import SkyCircleRegion, SkyRegionList
+
+    filename = gammapy_extra.filename('test_datasets/unbundled/hess/run_0023037_hard_eventlist.fits.gz')
+    event_list = EventList.read(filename, hdu='EVENTS')
+
+    pos = SkyCoord(81, 21, unit='deg', frame='icrs')
+    radius = Angle(1, 'deg')
+    circ = SkyCircleRegion(pos=pos, radius=radius)
+    region = SkyRegionList([circ])
+    filtered_list = event_list.select_circular_region(region)
+
+    assert_allclose(filtered_list[4]['RA'],81,rtol=1)
+    assert_allclose(filtered_list[2]['DEC'],21,rtol=1)
+
 @requires_data('gammapy-extra')
 def test_EventListDataset():
     filename = gammapy_extra.filename('test_datasets/unbundled/hess/run_0023037_hard_eventlist.fits.gz')

gammapy/image/__init__.py

@@ -5,3 +5,4 @@
 from .profile import *
 from .plotting import *
 from .catalog import *
+from .mask import *