gammapy.data.observations.Observations.__getitem__() - Code Metrics - Inspection of "Introducing the Observations class" - gammapy/gammapy - Measure and Improve Code Quality continuously with Scrutinizer

Passed

Pull Request — master (#1888)

unknown

created 2018-10-26 07:58 UTC

Observations.getitem() A

↳ Parent: gammapy.data.observations

Complexity

Conditions

Size

Total Lines	2
Code Lines	2

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
cc	1
eloc	2
nop	2
dl	0
loc	2
rs	10
c	0
b	0
f	0

# Licensed under a 3-clause BSD style license - see LICENSE.rst
from __future__ import absolute_import, division, print_function, unicode_literals
import logging
import numpy as np
from collections import OrderedDict
from astropy.coordinates import SkyCoord
from astropy.units import Quantity
from astropy.utils import lazyproperty
from astropy.time import Time
from .event_list import EventListChecker
from ..utils.testing import Checker
from ..utils.fits import earth_location_from_dict
from ..utils.table import table_row_to_dict
from ..utils.time import time_ref_from_dict

__all__ = ["ObservationCTA", "DataStoreObservation", "Observations"]

log = logging.getLogger(__name__)


class ObservationCTA(object):
    """Container class for an CTA observation

    Parameters follow loosely the "Required columns" here:
    http://gamma-astro-data-formats.readthedocs.io/en/latest/data_storage/obs_index/index.html#required-columns

    Parameters
    ----------
    obs_id : int
        Observation ID
    gti : `~gammapy.data.GTI`
        Good Time Intervals
    events : `~gammapy.data.EventList`
        Event list
    aeff : `~gammapy.irf.EffectiveAreaTable2D`
        Effective area
    edisp : `~gammapy.irf.EnergyDispersion2D`
        Energy dispersion
    psf : `~gammapy.irf.PSF3D` or `~gammapy.irf.EnergyDependentMultiGaussPSF` or `~gammapy.irf.PSFKing`
        Tabled Point Spread Function
    bkg : `~gammapy.irf.Background2D` or `~gammapy.irf.Background3D`
        Background rate
    pointing_radec : `~astropy.coordinates.SkyCoord`
        Pointing RA / DEC sky coordinates
    observation_live_time_duration : `~astropy.units.Quantity`
        Live-time duration in seconds

        The dead-time-corrected observation time.

        Computed as ``t_live = t_observation * (1 - f_dead)``
        where ``f_dead`` is the dead-time fraction.
    observation_dead_time_fraction : float
        Dead-time fraction

        Defined as dead-time over observation time.

        Dead-time is defined as the time during the observation
        where the detector did not record events:
        https://en.wikipedia.org/wiki/Dead_time
        https://adsabs.harvard.edu/abs/2004APh....22..285F

        The dead-time fraction is used in the live-time computation,
        which in turn is used in the exposure and flux computation.
    meta : `~collections.OrderedDict`
        Dictionary to store metadata

    Examples
    --------
    A minimal working example of how to create an observation from CTA's 1DC is given in
    examples/example_observation_cta.py

    """

    def __init__(
        self,
        obs_id=None,
        gti=None,
        events=None,
        aeff=None,
        edisp=None,
        psf=None,
        bkg=None,
        pointing_radec=None,
        observation_live_time_duration=None,
        observation_dead_time_fraction=None,
        meta=None,
    ):
        self.obs_id = obs_id
        self.gti = gti
        self.events = events
        self.aeff = aeff
        self.edisp = edisp
        self.psf = psf
        self.bkg = bkg
        self.pointing_radec = pointing_radec
        self.observation_live_time_duration = observation_live_time_duration
        self.observation_dead_time_fraction = observation_dead_time_fraction
        self.meta = meta or OrderedDict()

    def __str__(self):
        ss = "Info for OBS_ID = {}\n".format(self.obs_id)

        ss += "- Pointing pos: RA {:.2f} / Dec {:.2f}\n".format(
            self.pointing_radec.ra if self.pointing_radec else "None",
            self.pointing_radec.dec if self.pointing_radec else "None",
        )

        tstart = np.atleast_1d(self.gti.time_start.fits)[0] if self.gti else "None"
        ss += "- Start time: {}\n".format(tstart)
        ss += "- Observation duration: {}\n".format(
            self.gti.time_sum if self.gti else "None"
        )
        ss += "- Dead-time fraction: {:5.3f} %\n".format(
            100 * self.observation_dead_time_fraction
        )

        return ss


class DataStoreObservation(object):
    """IACT data store observation.

    Parameters
    ----------
    obs_id : int
        Observation ID
    data_store : `~gammapy.data.DataStore`
        Data store
    """

    def __init__(self, obs_id, data_store):
        # Assert that `obs_id` is available
        if obs_id not in data_store.obs_table["OBS_ID"]:
            raise ValueError("OBS_ID = {} not in obs index table.".format(obs_id))
        if obs_id not in data_store.hdu_table["OBS_ID"]:
            raise ValueError("OBS_ID = {} not in HDU index table.".format(obs_id))

        self.obs_id = obs_id
        self.data_store = data_store

    def __str__(self):
        ss = "Info for OBS_ID = {}\n".format(self.obs_id)
        ss += "- Start time: {:.2f}\n".format(self.tstart.mjd)
        ss += "- Pointing pos: RA {:.2f} / Dec {:.2f}\n".format(
            self.pointing_radec.ra, self.pointing_radec.dec
        )
        ss += "- Observation duration: {}\n".format(self.observation_time_duration)
        ss += "- Dead-time fraction: {:5.3f} %\n".format(
            100 * self.observation_dead_time_fraction
        )

        # TODO: Which target was observed?
        # TODO: print info about available HDUs for this observation ...
        return ss

    def location(self, hdu_type=None, hdu_class=None):
        """HDU location object.

        Parameters
        ----------
        hdu_type : str
            HDU type (see `~gammapy.data.HDUIndexTable.VALID_HDU_TYPE`)
        hdu_class : str
            HDU class (see `~gammapy.data.HDUIndexTable.VALID_HDU_CLASS`)

        Returns
        -------
        location : `~gammapy.data.HDULocation`
            HDU location
        """
        return self.data_store.hdu_table.hdu_location(
            obs_id=self.obs_id, hdu_type=hdu_type, hdu_class=hdu_class
        )

    def load(self, hdu_type=None, hdu_class=None):
        """Load data file as appropriate object.

        Parameters
        ----------
        hdu_type : str
            HDU type (see `~gammapy.data.HDUIndexTable.VALID_HDU_TYPE`)
        hdu_class : str
            HDU class (see `~gammapy.data.HDUIndexTable.VALID_HDU_CLASS`)

        Returns
        -------
        object : object
            Object depends on type, e.g. for `events` it's a `~gammapy.data.EventList`.
        """
        location = self.location(hdu_type=hdu_type, hdu_class=hdu_class)
        return location.load()

    @property
    def events(self):
        """Load `gammapy.data.EventList` object."""
        return self.load(hdu_type="events")

    @property
    def gti(self):
        """Load `gammapy.data.GTI` object."""
        return self.load(hdu_type="gti")

    @property
    def aeff(self):
        """Load effective area object."""
        return self.load(hdu_type="aeff")

    @property
    def edisp(self):
        """Load energy dispersion object."""
        return self.load(hdu_type="edisp")

    @property
    def psf(self):
        """Load point spread function object."""
        return self.load(hdu_type="psf")

    @property
    def bkg(self):
        """Load background object."""
        return self.load(hdu_type="bkg")

    @lazyproperty
    def obs_info(self):
        """Observation information (`~collections.OrderedDict`)."""
        row = self.data_store.obs_table.select_obs_id(obs_id=self.obs_id)[0]
        return table_row_to_dict(row)

    @lazyproperty
    def tstart(self):
        """Observation start time (`~astropy.time.Time`)."""
        met_ref = time_ref_from_dict(self.data_store.obs_table.meta)
        met = Quantity(self.obs_info["TSTART"].astype("float64"), "second")
        time = met_ref + met
        return time

    @lazyproperty
    def tstop(self):
        """Observation stop time (`~astropy.time.Time`)."""
        met_ref = time_ref_from_dict(self.data_store.obs_table.meta)
        met = Quantity(self.obs_info["TSTOP"].astype("float64"), "second")
        time = met_ref + met
        return time

    @lazyproperty
    def observation_time_duration(self):
        """Observation time duration in seconds (`~astropy.units.Quantity`).

        The wall time, including dead-time.
        """
        return Quantity(self.obs_info["ONTIME"], "second")

    @lazyproperty
    def observation_live_time_duration(self):
        """Live-time duration in seconds (`~astropy.units.Quantity`).

        The dead-time-corrected observation time.

        Computed as ``t_live = t_observation * (1 - f_dead)``
        where ``f_dead`` is the dead-time fraction.
        """
        return Quantity(self.obs_info["LIVETIME"], "second")

    @lazyproperty
    def observation_dead_time_fraction(self):
        """Dead-time fraction (float).

        Defined as dead-time over observation time.

        Dead-time is defined as the time during the observation
        where the detector didn't record events:
        https://en.wikipedia.org/wiki/Dead_time
        https://adsabs.harvard.edu/abs/2004APh....22..285F

        The dead-time fraction is used in the live-time computation,
        which in turn is used in the exposure and flux computation.
        """
        return 1 - self.obs_info["DEADC"]

    @lazyproperty
    def pointing_radec(self):
        """Pointing RA / DEC sky coordinates (`~astropy.coordinates.SkyCoord`)."""
        lon, lat = self.obs_info["RA_PNT"], self.obs_info["DEC_PNT"]
        return SkyCoord(lon, lat, unit="deg", frame="icrs")

    @lazyproperty
    def pointing_altaz(self):
        """Pointing ALT / AZ sky coordinates (`~astropy.coordinates.SkyCoord`)."""
        alt, az = self.obs_info["ALT_PNT"], self.obs_info["AZ_PNT"]
        return SkyCoord(az, alt, unit="deg", frame="altaz")

    @lazyproperty
    def pointing_zen(self):
        """Pointing zenith angle sky (`~astropy.units.Quantity`)."""
        return Quantity(self.obs_info["ZEN_PNT"], unit="deg")

    @lazyproperty
    def target_radec(self):
        """Target RA / DEC sky coordinates (`~astropy.coordinates.SkyCoord`)."""
        lon, lat = self.obs_info["RA_OBJ"], self.obs_info["DEC_OBJ"]
        return SkyCoord(lon, lat, unit="deg", frame="icrs")

    @lazyproperty
    def observatory_earth_location(self):
        """Observatory location (`~astropy.coordinates.EarthLocation`)."""
        return earth_location_from_dict(self.obs_info)

    @lazyproperty
    def muoneff(self):
        """Observation muon efficiency."""
        return self.obs_info["MUONEFF"]

    def peek(self):
        """Quick-look plots in a few panels."""
        raise NotImplementedError

    def to_observation_cta(self):
        """Convert to `~gammapy.data.ObservationCTA`.

        This loads all observation-related info from disk
        and stores it in the in-memory ``ObservationCTA``.

        Returns
        -------
        obs : `~gammapy.data.ObservationCTA`
            Observation
        """
        # maps the ObservationCTA class attributes to the DataStoreObservation properties
        props = {
            "obs_id": "obs_id",
            "gti": "gti",
            "events": "events",
            "aeff": "aeff",
            "edisp": "edisp",
            "psf": "psf",
            "bkg": "bkg",
            "pointing_radec": "pointing_radec",
            "observation_live_time_duration": "observation_live_time_duration",
            "observation_dead_time_fraction": "observation_dead_time_fraction",
        }

        for obs_cta_kwarg, ds_obs_prop in props.items():
            try:
                props[obs_cta_kwarg] = getattr(self, ds_obs_prop)
            except Exception as exception:
                log.warning(exception)
                props[obs_cta_kwarg] = None

        return ObservationCTA(**props)

    def check(self, checks="all"):
        """Run checks.

        This is a generator that yields a list of dicts.
        """
        checker = ObservationChecker(self)
        return checker.run(checks=checks)


class Observations(object):
    """Container class that holds a list of observations.

    Parameters:
    ----------
    obs_list : list
        A list of `~gammapy.data.DataStoreObservation`
    """

    def __init__(self, obs_list=None):
        self._obs_list = obs_list or []

    def __getitem__(self, key):
        return self._obs_list[key]

    def __len__(self):
        return len(self._obs_list)

    def __str__(self):
        s = self.__class__.__name__ + "\n"
        s += "Number of observations: {}\n".format(len(self))
        for obs in self:
            s += str(obs)
        return s


class ObservationChecker(Checker):
    """Check an observation.

    Checks data format and a bit about the content.
    """

    CHECKS = {
        "events": "check_events",
        "gti": "check_gti",
        "aeff": "check_aeff",
        "edisp": "check_edisp",
        "psf": "check_psf",
    }

    def __init__(self, obs):
        self.obs = obs

    def _record(self, level="info", msg=None):
        return {"level": level, "obs_id": self.obs.obs_id, "msg": msg}

    def check_events(self):
        yield self._record(level="debug", msg="Starting events check")

        try:
            events = self.obs.load("events")
        except Exception:
            yield self._record(level="warning", msg="Loading events failed")
            return

        for record in EventListChecker(events).run():
            yield record

    # TODO: split this out into a GTIChecker
    def check_gti(self):
        yield self._record(level="debug", msg="Starting gti check")

        try:
            gti = self.obs.load("gti")
        except Exception:
            yield self._record(level="warning", msg="Loading GTI failed")
            return

        if len(gti.table) == 0:
            yield self._record(level="error", msg="GTI table has zero rows")

        columns_required = ["START", "STOP"]
        for name in columns_required:
            if name not in gti.table.colnames:
                yield self._record(
                    level="error", msg="Missing table column: {!r}".format(name)
                )

        # TODO: Check that header keywords agree with table entries
        # TSTART, TSTOP, MJDREFI, MJDREFF

        # Check that START and STOP times are consecutive
        # times = np.ravel(self.table['START'], self.table['STOP'])
        # # TODO: not sure this is correct ... add test with a multi-gti table from Fermi.
        # if not np.all(np.diff(times) >= 0):
        #     yield 'GTIs are not consecutive or sorted.'

    # TODO: add reference times for all instruments and check for this
    # Use TELESCOP header key to check which instrument it is.
    def _check_times(self):
        """Check if various times are consistent.

        The headers and tables of the FITS EVENTS and GTI extension
        contain various observation and event time information.
        """
        # http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_Data/Time_in_ScienceTools.html
        # https://hess-confluence.desy.de/confluence/display/HESS/HESS+FITS+data+-+References+and+checks#HESSFITSdata-Referencesandchecks-Time
        telescope_met_refs = OrderedDict(
            FERMI=Time("2001-01-01T00:00:00"), HESS=Time("2001-01-01T00:00:00")
        )

        meta = self.dset.event_list.table.meta
        telescope = meta["TELESCOP"]

        if telescope in telescope_met_refs.keys():
            dt = self.time_ref - telescope_met_refs[telescope]
            if dt > self.accuracy["time"]:
                yield self._record(
                    level="error", msg="Reference time incorrect for telescope"
                )

    def check_aeff(self):
        yield self._record(level="debug", msg="Starting aeff check")

        try:
            aeff = self.obs.load("aeff")
        except Exception:
            yield self._record(level="warning", msg="Loading aeff failed")
            return

        # Check that thresholds are meaningful for aeff
        if (
            "LO_THRES" in aeff.meta
            and "HI_THRES" in aeff.meta
            and aeff.meta["LO_THRES"] >= aeff.meta["HI_THRES"]
        ):
            yield self._record(
                level="error", msg="LO_THRES >= HI_THRES in effective area meta data"
            )

        # Check that data isn't all null
        if np.max(aeff.data.data) <= 0:
            yield self._record(
                level="error", msg="maximum entry of effective area is <= 0"
            )

    def check_edisp(self):
        yield self._record(level="debug", msg="Starting edisp check")

        try:
            edisp = self.obs.load("edisp")
        except Exception:
            yield self._record(level="warning", msg="Loading edisp failed")
            return

        # Check that data isn't all null
        if np.max(edisp.data.data) <= 0:
            yield self._record(level="error", msg="maximum entry of edisp is <= 0")

    def check_psf(self):
        yield self._record(level="debug", msg="Starting psf check")

        try:
            psf = self.obs.load("psf")
        except Exception:
            yield self._record(level="warning", msg="Loading psf failed")
            return

        # TODO: implement some basic check
        # The following doesn't work, because EnergyDependentMultiGaussPSF
        # has no attribute `data`
        # Check that data isn't all null
        # if np.max(psf.data.data) <= 0:
        #     yield self._record(
        #         level="error", msg="maximum entry of psf is <= 0"
        #     )


1			# Licensed under a 3-clause BSD style license - see LICENSE.rst
2			from __future__ import absolute_import, division, print_function, unicode_literals
3			import logging
4			import numpy as np
5			from collections import OrderedDict
6			from astropy.coordinates import SkyCoord
7			from astropy.units import Quantity
8			from astropy.utils import lazyproperty
9			from astropy.time import Time
10			from .event_list import EventListChecker
11			from ..utils.testing import Checker
12			from ..utils.fits import earth_location_from_dict
13			from ..utils.table import table_row_to_dict
14			from ..utils.time import time_ref_from_dict
15
16			__all__ = ["ObservationCTA", "DataStoreObservation", "Observations"]
17
18			log = logging.getLogger(__name__)
19
20
21			class ObservationCTA(object):
22			"""Container class for an CTA observation
23
24			Parameters follow loosely the "Required columns" here:
25			http://gamma-astro-data-formats.readthedocs.io/en/latest/data_storage/obs_index/index.html#required-columns
26
27			Parameters
28			----------
29			obs_id : int
30			Observation ID
31			gti : `~gammapy.data.GTI`
32			Good Time Intervals
33			events : `~gammapy.data.EventList`
34			Event list
35			aeff : `~gammapy.irf.EffectiveAreaTable2D`
36			Effective area
37			edisp : `~gammapy.irf.EnergyDispersion2D`
38			Energy dispersion
39			psf : `~gammapy.irf.PSF3D` or `~gammapy.irf.EnergyDependentMultiGaussPSF` or `~gammapy.irf.PSFKing`
40			Tabled Point Spread Function
41			bkg : `~gammapy.irf.Background2D` or `~gammapy.irf.Background3D`
42			Background rate
43			pointing_radec : `~astropy.coordinates.SkyCoord`
44			Pointing RA / DEC sky coordinates
45			observation_live_time_duration : `~astropy.units.Quantity`
46			Live-time duration in seconds
47
48			The dead-time-corrected observation time.
49
50			Computed as ``t_live = t_observation * (1 - f_dead)``
51			where ``f_dead`` is the dead-time fraction.
52			observation_dead_time_fraction : float
53			Dead-time fraction
54
55			Defined as dead-time over observation time.
56
57			Dead-time is defined as the time during the observation
58			where the detector did not record events:
59			https://en.wikipedia.org/wiki/Dead_time
60			https://adsabs.harvard.edu/abs/2004APh....22..285F
61
62			The dead-time fraction is used in the live-time computation,
63			which in turn is used in the exposure and flux computation.
64			meta : `~collections.OrderedDict`
65			Dictionary to store metadata
66
67			Examples
68			--------
69			A minimal working example of how to create an observation from CTA's 1DC is given in
70			examples/example_observation_cta.py
71
72			"""
73
74			def __init__(
75			self,
76			obs_id=None,
77			gti=None,
78			events=None,
79			aeff=None,
80			edisp=None,
81			psf=None,
82			bkg=None,
83			pointing_radec=None,
84			observation_live_time_duration=None,
85			observation_dead_time_fraction=None,
86			meta=None,
87			):
88			self.obs_id = obs_id
89			self.gti = gti
90			self.events = events
91			self.aeff = aeff
92			self.edisp = edisp
93			self.psf = psf
94			self.bkg = bkg
95			self.pointing_radec = pointing_radec
96			self.observation_live_time_duration = observation_live_time_duration
97			self.observation_dead_time_fraction = observation_dead_time_fraction
98			self.meta = meta or OrderedDict()
99
100			def __str__(self):
101			ss = "Info for OBS_ID = {}\n".format(self.obs_id)
102
103			ss += "- Pointing pos: RA {:.2f} / Dec {:.2f}\n".format(
104			self.pointing_radec.ra if self.pointing_radec else "None",
105			self.pointing_radec.dec if self.pointing_radec else "None",
106			)
107
108			tstart = np.atleast_1d(self.gti.time_start.fits)[0] if self.gti else "None"
109			ss += "- Start time: {}\n".format(tstart)
110			ss += "- Observation duration: {}\n".format(
111			self.gti.time_sum if self.gti else "None"
112			)
113			ss += "- Dead-time fraction: {:5.3f} %\n".format(
114			100 * self.observation_dead_time_fraction
115			)
116
117			return ss
118
119
120			class DataStoreObservation(object):
121			"""IACT data store observation.
122
123			Parameters
124			----------
125			obs_id : int
126			Observation ID
127			data_store : `~gammapy.data.DataStore`
128			Data store
129			"""
130
131			def __init__(self, obs_id, data_store):
132			# Assert that `obs_id` is available
133			if obs_id not in data_store.obs_table["OBS_ID"]:
134			raise ValueError("OBS_ID = {} not in obs index table.".format(obs_id))
135			if obs_id not in data_store.hdu_table["OBS_ID"]:
136			raise ValueError("OBS_ID = {} not in HDU index table.".format(obs_id))
137
138			self.obs_id = obs_id
139			self.data_store = data_store
140
141			def __str__(self):
142			ss = "Info for OBS_ID = {}\n".format(self.obs_id)
143			ss += "- Start time: {:.2f}\n".format(self.tstart.mjd)
144			ss += "- Pointing pos: RA {:.2f} / Dec {:.2f}\n".format(
145			self.pointing_radec.ra, self.pointing_radec.dec
146			)
147			ss += "- Observation duration: {}\n".format(self.observation_time_duration)
148			ss += "- Dead-time fraction: {:5.3f} %\n".format(
149			100 * self.observation_dead_time_fraction
150			)
151
152			# TODO: Which target was observed?
153			# TODO: print info about available HDUs for this observation ...
154			return ss
155
156			def location(self, hdu_type=None, hdu_class=None):
157			"""HDU location object.
158
159			Parameters
160			----------
161			hdu_type : str
162			HDU type (see `~gammapy.data.HDUIndexTable.VALID_HDU_TYPE`)
163			hdu_class : str
164			HDU class (see `~gammapy.data.HDUIndexTable.VALID_HDU_CLASS`)
165
166			Returns
167			-------
168			location : `~gammapy.data.HDULocation`
169			HDU location
170			"""
171			return self.data_store.hdu_table.hdu_location(
172			obs_id=self.obs_id, hdu_type=hdu_type, hdu_class=hdu_class
173			)
174
175			def load(self, hdu_type=None, hdu_class=None):
176			"""Load data file as appropriate object.
177
178			Parameters
179			----------
180			hdu_type : str
181			HDU type (see `~gammapy.data.HDUIndexTable.VALID_HDU_TYPE`)
182			hdu_class : str
183			HDU class (see `~gammapy.data.HDUIndexTable.VALID_HDU_CLASS`)
184
185			Returns
186			-------
187			object : object
188			Object depends on type, e.g. for `events` it's a `~gammapy.data.EventList`.
189			"""
190			location = self.location(hdu_type=hdu_type, hdu_class=hdu_class)
191			return location.load()
192
193			@property
194			def events(self):
195			"""Load `gammapy.data.EventList` object."""
196			return self.load(hdu_type="events")
197
198			@property
199			def gti(self):
200			"""Load `gammapy.data.GTI` object."""
201			return self.load(hdu_type="gti")
202
203			@property
204			def aeff(self):
205			"""Load effective area object."""
206			return self.load(hdu_type="aeff")
207
208			@property
209			def edisp(self):
210			"""Load energy dispersion object."""
211			return self.load(hdu_type="edisp")
212
213			@property
214			def psf(self):
215			"""Load point spread function object."""
216			return self.load(hdu_type="psf")
217
218			@property
219			def bkg(self):
220			"""Load background object."""
221			return self.load(hdu_type="bkg")
222
223			@lazyproperty
224			def obs_info(self):
225			"""Observation information (`~collections.OrderedDict`)."""
226			row = self.data_store.obs_table.select_obs_id(obs_id=self.obs_id)[0]
227			return table_row_to_dict(row)
228
229			@lazyproperty
230			def tstart(self):
231			"""Observation start time (`~astropy.time.Time`)."""
232			met_ref = time_ref_from_dict(self.data_store.obs_table.meta)
233			met = Quantity(self.obs_info["TSTART"].astype("float64"), "second")
234			time = met_ref + met
235			return time
236
237			@lazyproperty
238			def tstop(self):
239			"""Observation stop time (`~astropy.time.Time`)."""
240			met_ref = time_ref_from_dict(self.data_store.obs_table.meta)
241			met = Quantity(self.obs_info["TSTOP"].astype("float64"), "second")
242			time = met_ref + met
243			return time
244
245			@lazyproperty
246			def observation_time_duration(self):
247			"""Observation time duration in seconds (`~astropy.units.Quantity`).
248
249			The wall time, including dead-time.
250			"""
251			return Quantity(self.obs_info["ONTIME"], "second")
252
253			@lazyproperty
254			def observation_live_time_duration(self):
255			"""Live-time duration in seconds (`~astropy.units.Quantity`).
256
257			The dead-time-corrected observation time.
258
259			Computed as ``t_live = t_observation * (1 - f_dead)``
260			where ``f_dead`` is the dead-time fraction.
261			"""
262			return Quantity(self.obs_info["LIVETIME"], "second")
263
264			@lazyproperty
265			def observation_dead_time_fraction(self):
266			"""Dead-time fraction (float).
267
268			Defined as dead-time over observation time.
269
270			Dead-time is defined as the time during the observation
271			where the detector didn't record events:
272			https://en.wikipedia.org/wiki/Dead_time
273			https://adsabs.harvard.edu/abs/2004APh....22..285F
274
275			The dead-time fraction is used in the live-time computation,
276			which in turn is used in the exposure and flux computation.
277			"""
278			return 1 - self.obs_info["DEADC"]
279
280			@lazyproperty
281			def pointing_radec(self):
282			"""Pointing RA / DEC sky coordinates (`~astropy.coordinates.SkyCoord`)."""
283			lon, lat = self.obs_info["RA_PNT"], self.obs_info["DEC_PNT"]
284			return SkyCoord(lon, lat, unit="deg", frame="icrs")
285
286			@lazyproperty
287			def pointing_altaz(self):
288			"""Pointing ALT / AZ sky coordinates (`~astropy.coordinates.SkyCoord`)."""
289			alt, az = self.obs_info["ALT_PNT"], self.obs_info["AZ_PNT"]
290			return SkyCoord(az, alt, unit="deg", frame="altaz")
291
292			@lazyproperty
293			def pointing_zen(self):
294			"""Pointing zenith angle sky (`~astropy.units.Quantity`)."""
295			return Quantity(self.obs_info["ZEN_PNT"], unit="deg")
296
297			@lazyproperty
298			def target_radec(self):
299			"""Target RA / DEC sky coordinates (`~astropy.coordinates.SkyCoord`)."""
300			lon, lat = self.obs_info["RA_OBJ"], self.obs_info["DEC_OBJ"]
301			return SkyCoord(lon, lat, unit="deg", frame="icrs")
302
303			@lazyproperty
304			def observatory_earth_location(self):
305			"""Observatory location (`~astropy.coordinates.EarthLocation`)."""
306			return earth_location_from_dict(self.obs_info)
307
308			@lazyproperty
309			def muoneff(self):
310			"""Observation muon efficiency."""
311			return self.obs_info["MUONEFF"]
312
313			def peek(self):
314			"""Quick-look plots in a few panels."""
315			raise NotImplementedError
316
317			def to_observation_cta(self):
318			"""Convert to `~gammapy.data.ObservationCTA`.
319
320			This loads all observation-related info from disk
321			and stores it in the in-memory ``ObservationCTA``.
322
323			Returns
324			-------
325			obs : `~gammapy.data.ObservationCTA`
326			Observation
327			"""
328			# maps the ObservationCTA class attributes to the DataStoreObservation properties
329			props = {
330			"obs_id": "obs_id",
331			"gti": "gti",
332			"events": "events",
333			"aeff": "aeff",
334			"edisp": "edisp",
335			"psf": "psf",
336			"bkg": "bkg",
337			"pointing_radec": "pointing_radec",
338			"observation_live_time_duration": "observation_live_time_duration",
339			"observation_dead_time_fraction": "observation_dead_time_fraction",
340			}
341
342			for obs_cta_kwarg, ds_obs_prop in props.items():
343			try:
344			props[obs_cta_kwarg] = getattr(self, ds_obs_prop)
345			except Exception as exception:
346			log.warning(exception)
347			props[obs_cta_kwarg] = None
348
349			return ObservationCTA(**props)
350
351			def check(self, checks="all"):
352			"""Run checks.
353
354			This is a generator that yields a list of dicts.
355			"""
356			checker = ObservationChecker(self)
357			return checker.run(checks=checks)
358
359
360			class Observations(object):
361			"""Container class that holds a list of observations.
362
363			Parameters:
364			----------
365			obs_list : list
366			A list of `~gammapy.data.DataStoreObservation`
367			"""
368
369			def __init__(self, obs_list=None):
370			self._obs_list = obs_list or []
371
372			def __getitem__(self, key):
373			return self._obs_list[key]
374
375			def __len__(self):
376			return len(self._obs_list)
377
378			def __str__(self):
379			s = self.__class__.__name__ + "\n"
380			s += "Number of observations: {}\n".format(len(self))
381			for obs in self:
382			s += str(obs)
383			return s
384
385
386			class ObservationChecker(Checker):
387			"""Check an observation.
388
389			Checks data format and a bit about the content.
390			"""
391
392			CHECKS = {
393			"events": "check_events",
394			"gti": "check_gti",
395			"aeff": "check_aeff",
396			"edisp": "check_edisp",
397			"psf": "check_psf",
398			}
399
400			def __init__(self, obs):
401			self.obs = obs
402
403			def _record(self, level="info", msg=None):
404			return {"level": level, "obs_id": self.obs.obs_id, "msg": msg}
405
406			def check_events(self):
407			yield self._record(level="debug", msg="Starting events check")
408
409			try:
410			events = self.obs.load("events")
411			except Exception:
412			yield self._record(level="warning", msg="Loading events failed")
413			return
414
415			for record in EventListChecker(events).run():
416			yield record
417
418			# TODO: split this out into a GTIChecker
419			def check_gti(self):
420			yield self._record(level="debug", msg="Starting gti check")
421
422			try:
423			gti = self.obs.load("gti")
424			except Exception:
425			yield self._record(level="warning", msg="Loading GTI failed")
426			return
427
428			if len(gti.table) == 0:
429			yield self._record(level="error", msg="GTI table has zero rows")
430
431			columns_required = ["START", "STOP"]
432			for name in columns_required:
433			if name not in gti.table.colnames:
434			yield self._record(
435			level="error", msg="Missing table column: {!r}".format(name)
436			)
437
438			# TODO: Check that header keywords agree with table entries
439			# TSTART, TSTOP, MJDREFI, MJDREFF
440
441			# Check that START and STOP times are consecutive
442			# times = np.ravel(self.table['START'], self.table['STOP'])
443			# # TODO: not sure this is correct ... add test with a multi-gti table from Fermi.
444			# if not np.all(np.diff(times) >= 0):
445			# yield 'GTIs are not consecutive or sorted.'
446
447			# TODO: add reference times for all instruments and check for this
448			# Use TELESCOP header key to check which instrument it is.
449			def _check_times(self):
450			"""Check if various times are consistent.
451
452			The headers and tables of the FITS EVENTS and GTI extension
453			contain various observation and event time information.
454			"""
455			# http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_Data/Time_in_ScienceTools.html
456			# https://hess-confluence.desy.de/confluence/display/HESS/HESS+FITS+data+-+References+and+checks#HESSFITSdata-Referencesandchecks-Time
457			telescope_met_refs = OrderedDict(
458			FERMI=Time("2001-01-01T00:00:00"), HESS=Time("2001-01-01T00:00:00")
459			)
460
461			meta = self.dset.event_list.table.meta
462			telescope = meta["TELESCOP"]
463
464			if telescope in telescope_met_refs.keys():
465			dt = self.time_ref - telescope_met_refs[telescope]
466			if dt > self.accuracy["time"]:
467			yield self._record(
468			level="error", msg="Reference time incorrect for telescope"
469			)
470
471			def check_aeff(self):
472			yield self._record(level="debug", msg="Starting aeff check")
473
474			try:
475			aeff = self.obs.load("aeff")
476			except Exception:
477			yield self._record(level="warning", msg="Loading aeff failed")
478			return
479
480			# Check that thresholds are meaningful for aeff
481			if (
482			"LO_THRES" in aeff.meta
483			and "HI_THRES" in aeff.meta
484			and aeff.meta["LO_THRES"] >= aeff.meta["HI_THRES"]
485			):
486			yield self._record(
487			level="error", msg="LO_THRES >= HI_THRES in effective area meta data"
488			)
489
490			# Check that data isn't all null
491			if np.max(aeff.data.data) <= 0:
492			yield self._record(
493			level="error", msg="maximum entry of effective area is <= 0"
494			)
495
496			def check_edisp(self):
497			yield self._record(level="debug", msg="Starting edisp check")
498
499			try:
500			edisp = self.obs.load("edisp")
501			except Exception:
502			yield self._record(level="warning", msg="Loading edisp failed")
503			return
504
505			# Check that data isn't all null
506			if np.max(edisp.data.data) <= 0:
507			yield self._record(level="error", msg="maximum entry of edisp is <= 0")
508
509			def check_psf(self):
510			yield self._record(level="debug", msg="Starting psf check")
511
512			try:
513			psf = self.obs.load("psf")
514			except Exception:
515			yield self._record(level="warning", msg="Loading psf failed")
516			return
517
518			# TODO: implement some basic check
519			# The following doesn't work, because EnergyDependentMultiGaussPSF
520			# has no attribute `data`
521			# Check that data isn't all null
522			# if np.max(psf.data.data) <= 0:
523			# yield self._record(
524			# level="error", msg="maximum entry of psf is <= 0"
525			# )
526

gammapy / gammapy

Pull Request — master (#1888)

Observations.__getitem__() A

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like

Observations.getitem() A