Merge pull request #1173 from lmohrmann/fix-obs-index-times

Fix readout of times from obs-index file

gammapy/data/data_store.py

@@ -9,10 +9,10 @@
 from astropy.table import Table
 from astropy.utils import lazyproperty
 from astropy.units import Quantity
-from astropy.time import Time
 from astropy.coordinates import SkyCoord
 from ..utils.scripts import make_path
 from ..utils.energy import Energy
+from ..utils.time import time_ref_from_dict
 from .obs_table import ObservationTable
 from .hdu_index_table import HDUIndexTable
 from .utils import _earth_location_from_dict
@@ -468,6 +468,13 @@ class DataStoreObservation(object):
     """IACT data store observation.
 
     See :ref:`data_store`
+
+    Parameters
+    ----------
+    obs_id : int
+        Observation ID
+    data_store : `~gammapy.data.DataStore`
+        Data store
     """
 
     def __init__(self, obs_id, data_store):
@@ -480,6 +487,18 @@ def __init__(self, obs_id, data_store):
         self.obs_id = obs_id
         self.data_store = data_store
 
+    def __str__(self):
+        """Generate summary info string."""
+        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.
 
@@ -549,75 +568,35 @@ def bkg(self):
         """Load background object (lazy property)."""
         return self.load(hdu_type='bkg')
 
-    # TODO: maybe the obs table row info should be put in a separate object?
     @lazyproperty
-    def _obs_info(self):
-        """Observation information."""
+    def obs_info(self):
+        """Observation information (`~collections.OrderedDict`)."""
         row = self.data_store.obs_table.select_obs_id(obs_id=self.obs_id)[0]
-        data = OrderedDict(zip(row.colnames, row.as_void()))
-        return data
-
-    @lazyproperty
-    def pointing_radec(self):
-        """Pointing RA / DEC sky coordinates (`~astropy.coordinates.SkyCoord`)."""
-        info = self._obs_info
-        lon, lat = info['RA_PNT'], info['DEC_PNT']
-        return SkyCoord(lon, lat, unit='deg', frame='icrs')
+        return OrderedDict(zip(row.colnames, row.as_void()))
 
     @lazyproperty
     def tstart(self):
         """Observation start time (`~astropy.time.Time`)."""
