Add Target class

docs/region/make_reflected_regions.py

@@ -1,7 +1,8 @@
 from astropy.coordinates import SkyCoord, Angle
 from astropy.wcs import WCS
 
-from gammapy.region import find_reflected_regions, SkyCircleRegion
+from gammapy.region import find_reflected_regions
+from gammapy.extern.regions.shapes import CircleSkyRegion
 from gammapy.image import ExclusionMask
 from gammapy.utils.testing import requires_data
 
@@ -11,7 +12,7 @@
 
 pos = SkyCoord(80.2, 23.5, unit='deg', frame='icrs')
 radius = Angle(0.4, 'deg')
-test_region = SkyCircleRegion(pos=pos, radius=radius)
+test_region = CircleSkyRegion(pos, radius)
 center = SkyCoord(82.8, 22.5, unit='deg', frame='icrs')
 regions = find_reflected_regions(test_region, center, mask)
 

gammapy/background/__init__.py

@@ -3,6 +3,7 @@
 Background estimation and modeling methods.
 """
 from .energy_offset_array import *
+from .background_estimate import *
 from .fov import *
 from .cube import *
 from .on_off import *

gammapy/background/background_estimate.py

@@ -0,0 +1,85 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+from __future__ import absolute_import, division, print_function, unicode_literals
+
+from .ring import ring_area_factor
+from ..region import find_reflected_regions
+
+__all__ = [
+    'BackgroundEstimate',
+    'ring_background_estimate',
+    'reflected_regions_background_estimate',
+]
+
+
+class BackgroundEstimate(object):
+    """Container class for background estimate
+
+    This container holds the result from a region based background estimation
+    for one observation. Currently, it is filled by the functions
+    :func:`~gammapy.background.ring_background_estimate` and
+    :func:`~gammapy.background.reflected_regions_background_estimate`. At some
+    points this should be replaced by a more elaborate analysis class.
+
+    Parameters
+    ----------
+    off_events : `~gammapy.data.EventList`
+        Background events
+    off_region : `~gammapy.extern.regions.SkyRegion`
+        Background extraction region
+    alpha : float
+        Background scaling factor
+    tag : str
+        Background estimation method
+    """
+    def __init__(self, off_region, off_events, alpha, tag='default'):
+        self.off_region = off_region
+        self.off_events = off_events
+        self.alpha = alpha
+        self.tag = tag
+
+
+def ring_background_estimate(pos, on_radius, inner_radius, outer_radius, events):
+    """Simple ring background estimate
+
+    No acceptance correction is applied
+
+    TODO : Replace with AnnulusSkyRegion
+
+    Parameters
+    ----------
+    pos : `~astropy.coordinates.SkyCoord`
+        Target position
+    on_radius : `~astropy.coordinates.Angle`
+        On region radius
+    inner_radius : `~astropy.coordinates.Angle`
+        Inner ring radius
+    outer_radius : `~astropy.coordinates.Angle`
+        Outer ring radius
+    events : `~gammapy.data.EventList`
+        Events
+    """
+    off_events = events.select_sky_ring(pos, inner_radius, outer_radius)
+    off_region = dict(inner=inner_radius, outer=outer_radius)
+    alpha = ring_area_factor(on_radius, inner_radius, outer_radius)
+
+    return BackgroundEstimate(off_region, off_events, alpha, tag='ring')
+
+def reflected_regions_background_estimate(on_region, pointing, exclusion, events):
+    """Reflected regions background estimate
+
+    Parameters
+    ----------
+    on_region : `~gammapy.extern.regions.CircleSkyRegion`
+        On region
+    pointing : `~astropy.coordinates.SkyCoord`
+        Pointing position
+    exclusion : `~gammapy.image.SkyMap`
+        Exclusion mask
+    events : `gammapy.data.EventList`
+        Events
+    """
+    off_region = find_reflected_regions(on_region, pointing, exclusion)
+    off_events = events.select_circular_region(off_region)
+    alpha = len(off_region)
+
+    return BackgroundEstimate(off_region, off_events, alpha, tag='reflected')

gammapy/data/__init__.py

@@ -6,6 +6,7 @@
 class InvalidDataError(Exception):
     """Invalid data found."""
 
+from .target import *
 from .data_manager import *
 from .data_store import *
 from .event_list import *

gammapy/data/target.py

@@ -0,0 +1,217 @@
+# 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 ..stats import Stats
+from astropy.table import vstack as table_vstack
+from astropy.coordinates import SkyCoord
+from ..extern.regions.shapes import CircleSkyRegion
+import astropy.units as u
+
+__all__ = [
+    'Target',
+    'TargetSummary',
+]
+
+class Target(object):
+    """Observation Target.
+
+    This class represents an observation target. It serves as input for several
+    analysis classes, e.g. `~gammapy.spectrum.SpectrumExtraction` and
+    `~gammapy.data.TargetSummary`. Some analyses, e.g. background estimation,
+    are functions on the ``Target`` class, but should be refactored as separate
+    analysis classes. Each analysis class can add attributes to the ``Target``
+    instance in order to make some analysis steps, e.g. background estimation
+    reusable.
+
+    Parameters
+    ----------
+    position : `~astropy.coordinates.SkyCoord`
+        Target position
+    on_region : `~astropy.regions.SkyRegion`
+        Signal extraction region
+    obs_id : int, optional
+        Observatinos for this target
+    name : str, optional
+        Target name
+    tag : str, optional
+        Target identifier
+
+    Examples
+    --------
+    Initialize target and define observations
+
+    >>> import astropy.units as u
+    >>> from astropy.coordinates import SkyCoord
+    >>> from gammapy.extern.regions.shapes import CircleSkyRegion
+    >>> from gammapy.data import Target
+    >>> pos = SkyCoord(83.63 * u.deg, 22.01 * u.deg, frame='icrs')
+    >>> on_size = 0.3 * u.deg
+    >>> on_region = CircleSkyRegion(pos, on_size)
+    >>> target = Target(pos, on_region, name='Crab Nebula', tag='crab')
+    >>> print(target)
+    Target: Crab Nebula
+    Tag: crab
+    Position: <SkyCoord (ICRS): (ra, dec) in deg
+        (83.63, 22.01)>
+    On region: CircleSkyRegion
+        Center:<SkyCoord (ICRS): (ra, dec) in deg
+            (83.63, 22.01)>
+        Radius:0.3 deg
+
+    Add data and background estimate
+
+    >>> from gammapy.data import DataStore
+    >>> store_dir = "$GAMMAPY_EXTRA/datasets/hess-crab4-hd-hap-prod2")
+    >>> data_store = DataStore.from_dir(store_dir)
+    >>> target.obs_id = [23523, 23592]
+    >>> target.add_obs_from_store(data_store)
+    >>> print(target.obs[0])
+    Info for OBS_ID = 23523
+    - Start time: 53343.92 s
+    - Pointing pos: RA 83.63 deg / Dec 21.51 deg
+    - Observation duration: 1687.0 s
+    - Dead-time fraction: 6.240 %
+    >>> target.estimate_background(method='ring', inner_radius=inner_radius, outer_radius=outer_radius)
+    >>> print(len(target.background[0].off_events))
+    37
+
+    Create `~gammapy.data.TargetSummary` and `~gammapy.stats.data.Stats`
+
+    >>> summary = target.summary
+    >>> type(summary)
+    <class 'gammapy.data.target.TargetSummary'>
+    >>> stats = summary.stats
+    >>> type(stats)
+    <class 'gammapy.stats.data.Stats'>
+
+    """
+    def __init__(self, position, on_region, obs_id=None, name='Target', tag='target'):
+        self.position = position
+        self.on_region = on_region
+        self.name = name
+        self.tag = tag
+        self.obs_id = obs_id
+
+    def __str__(self):
+        """String representation"""
+        ss = "Target: {}\n".format(self.name)
+        ss += "Tag: {}\n".format(self.tag)
+        ss += "Position: {}\n".format(self.position)
+        ss += "On region: {}\n".format(self.on_region)
+        return ss
+
+    @property
+    def summary(self):
+        """`~gammapy.data.TargetSummary`"""
+        return TargetSummary(self)
+
+    @classmethod
+    def from_config(cls, config):
+        """Initialize target from config
+
+        The config dict is stored as attribute for later use by other analysis
+        classes.
+        """
+        obs_id = config['obs']
+        if not isinstance(obs_id, list):
+
+            from . import ObservationTable
+            obs_table = ObservationTable.read(obs_id)
+            obs_id = obs_table['OBS_ID'].data
+        # TODO : This should also accept also Galactic coordinates
+        pos = SkyCoord(config['ra'], config['dec'], unit='deg')
+        on_radius = config['on_size'] * u.deg
+        on_region = CircleSkyRegion(pos, on_radius)
+        target = cls(pos, on_region, obs_id=obs_id,
+                     name=config['name'], tag=config['tag'])
+        target.config = config
+        return target
+
+    def add_obs_from_store(self, data_store):
+        """Add a list of `~gammapy.data.DataStoreObservations`"""
+        if self.obs_id is None:
+            raise ValueError("No observations specified")
+        self.obs = [data_store.obs(_) for _ in self.obs_id]
+
+    # TODO: This should probably go into a separat BackgroundEstimator class
+    def estimate_background(self, method='ring', **kwargs):
+        """Wrapper for different background estimators"""
+
+        if method == 'ring':
+            from gammapy.background import ring_background_estimate as ring
+            pos = self.position
+            on_radius = self.on_region.radius
+            inner_radius = kwargs['inner_radius']
+            outer_radius = kwargs['outer_radius']
+            self.background = [ring(pos, on_radius, inner_radius, outer_radius,
+                                     _.events) for _ in self.obs]
+        elif method == 'reflected':
+            on_region = self.on_region
+            exclusion = kwargs['exclusion']
+            from gammapy.background import reflected_regions_background_estimate as refl
+            bkg = [refl(on_region, _.pointing_radec, exclusion,
+                        _.events) for _ in self.obs]
+            self.background = bkg
+        else:
+            raise NotImplementedError('{}'.format(method))
+
+    def run_spectral_analysis(self, outdir=None):
+        """Run spectral analysis
+
+        This runs a spectral analysis with the parameters attached as config
+        dict
+
+        Parameters
+        ----------
+        outdir : Path
+            Analysis dir
+        """
+        from . import DataStore
+        from ..image import ExclusionMask
+        from ..spectrum import SpectrumExtraction
+
+        conf = self.config
+        data_store = DataStore.from_all(conf['datastore'])
+        self.add_obs_from_store(data_store)
+        exclusion = ExclusionMask.read(conf['exclusion_mask']) or None
+        conf.update(exclusion=exclusion)
+        self.estimate_background(method=conf['background_method'], **conf)
+        # Use default energy binning
+        self.extraction = SpectrumExtraction(self)
+        self.extraction.run(outdir=outdir)
+
+
+class TargetSummary(object):
+    """Summary Info for an observation Target
+
+    For an examples see `~gammapy.data.Target`
+
+    Parameters
+    ----------
+    target : `~gammapy.data.Target`
+        Observation target
+    """
+
+    def __init__(self, target):
+        self.target = target
+
+    @property
+    def stats(self):
+        """Total `~gammapy.stats.Stats`"""
+        if self.target.background is None:
+            raise ValueError("Need Background estimate for target")
+
+        idx = self.target.on_region.contains(self.events.radec)
+        on_events = self.events[idx]
+        n_on = len(on_events)
+
+        n_off = np.sum([len(_.off_events) for _ in self.target.background])
+        # FIXME : This is only true for the ring bg
+        alpha = self.target.background[0].alpha
+        return Stats(n_on, n_off, 1, alpha)
+
+    @property
+    def events(self):
+        """All events"""
+        return table_vstack([_.events for _ in self.target.obs])

gammapy/data/tests/test_target.py

@@ -0,0 +1,28 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+from __future__ import absolute_import, division, print_function, unicode_literals
+from .. import Target, TargetSummary
+from astropy.coordinates import SkyCoord
+from astropy.tests.helper import pytest
+import astropy.units as u
+from ...extern.regions.shapes import CircleSkyRegion
+from ...utils.testing import data_manager, requires_data, requires_dependency
+
+@requires_dependency('yaml')
+@requires_data('gammapy-extra')
+def test_targetsummary(data_manager):
+    pos = SkyCoord(83.63 * u.deg, 22.01 * u.deg, frame='icrs')
+    on_size = 0.3 * u.deg
+    on_region = CircleSkyRegion(pos, on_size)
+    target = Target(pos, on_region, name='Test Target', obs_id=[23523, 23592])
+
+    data_store = data_manager['hess-crab4-hd-hap-prod2']
+    target.add_obs_from_store(data_store) 
+
+    irad = 0.5 * u.deg
+    orad = 0.7 * u.deg
+    target.estimate_background(method='ring', inner_radius=irad, outer_radius=orad)
+
+    summary = TargetSummary(target)
+
+    stats = summary.stats
+    assert stats.n_on == 432

gammapy/extern/regions/__init__.py

@@ -3,15 +3,6 @@
 """
 This is an experimental package for representing regions
 """
-
-# Affiliated packages may add whatever they like to this file, but
-# should keep this content at the top.
-# ----------------------------------------------------------------------------
-from ._astropy_init import *
-# ----------------------------------------------------------------------------
-
-# For egg_info test builds to pass, put package imports here.
-if not _ASTROPY_SETUP_:
-    from . import io
-    from . import shapes
-    from . import core
+from . import io
+from . import shapes
+from . import core

gammapy/irf/effective_area.py

@@ -107,7 +107,6 @@ def to_table(self):
 
         http://gamma-astro-data-formats.readthedocs.io/en/latest/ogip/index.html#arf-file 
         """
-
         ener_lo = self.energy.data[:-1]
         ener_hi = self.energy.data[1:]
         names = ['ENERG_LO', 'ENERG_HI', 'SPECRESP']