Take observation time from GTI table

gammapy/data/obs_table.py

@@ -1,6 +1,6 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 from __future__ import absolute_import, division, print_function, unicode_literals
-from collections import namedtuple
+from collections import namedtuple, OrderedDict
 import numpy as np
 from astropy.table import Table
 from astropy.units import Unit, Quantity
@@ -10,6 +10,7 @@
 from ..utils.scripts import make_path
 from ..utils.time import time_relative_to_ref
 from ..utils.testing import Checker
+from .gti import GTI
 
 __all__ = ["ObservationTable"]
 
@@ -351,6 +352,41 @@ def select_observations(self, selection=None):
         else:
             raise ValueError("Invalid selection type: {}".format(selection["type"]))
 
+    def create_gti(self, obs_id):
+        """Returns a GTI table containing TSTART and TSTOP from the table.
+
+        TODO: This method can be removed once GTI tables are explicitly required in Gammapy.
+
+        Parameters
+        ----------
+        obs_id : int
+            ID of the observation for which the GTI table will be created
+
+        Returns
+        -------
+        gti : `~gammapy.data.GTI`
+            GTI table containing one row (TSTART and TSTOP of the observation with `obs_id`)
+        """
+        header = OrderedDict(
+            EXTNAME="GTI",
+            HDUCLASS="GADF",
+            HDUDOC="https://github.com/open-gamma-ray-astro/gamma-astro-data-formats",
+            HDUVERS="0.2",
+            HDUCLAS1="GTI",
+            MJDREFI=self.meta["MJDREFI"],
+            MJDREFF=self.meta["MJDREFF"],
+            TIMEUNIT=self.meta.get("TIMEUNIT", "s"),
+            TIMESYS=self.meta["TIMESYS"],
+            TIMEREF=self.meta.get("TIMEREF", "LOCAL"),
+        )
+
+        idx = self.get_obs_idx(obs_id)[0]
+        start = self["TSTART"][idx].astype("float64")
+        stop = self["TSTOP"][idx].astype("float64")
+        gti_table = Table([[start], [stop]], names=("START", "STOP"), meta=header)
+
+        return GTI(gti_table)
+
 
 class ObservationTableChecker(Checker):
     """Event list checker.

gammapy/data/observations.py

@@ -198,7 +198,10 @@ def events(self):
     @property
     def gti(self):
         """Load `gammapy.data.GTI` object."""
-        return self.load(hdu_type="gti")
+        try:
+            return self.load(hdu_type="gti")
+        except IndexError:  # HDU index file does not contain the GTI table. We catch this for backward compatibility.
+            return self.data_store.obs_table.create_gti(obs_id=self.obs_id)
 
     @property
     def aeff(self):
@@ -220,37 +223,37 @@ def bkg(self):
         """Load background object."""
         return self.load(hdu_type="bkg")
 
-    @lazyproperty
+    @property
     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
+    @property
     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
+    @property
     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
+    @property
     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")
+        return self.gti.time_sum
 
-    @lazyproperty
+    @property
     def observation_live_time_duration(self):
         """Live-time duration in seconds (`~astropy.units.Quantity`).
 
@@ -259,9 +262,9 @@ def observation_live_time_duration(self):
         Computed as ``t_live = t_observation * (1 - f_dead)``
         where ``f_dead`` is the dead-time fraction.
         """
-        return Quantity(self.obs_info["LIVETIME"], "second")
+        return self.gti.time_sum * (1 - self.observation_dead_time_fraction)
 
-    @lazyproperty
+    @property
     def observation_dead_time_fraction(self):
         """Dead-time fraction (float).
 