-        info = self._obs_info
-        return Time(info['TSTART'], format='mjd')
+        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`)."""
-        info = self._obs_info
-        return Time(info['TSTOP'], format='mjd')
-
-    @lazyproperty
-    def muoneff(self):
-        """Observation muon efficiency."""
-        info = self._obs_info
-        return info['MUONEFF']
-
-    @lazyproperty
-    def pointing_altaz(self):
-        """Pointing ALT / AZ sky coordinates (`~astropy.coordinates.SkyCoord`)."""
-        info = self._obs_info
-        alt, az = info['ALT_PNT'], info['AZ_PNT']
-        return SkyCoord(az, alt, unit='deg', frame='altaz')
-
-    @lazyproperty
-    def pointing_zen(self):
-        """Pointing zenith angle sky (`~astropy.units.Quantity`)."""
-        info = self._obs_info
-        zen = info['ZEN_PNT']
-        return Quantity(zen, unit='deg')
-        return Quantity(info['TSTART'], 'second')
-
-    @lazyproperty
-    def target_radec(self):
-        """Target RA / DEC sky coordinates (`~astropy.coordinates.SkyCoord`)."""
-        info = self._obs_info
-        lon, lat = info['RA_OBJ'], info['DEC_OBJ']
-        return SkyCoord(lon, lat, unit='deg', frame='icrs')
-
-    @lazyproperty
-    def observatory_earth_location(self):
-        """Observatory location (`~astropy.coordinates.EarthLocation`)."""
-        info = self._obs_info
-        return _earth_location_from_dict(info)
+        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.
         """
-        info = self._obs_info
-        return Quantity(info['ONTIME'], 'second')
+        return Quantity(self.obs_info['ONTIME'], 'second')
 
     @lazyproperty
     def observation_live_time_duration(self):
@@ -628,8 +607,7 @@ def observation_live_time_duration(self):
         Computed as ``t_live = t_observation * (1 - f_dead)``
         where ``f_dead`` is the dead-time fraction.
         """
-        info = self._obs_info
-        return Quantity(info['LIVETIME'], 'second')
+        return Quantity(self.obs_info['LIVETIME'], 'second')
 
     @lazyproperty
     def observation_dead_time_fraction(self):
@@ -645,30 +623,43 @@ def observation_dead_time_fraction(self):
         The dead-time fraction is used in the live-time computation,
         which in turn is used in the exposure and flux computation.
         """
-        info = self._obs_info
-        return 1 - info['DEADC']
+        return 1 - self.obs_info['DEADC']
 
-    def __str__(self):
-        """Generate summary info string."""
-        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)
+    @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')
 
-        # TODO: Which target was observed?
-        # TODO: print info about available HDUs for this observation ...
-        return ss
+    @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')
 
-    def peek(self):
-        """Quick-look plots in a few panels."""
-        raise NotImplementedError
+    @lazyproperty
+    def pointing_zen(self):
+        """Pointing zenith angle sky (`~astropy.units.Quantity`)."""
+        return Quantity(self.obs_info['ZEN_PNT'], unit='deg')
 
-    def make_exposure_image(self, fov, energy_range):
-        """Make exposure image.
+    @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')
 
-        TODO: Do we want such methods here or as standalone functions that work with obs objects?
-        """
+    @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 make_psf(self, position, energy=None, rad=None):

gammapy/data/tests/test_data_store.py

@@ -7,8 +7,10 @@
 from astropy.coordinates import Angle, SkyCoord
 from astropy.units import Quantity
 import astropy.units as u
+from astropy.time import Time
 from ...data import DataStore, DataManager, ObservationList
 from ...utils.testing import data_manager, requires_data, requires_dependency
+from ...utils.testing import assert_time_allclose, assert_skycoord_allclose
 from ...utils.energy import Energy
 from ...datasets import gammapy_extra
 from ...utils.energy import EnergyBounds
@@ -75,6 +77,7 @@ def test_datastore_load_all(data_manager):
     assert_allclose(event_lists[0].table['ENERGY'][0], 1.1156039)
     assert_allclose(event_lists[-1].table['ENERGY'][0], 1.0204216)
 
+
 @requires_data('gammapy-extra')
 @requires_dependency('yaml')
 def test_datastore_obslist(data_manager):
@@ -89,6 +92,7 @@ def test_datastore_obslist(data_manager):
     obslist = data_store.obs_list([11111, 23523], skip_missing=True)
     assert obslist[0].obs_id == 23523
 
+
 @requires_data('gammapy-extra')
 @requires_dependency('yaml')
 def test_datastore_subset(tmpdir, data_manager):
@@ -132,6 +136,27 @@ def test_data_summary(data_manager):
     assert t[0]['psf_3gauss'] == 6042
 
 
+
+@requires_data('gammapy-extra')
+@requires_dependency('yaml')
+def test_data_store_observation():
+    """Test DataStoreObservation class"""
+    data_store = DataStore.from_dir('$GAMMAPY_EXTRA/datasets/hess-crab4-hd-hap-prod2/')
+    obs = data_store.obs(23523)
+
+    assert_time_allclose(obs.tstart, Time(51545.11740650318, scale='tt', format='mjd'))
+    assert_time_allclose(obs.tstop, Time(51545.11740672924, scale='tt', format='mjd'))
+
+    c = SkyCoord(83.63333129882812, 21.51444435119629, unit='deg')
+    assert_skycoord_allclose(obs.pointing_radec, c)
+
+    c = SkyCoord(26.533863067626953, 40.60616683959961, unit='deg')
+    assert_skycoord_allclose(obs.pointing_altaz, c)
+
+    c = SkyCoord(83.63333129882812, 22.01444435119629, unit='deg')
+    assert_skycoord_allclose(obs.target_radec, c)
+
+
 @requires_dependency('scipy')
 @requires_data('gammapy-extra')
 @pytest.mark.parametrize("pars,result", [
@@ -179,7 +204,7 @@ def test_make_psf(pars, result):
 
 @requires_dependency('scipy')
 @requires_data('gammapy-extra')
-def test_make_mean_edisp(tmpdir):
+def test_make_mean_edisp():
     position = SkyCoord(83.63, 22.01, unit='deg')
     store = gammapy_extra.filename("datasets/hess-crab4-hd-hap-prod2")
     data_store = DataStore.from_dir(store)

gammapy/spectrum/extract.py

@@ -148,7 +148,7 @@ def make_empty_vectors(self, obs, bkg):
         """
         log.info('Update observation meta info')
         # Copy over existing meta information
-        meta = OrderedDict(obs._obs_info)
+        meta = OrderedDict(obs.obs_info)
         offset = obs.pointing_radec.separation(bkg.on_region.center)
         log.info('Offset : {}\n'.format(offset))
         meta['OFFSET'] = offset.deg

gammapy/utils/testing.py

@@ -4,7 +4,8 @@
 
 import os
 import pytest
-from astropy.coordinates import Angle
+from astropy.coordinates import Angle, SkyCoord
+from astropy.time import Time
 from numpy.testing import assert_array_less, assert_allclose
 
 from ..data import DataManager
@@ -15,6 +16,7 @@
     'requires_data',
     'assert_wcs_allclose',
     'assert_skycoord_allclose',
+    'assert_time_allclose',
 ]
 
 # Cache for `requires_dependency`
@@ -134,12 +136,21 @@ def assert_wcs_allclose(wcs1, wcs2):
     assert_allclose(wcs1.wcs.cdelt, wcs2.wcs.cdelt)
 
 
-def assert_skycoord_allclose(skycoord1, skycoord2, atol='1 arcsec'):
+def assert_skycoord_allclose(actual, desired):
     """Assert all-close for `~astropy.coordinates.SkyCoord`.
 
-    - Checks if separation on the sky is within ``atol``.
     - Frames can be different, aren't checked at the moment.
     """
-    atol = Angle(atol)
-    sep = skycoord1.separation(skycoord2).deg
-    assert_array_less(sep.deg, atol.deg)
+    assert isinstance(actual, SkyCoord)
+    assert isinstance(desired, SkyCoord)
+    assert_allclose(actual.data.lon.value, desired.data.lon.value)
+    assert_allclose(actual.data.lat.value, desired.data.lat.value)
+
+
+def assert_time_allclose(actual, desired):
+    """Assert that two `astropy.time.Time` objects are almost the same."""
+    assert isinstance(actual, Time)
+    assert isinstance(desired, Time)
+    assert_allclose(actual.value, desired.value)
+    assert actual.scale == desired.scale
+    assert actual.format == desired.format