@@ -277,35 +280,35 @@ def observation_dead_time_fraction(self):
         """
         return 1 - self.obs_info["DEADC"]
 
-    @lazyproperty
+    @property
     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
+    @property
     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
+    @property
     def pointing_zen(self):
         """Pointing zenith angle sky (`~astropy.units.Quantity`)."""
         return Quantity(self.obs_info["ZEN_PNT"], unit="deg")
 
-    @lazyproperty
+    @property
     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
+    @property
     def observatory_earth_location(self):
         """Observatory location (`~astropy.coordinates.EarthLocation`)."""
         return earth_location_from_dict(self.obs_info)
 
-    @lazyproperty
+    @property
     def muoneff(self):
         """Observation muon efficiency."""
         return self.obs_info["MUONEFF"]

gammapy/data/tests/test_obs_table.py

@@ -4,7 +4,12 @@
 from astropy.units import Quantity
 from astropy.coordinates import Angle, SkyCoord, AltAz
 from astropy.time import Time, TimeDelta
-from ...utils.testing import requires_data
+from ...data import GTI
+from ...utils.testing import (
+    requires_data,
+    assert_time_allclose,
+    assert_quantity_allclose,
+)
 from ...utils.time import time_ref_from_dict, time_relative_to_ref
 from ...utils.random import sample_sphere, get_random_state
 from ...catalog import skycoord_from_table
@@ -78,6 +83,9 @@ def make_test_observation_table(
     obs_table.meta["OBSERVATORY_NAME"] = observatory_name
     obs_table.meta["MJDREFI"] = dateref_mjd_int
     obs_table.meta["MJDREFF"] = dateref_mjd_fra
+    obs_table.meta["TIMESYS"] = "TT"
+    obs_table.meta["TIMEUNIT"] = "s"
+    obs_table.meta["TIMEREF"] = "LOCAL"
     if use_abs_time:
         # show the observation times in UTC
         obs_table.meta["TIME_FORMAT"] = "absolute"
@@ -396,6 +404,30 @@ def test_select_sky_regions():
     common_sky_region_select_test_routines(obs_table, selection)
 
 
+def test_create_gti():
+    """Test the `~gammapy.data.ObservationTable.create_gti()` method.
+    """
+    date_start = Time("2012-01-01T00:30:00")
+    date_end = Time("2012-01-01T02:30:00")
+    random_state = np.random.RandomState(seed=0)
+    obs_table = make_test_observation_table(
+        n_obs=1, date_range=(date_start, date_end), random_state=random_state
+    )
+
+    gti = obs_table.create_gti(obs_id=1)
+
+    met_ref = time_ref_from_dict(obs_table.meta)
+    time_start = met_ref + Quantity(obs_table[0]["TSTART"].astype("float64"), "second")
+    time_stop = met_ref + Quantity(obs_table[0]["TSTOP"].astype("float64"), "second")
+
+    assert type(gti) == GTI
+    assert_time_allclose(gti.time_start, time_start)
+    assert_time_allclose(gti.time_stop, time_stop)
+    assert_quantity_allclose(
+        gti.time_sum, Quantity(obs_table[0]["ONTIME"].astype("float64"), "second")
+    )
+
+
 @requires_data("gammapy-extra")
 def test_observation_table_checker():
     path = "$GAMMAPY_EXTRA/datasets/cta-1dc/index/gps/obs-index.fits.gz"

gammapy/data/tests/test_observations.py

@@ -8,11 +8,7 @@
 from ...data import DataStore, EventList, GTI, ObservationCTA
 from ...irf import EffectiveAreaTable2D, EnergyDispersion2D, PSF3D
 from ...utils.testing import requires_data
-from ...utils.testing import (
-    assert_quantity_allclose,
-    assert_time_allclose,
-    assert_skycoord_allclose,
-)
+from ...utils.testing import assert_time_allclose, assert_skycoord_allclose
 
 
 @pytest.fixture(scope="session")
@@ -38,6 +34,26 @@ def test_data_store_observation(data_store):
     assert_skycoord_allclose(obs.target_radec, c)
 
 
+@requires_data("gammapy-extra")
+def test_missing_gti(data_store):
+    """This tests the case when a GTI table is missing in the HDU index file.
+
+    For backward compatibility we create a GTI table on-the-fly by means of the Obs Index file
+    TODO: This method can be removed once GTI tables are explicitly required in Gammapy.
+    """
+    obs = data_store.obs(20136)
+    hdu_bkp = data_store.hdu_table.copy()
+    # remove the row with the GTI file location for obs_id 20136
+    data_store.hdu_table.remove_row(1)
+
+    try:
+        gti = obs.gti
+    finally:
+        data_store.hdu_table = hdu_bkp
+
+    assert type(gti) == GTI
+
+
 @requires_data("gammapy-extra")
 def test_data_store_observation_to_observation_cta(data_store):
     """Test the DataStoreObservation.to_observation_cta conversion method"""