diff --git a/docs/dataformats/file_formats.rst b/docs/dataformats/file_formats.rst index 58ff09e7b4..57d21f9ba2 100644 --- a/docs/dataformats/file_formats.rst +++ b/docs/dataformats/file_formats.rst @@ -145,6 +145,8 @@ INI files ... to use it in your code Unfortunately INI files are not standardised, so there's only conventions and tons of variants. +.. _CSV_files: + CSV +++ @@ -194,6 +196,9 @@ XML **XML** files (see `Wikipedia `__) +GammaLib / ctools uses an "observation definition" XML format described +`here `__. + TODO: describe FITS diff --git a/docs/dataformats/index.rst b/docs/dataformats/index.rst index 57c9b15426..5b03bd490f 100644 --- a/docs/dataformats/index.rst +++ b/docs/dataformats/index.rst @@ -19,6 +19,21 @@ we simply use general-purpose file formats define our own semantics as described in the :ref:`dataformats_file_formats` section. +.. _time_gammapy: + +*************** +Time in Gammapy +*************** + +Gammapy adopts the exact same time standard as the Fermi Science Tools as described on +this page: `Cicerone: Data - Time in Fermi Data Analysis +`_. +Every time should be defined as a `~astropy.time.Time` object, representing UTC seconds +after the "Mission elapsed time" [MET]_ origin. Each experiment may have a different +[MET]_ origin that should be included in the header of the data files. + + .. toctree:: :maxdepth: 1 diff --git a/docs/dataformats/observation_lists.rst b/docs/dataformats/observation_lists.rst index a5cb0b3239..02d593f96e 100644 --- a/docs/dataformats/observation_lists.rst +++ b/docs/dataformats/observation_lists.rst @@ -10,38 +10,90 @@ the information in the run list. These filenames can either be given directly or there can be naming conventions, e.g. the event list for run ``42`` could be stored at ``$TEV_DATA/RUN_000042/Events.fits``. +A run list must have at least a column called ``OBS_ID``:: + + OBS_ID + 42 + 43 -CSV format ----------- +Usually it has many more columns with information about each observation. A list of +Gammapy supported columns is: -We use the `CSV `_ (comma-separated values) format for run lists. -This has the advantage that everyone can work with whatever tool they like. Here's some good options: +================ =========================================================================== ========= +column name description required? +================ =========================================================================== ========= +OBS_ID observation ID as an integer yes +RA pointing position right ascension in equatorial (ICRS) coordinates yes? +DEC pointing position declination in equatorial (ICRS) coordinates yes? +AZ average azimuth angle during the observation no +ALT average altitude angle during the observation no +MUON_EFFICIENCY average muon efficiency of the telescopes yes? +TIME_START start time of the observation stored as number of seconds after [MET]_ yes? +TIME_STOP end time of the observation in the same format as TIME_START no +TIME_OBSERVATION duration of the observation no +TIME_LIVE duration of the observation without dead time no +TRIGGER_RATE average trigger rate of the system no +MEAN_TEMPERATURE average temperature of the atmosphere no +N_TELS number of telescopes participating in the observation yes? +TEL_LIST string with a [CSV]_ list of telescope IDs participating in the observation yes? +QUALITY data quality; recommended: "spectral" or "detection" no +================ =========================================================================== ========= -* `LibreOffice Calc `_ -* `TOPCAT `_ or `STILTS `_ -* `csvkit `_ -* Your own script using e.g the `Python standard library CSV module `_ or `pandas `_ +Extra user-defined columns are allowed; Gammapy will ignore them. -A run list must have at least a column called ``Run``:: - - Run - 42 - 43 +In order for the extra columns to have full meaning the following is needed: -Usually it has many more columns with information about each run:: - - Run,Telescope_Pattern,Date,Duration,GLON,GLAT,Target,Zenith - 1234,24,2013-03-22 14:32,1832,83.7,-5.2,Crab Nebula,32 + * Extra row right after the column name, specifying the unit of the quantity listed on each column. + * A header with at least the following keywords: + +================ =========================================================================== ========= +keyword description required? +================ =========================================================================== ========= +OBSERVATORY_NAME name of the observatory where the observations were taken. This is no + important for instance for coordinate transformations between celestial + (i.e. RA/dec) and terrestrial (i.e. az/alt) coordinate systems. +MJDREFI reference time: integer value in mean julian days; details in yes? + :ref:`time_gammapy`. +MJDREFF reference time: fraction of integer value defined in MJDREFI; details in yes? + :ref:`time_gammapy`. +================ =========================================================================== ========= + +Extra user-defined header entries are allowed; Gammapy will ignore them. + + +Example +------- +The tool `~gammapy.datasets.make_test_observation_table` can generate a `~gammapy.obs.ObservationTable` +with dummy values. -Special column names that the Gammapy analysis tools understand: +Header (metadata):: -* ``Run`` --- Run number (int) -* ``Telescope_Pattern`` --- Binary pattern describing which telescopes participated in the run -* ``Date`` --- Date of observation -* ``RA``, ``DEC`` or ``GLON``, ``GLAT`` -- Pointing position in Equatorial (ICRS) or Galactic coordinates + {u'MJDREFI': 55197.0, u'OBSERVATORY_NAME': u'HESS', u'MJDREFF': 0.0} -XML format ----------- +Table: -GammaLib / ctools uses an "observation definition" XML format described -`here `__. ++------+----------------+---------+-------------+-------------+-------------+-------------+-------------+--------------+------+---------------+ +|OBS_ID|TIME_OBSERVATION|TIME_LIVE| TIME_START | TIME_STOP | AZ | ALT | RA | DEC |N_TELS|MUON_EFFICIENCY| ++------+----------------+---------+-------------+-------------+-------------+-------------+-------------+--------------+------+---------------+ +| | s | s | s | s | deg | deg | deg | deg | | | ++======+================+=========+=============+=============+=============+=============+=============+==============+======+===============+ +| 1| 1800.0| 1500.0|118890428.352|118892228.352|130.926874751|49.6209457026|96.3849089136|-43.6914197077| 3| 0.814535992712| ++------+----------------+---------+-------------+-------------+-------------+-------------+-------------+--------------+------+---------------+ +| 2| 1800.0| 1500.0|112242990.559|112244790.559|272.213179564|70.2673929472| 339.00128923|-21.1698098192| 3| 0.976469816749| ++------+----------------+---------+-------------+-------------+-------------+-------------+-------------+--------------+------+---------------+ +| 3| 1800.0| 1500.0| 66444741.854| 66446541.854|346.848818939| 46.138583188|162.086175054| 19.6398873974| 4| 0.920096961383| ++------+----------------+---------+-------------+-------------+-------------+-------------+-------------+--------------+------+---------------+ +| 4| 1800.0| 1500.0|22388716.9244|22390516.9244|300.850748052|55.1330124055|32.9474858892|-3.19910057294| 3| 0.678431411337| ++------+----------------+---------+-------------+-------------+-------------+-------------+-------------+--------------+------+---------------+ +| 5| 1800.0| 1500.0|135469063.888|135470863.888|355.160931662| 48.734744852|197.123663537| 17.9411145072| 4| 0.77879533822| ++------+----------------+---------+-------------+-------------+-------------+-------------+-------------+--------------+------+---------------+ +| 6| 1800.0| 1500.0|21857681.6916|21859481.6916|124.846967209| 78.859585347| 14.162859563|-29.3419432185| 4| 0.709642622408| ++------+----------------+---------+-------------+-------------+-------------+-------------+-------------+--------------+------+---------------+ +| 7| 1800.0| 1500.0| 57554741.106| 57556541.106|268.541714486|48.8489560299|64.8265458802|-18.2634404823| 3| 0.908426763354| ++------+----------------+---------+-------------+-------------+-------------+-------------+-------------+--------------+------+---------------+ +| 8| 1800.0| 1500.0|19181027.9045|19182827.9045|120.558129392| 49.663761361| 24.791511978|-37.1789681608| 4| 0.980162662473| ++------+----------------+---------+-------------+-------------+-------------+-------------+-------------+--------------+------+---------------+ +| 9| 1800.0| 1500.0|120447694.583|120449494.583| 132.10271454|78.7455993174|89.7950895353|-30.5128854184| 3| 0.807695978946| ++------+----------------+---------+-------------+-------------+-------------+-------------+-------------+--------------+------+---------------+ +| 10| 1800.0| 1500.0|144207430.361|144209230.361|323.562657045|45.4005803262|324.596045439| 13.6761217326| 3| 0.694201696626| ++------+----------------+---------+-------------+-------------+-------------+-------------+-------------+--------------+------+---------------+ diff --git a/docs/references.rst b/docs/references.rst index 64cbf273d9..1999b7cc3f 100644 --- a/docs/references.rst +++ b/docs/references.rst @@ -3,6 +3,11 @@ References ========== +.. _publications: + +Publications +------------ + This is the bibliography containing the literature references for the implemented methods referenced from the Gammapy docs. @@ -58,3 +63,12 @@ Software references: .. [FSSC2013] `Fermi LAT Collaboration (2013) `_ "Science Tools: LAT Data Analysis Tools" + + +.. _glossary: + +Glossary +-------- + +.. [CSV] comma-separated values; see also :ref:`CSV_files`. +.. [MET] mission elapsed time; see also :ref:`time_gammapy`. diff --git a/docs/utils/index.rst b/docs/utils/index.rst index 1fab483966..d9c37ba28c 100644 --- a/docs/utils/index.rst +++ b/docs/utils/index.rst @@ -56,4 +56,6 @@ Reference/API .. automodapi:: gammapy.utils.testing :no-inheritance-diagram: - \ No newline at end of file + +.. automodapi:: gammapy.utils.time + :no-inheritance-diagram: diff --git a/gammapy/data/event_list.py b/gammapy/data/event_list.py index d5c46a0d27..56766ab060 100644 --- a/gammapy/data/event_list.py +++ b/gammapy/data/event_list.py @@ -13,6 +13,7 @@ from ..image import wcs_histogram2d from ..data import GoodTimeIntervals, TelescopeArray from ..data import InvalidDataError +from ..utils.time import time_ref_from_dict from . import utils __all__ = ['EventList', @@ -78,7 +79,7 @@ def time(self): With 32-bit floats times will be incorrect by a few seconds when e.g. adding them to the reference time. """ - met_ref = utils._time_ref_from_dict(self.meta) + met_ref = time_ref_from_dict(self.meta) met = Quantity(self['TIME'].astype('float64'), 'second') time = met_ref + met return time @@ -619,7 +620,7 @@ def check_times(self): ) telescope = self.dset.event_list.meta['TELESCOP'] - met_ref = utils._time_ref_from_dict(self.dset.event_list.meta) + met_ref = time_ref_from_dict(self.dset.event_list.meta) if telescope in telescope_met_refs.keys(): dt = (met_ref - telescope_met_refs[telescope]) diff --git a/gammapy/data/gti.py b/gammapy/data/gti.py index 810053d55c..d6a4b2a437 100644 --- a/gammapy/data/gti.py +++ b/gammapy/data/gti.py @@ -4,6 +4,7 @@ import numpy as np from astropy.units import Quantity from astropy.table import Table +from ..utils.time import time_ref_from_dict from . import utils __all__ = ['GoodTimeIntervals'] @@ -44,13 +45,13 @@ def time_sum(self): @property def time_start(self): """GTI start times (`~astropy.time.Time`).""" - met_ref = utils._time_ref_from_dict(self.meta) + met_ref = time_ref_from_dict(self.meta) met = Quantity(self['START'].astype('float64'), 'second') return met_ref + met @property def time_stop(self): """GTI end times (`~astropy.time.Time`).""" - met_ref = utils._time_ref_from_dict(self.meta) + met_ref = time_ref_from_dict(self.meta) met = Quantity(self['STOP'].astype('float64'), 'second') return met_ref + met diff --git a/gammapy/data/utils.py b/gammapy/data/utils.py index cf8835a1c4..ef46c18359 100644 --- a/gammapy/data/utils.py +++ b/gammapy/data/utils.py @@ -2,7 +2,6 @@ """Misc utility functions.""" from __future__ import (absolute_import, division, print_function, unicode_literals) -from astropy.time import Time from astropy.coordinates import Angle, EarthLocation from astropy.units import Quantity @@ -12,13 +11,3 @@ def _earth_location_from_dict(meta): lat = Angle(meta['GEOLAT'], 'deg') height = Quantity(meta['ALTITUDE'], 'meter') return EarthLocation(lon=lon, lat=lat, height=height) - - -def _time_ref_from_dict(meta): - # Note: the `float` call here is to make sure we use 64-bit - mjd = float(meta['MJDREFI']) + float(meta['MJDREFF']) - # TODO: Is 'tt' a default we should put here? - scale = meta.get('TIMESYS', 'tt').lower() - # Note: we could call .copy('iso') or .replicate('iso') - # here if we prefer 'iso' over 'mjd' format in most places. - return Time(mjd, format='mjd', scale=scale) diff --git a/gammapy/datasets/make.py b/gammapy/datasets/make.py index 536303323c..0293e0900c 100644 --- a/gammapy/datasets/make.py +++ b/gammapy/datasets/make.py @@ -5,9 +5,17 @@ unicode_literals) import numpy as np from astropy.units import Quantity +from astropy.time import Time, TimeDelta +from astropy.table import Table +from astropy.coordinates import SkyCoord, AltAz, FK5, Angle from ..irf import EnergyDependentMultiGaussPSF +from ..obs import ObservationTable, observatory_locations +from ..utils.random import sample_sphere +from ..utils.time import time_ref_from_dict, time_relative_to_ref -__all__ = ['make_test_psf'] +__all__ = ['make_test_psf', + 'make_test_observation_table', + ] def make_test_psf(energy_bins=15, theta_bins=12): @@ -64,3 +72,144 @@ def sigma_energy_theta(energy, theta, sigma): zenith=zenith_hi) return psf + + +def make_test_observation_table(observatory_name, n_obs, debug=False): + """Make a test observation table. + + For the moment, only random observation tables are created. + + Parameters + ---------- + observatory_name : string + name of the observatory; a list of choices is given in `~gammapy.obs.observatory_locations` + n_obs : integer + number of observations for the obs table + debug : bool + show UTC times instead of seconds after the reference + + Returns + ------- + obs_table : `~gammapy.obs.ObservationTable` + observation table + """ + n_obs_start = 1 + + obs_table = ObservationTable() + + # build a time reference as the start of 2010 + dateref = Time('2010-01-01 00:00:00', format='iso', scale='utc') + dateref_mjd_fra, dateref_mjd_int = np.modf(dateref.mjd) + + # header + header = {'OBSERVATORY_NAME': observatory_name, + 'MJDREFI': dateref_mjd_int, 'MJDREFF': dateref_mjd_fra} + ObservationTable.meta = header + + # obs id + obs_id = np.arange(n_obs_start, n_obs_start + n_obs) + obs_table['OBS_ID'] = obs_id + + # obs time: 30 min + time_observation = Quantity( + 30. * np.ones_like(obs_id), 'minute').to('second') + obs_table['TIME_OBSERVATION'] = time_observation + + # livetime: 25 min + time_live = Quantity(25. * np.ones_like(obs_id), 'minute').to('second') + obs_table['TIME_LIVE'] = time_live + + # start time + # random points between the start of 2010 and the end of 2014 + # using the start of 2010 as a reference time for the header of the table + # observations restrict to night time + # considering start of astronomical day at midday: implicit in setting the + # start of the night, when generating random night hours + datestart = Time('2010-01-01 00:00:00', format='iso', scale='utc') + dateend = Time('2015-01-01 00:00:00', format='iso', scale='utc') + time_start = Time((dateend.mjd - datestart.mjd) * + np.random.random(len(obs_id)) + datestart.mjd, format='mjd', scale='utc') + + # keep only the integer part (i.e. the day, not the fraction of the day) + time_start_f, time_start_i = np.modf(time_start.mjd) + time_start = Time(time_start_i, format='mjd', scale='utc') + + # random generation of night hours: 6 h (from 22 h to 4 h), leaving 1/2 h + # time for the last run to finish + night_start = Quantity(22., 'hour') + night_duration = Quantity(5.5, 'hour') + hour_start = night_start + night_duration * np.random.random(len(obs_id)) + + # add night hour to integer part of MJD + time_start += hour_start + + if debug: + # show the observation times in UTC + time_start = time_start.iso + else: + # show the observation times in seconds after the reference + time_start = time_relative_to_ref(time_start, header) + # converting to quantity (beter treatment of units) + time_start = Quantity(time_start.sec, 'second') + + obs_table['TIME_START'] = time_start + + # stop time + # calculated as TIME_START + TIME_OBSERVATION + if debug: + time_stop = Time(obs_table['TIME_START']) + \ + TimeDelta(obs_table['TIME_OBSERVATION']) + else: + time_stop = TimeDelta( + obs_table['TIME_START']) + TimeDelta(obs_table['TIME_OBSERVATION']) + # converting to quantity (beter treatment of units) + time_stop = Quantity(time_stop.sec, 'second') + + obs_table['TIME_STOP'] = time_stop + + # az, alt + # random points in a sphere above 45 deg altitude + az, alt = sample_sphere(len(obs_id), (0, 360), (45, 90), 'degree') + az = Angle(az, 'degree') + alt = Angle(alt, 'degree') + obs_table['AZ'] = az + obs_table['ALT'] = alt + + # RA, dec + # derive from az, alt taking into account that alt, az represent the values + # at the middle of the observation, i.e. at time_ref + (TIME_START + TIME_STOP)/2 + # (or better: time_ref + TIME_START + (TIME_OBSERVATION/2)) + # in debug modus, the time_ref should not be added, since it's already included + # in TIME_START and TIME_STOP + az = Angle(obs_table['AZ']) + alt = Angle(obs_table['ALT']) + if debug: + obstime = Time(obs_table['TIME_START']) + \ + TimeDelta(obs_table['TIME_OBSERVATION']) / 2. + else: + obstime = time_ref_from_dict(obs_table.meta) + TimeDelta( + obs_table['TIME_START']) + TimeDelta(obs_table['TIME_OBSERVATION']) / 2. + location = observatory_locations[observatory_name] + alt_az_coord = AltAz(az=az, alt=alt, obstime=obstime, location=location) + sky_coord = alt_az_coord.transform_to(FK5) + obs_table['RA'] = sky_coord.ra + obs_table['DEC'] = sky_coord.dec + + # positions + + # number of telescopes + # random integers between 3 and 4 + n_tels_min = 3 + n_tels_max = 4 + n_tels = np.random.randint(n_tels_min, n_tels_max + 1, len(obs_id)) + obs_table['N_TELS'] = n_tels + + # muon efficiency + # random between 0.6 and 1.0 + muon_efficiency_min = 0.6 + muon_efficiency_max = 1.0 + muon_efficiency = np.random.random( + len(obs_id)) * (muon_efficiency_max - muon_efficiency_min) + muon_efficiency_min + obs_table['MUON_EFFICIENCY'] = muon_efficiency + + return obs_table diff --git a/gammapy/datasets/tests/test_make.py b/gammapy/datasets/tests/test_make.py index 41f116d522..447124138c 100644 --- a/gammapy/datasets/tests/test_make.py +++ b/gammapy/datasets/tests/test_make.py @@ -4,7 +4,8 @@ from astropy.coordinates import Angle from astropy.units import Quantity from ...utils.testing import assert_quantity -from ...datasets import make_test_psf +from ...datasets import make_test_psf, make_test_observation_table +from ...obs import ObservationTable def test_make_test_psf_fits_table(): @@ -15,3 +16,30 @@ def test_make_test_psf_fits_table(): psf2 = psf.psf_at_energy_and_theta(energy, theta) fraction = psf2.containment_fraction(0.1) assert fraction == 0.44485638998490795 + + +def test_make_test_observation_table(): + observatory_name='HESS' + n_obs = 10 + obs_table = make_test_observation_table(observatory_name, n_obs) + + # test: assert if the length of the table is n_obs: + assert len(obs_table) == n_obs + + # test: assert if the TIME_START > 0: + assert (obs_table['TIME_START'] > 0).all() + + # test: assert if TIME_STOP > TIME_START: + assert (obs_table['TIME_STOP'] > obs_table['TIME_START']).all() + + # test: assert if RA is in the interval (0, 360) deg: + ra_min = Angle(0, 'degree') + ra_max = Angle(360, 'degree') + assert (ra_min < obs_table['RA']).all() + assert (obs_table['RA'] < ra_max).all() + + # test: assert if dec is inthe interval (-90, 90) deg: + dec_min = Angle(-90, 'degree') + dec_max = Angle(90, 'degree') + assert (dec_min < obs_table['DEC']).all() + assert (obs_table['DEC'] < dec_max).all() diff --git a/gammapy/obs/observation.py b/gammapy/obs/observation.py index 23d102f462..97ac78d961 100644 --- a/gammapy/obs/observation.py +++ b/gammapy/obs/observation.py @@ -3,6 +3,9 @@ unicode_literals) import numpy as np from astropy.table import Table +from astropy.units import Quantity +from astropy.time import Time +from ..utils.time import time_ref_from_dict __all__ = [ # 'Observation', @@ -11,8 +14,9 @@ class Observation(object): - """Observation (a.k.a. Run). + """Observation. + An observation is a.k.a. run. TODO: not clear if this class is useful. Parameters @@ -38,24 +42,24 @@ class ObservationTable(Table): """Observation table (a.k.a. run list). This is an `~astropy.table.Table` sub-class, with a few - convenience methods and the following columns: - - * ``OBS_ID`` - * ``ONTIME`` - * ``LIVETIME`` - * ... - + convenience methods. The format of the observation table + is described in :ref:`dataformats_observation_lists`. """ def info(self): ss = 'Observation table:\n' + obs_name = self.meta['OBSERVATORY_NAME'] + ss += 'Observatory name: {}\n'.format(obs_name) ss += 'Number of observations: {}\n'.format(len(self)) - ontime = self['ONTIME'].sum() + ontime = Quantity(self['TIME_OBSERVATION'].sum(), self['TIME_OBSERVATION'].unit) ss += 'Total observation time: {}\n'.format(ontime) - livetime = self['LIVETIME'].sum() + livetime = Quantity(self['TIME_LIVE'].sum(), self['TIME_LIVE'].unit) ss += 'Total live time: {}\n'.format(livetime) dtf = 100. * (1 - livetime / ontime) - ss += 'Average dead time fraction: {:5.2f}%'.format(dtf) + ss += 'Average dead time fraction: {:5.2f}%\n'.format(dtf) + time_ref = time_ref_from_dict(self.meta) + time_ref_unit = time_ref_from_dict(self.meta).format + ss += 'Time reference: {} {}'.format(time_ref, time_ref_unit) return ss def select_linspace_subset(self, num): @@ -82,4 +86,3 @@ def select_linspace_subset(self, num): # Round down to nearest integer indices = indices.astype('int') return self[indices] - diff --git a/gammapy/utils/random.py b/gammapy/utils/random.py index 672796c49b..2f0c523cba 100644 --- a/gammapy/utils/random.py +++ b/gammapy/utils/random.py @@ -30,7 +30,7 @@ def check_random_state(seed): ' instance'.format(seed)) -def sample_sphere(size, lon_range=None, lat_range=None, unit='radians'): +def sample_sphere(size, lon_range=None, lat_range=None, unit='rad'): """Sample random points on the sphere. Reference: http://mathworld.wolfram.com/SpherePointPicking.html @@ -42,7 +42,7 @@ def sample_sphere(size, lon_range=None, lat_range=None, unit='radians'): lon_range : tuple Longitude range (min, max) tuple in range (0, 360) lat_range : tuple - Latitude range (min, max) tuple in range (-180, 180) + Latitude range (min, max) tuple in range (-90, 90) units : {'rad', 'deg'} Units of input range and returned angles @@ -51,6 +51,14 @@ def sample_sphere(size, lon_range=None, lat_range=None, unit='radians'): lon, lat: arrays Longitude and latitude coordinate arrays """ + # Correctly interpret units at an early stage + if unit in ['rad', 'radian', 'radians']: + unit = 'rad' + elif unit in ['deg', 'degree', 'degrees']: + unit = 'deg' + else: + raise ValueError('Invalid unit: {0}'.format(unit)) + # Convert inputs to internal format (all radians) size = int(size) @@ -76,12 +84,10 @@ def sample_sphere(size, lon_range=None, lat_range=None, unit='radians'): lat = np.arcsin(z) # Return result - if unit in ['rad', 'radian', 'radians']: + if unit == 'rad': return lon, lat - elif unit in ['deg', 'degree', 'degrees']: + elif unit == 'deg': return np.degrees(lon), np.degrees(lat) - else: - raise ValueError('Invalid unit: {0}'.format(unit)) def sample_powerlaw(x_min, x_max, gamma, size=None): diff --git a/gammapy/utils/tests/test_time.py b/gammapy/utils/tests/test_time.py new file mode 100644 index 0000000000..bc3a9d4cd0 --- /dev/null +++ b/gammapy/utils/tests/test_time.py @@ -0,0 +1,23 @@ +# 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_almost_equal +from astropy.time import Time, TimeDelta +from ..time import time_ref_from_dict, time_relative_to_ref + + +def test_time_ref_from_dict(): + mjd_int = 500 + mjd_frac = 0.5 + time_ref_dict = dict(MJDREFI=mjd_int, MJDREFF=mjd_frac) + time_ref = time_ref_from_dict(time_ref_dict) + assert_almost_equal(time_ref.mjd, mjd_int + mjd_frac, decimal=4) + + +def test_time_relative_to_ref(): + time_ref_dict = dict(MJDREFI=500, MJDREFF=0.5) + time_ref = time_ref_from_dict(time_ref_dict) + delta_time_1sec = TimeDelta(1., format='sec') + time = time_ref + delta_time_1sec + delta_time = time_relative_to_ref(time, time_ref_dict) + assert_almost_equal(delta_time.sec, delta_time_1sec.sec, decimal=4) diff --git a/gammapy/utils/time.py b/gammapy/utils/time.py new file mode 100644 index 0000000000..0b0c12ee31 --- /dev/null +++ b/gammapy/utils/time.py @@ -0,0 +1,59 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +"""Time related utility functions.""" +from __future__ import (absolute_import, division, print_function, + unicode_literals) +from astropy.time import Time, TimeDelta + +__all__ = ['time_ref_from_dict', + 'time_relative_to_ref', + ] + + +def time_ref_from_dict(meta): + """Calculate the time reference from metadata. + + The time reference is built as MJDREFI + MJDREFF in units of MJD. + All other times should be interpreted as seconds after the reference. + + Parameters + ---------- + meta : `~dict` + dictionary with the keywords `MJDREFI` and `MJDREFF` + + Returns + ------- + time : `~astropy.time.Time` + reference time with ``format='MJD'`` + """ + # Note: the `float` call here is to make sure we use 64-bit + mjd = float(meta['MJDREFI']) + float(meta['MJDREFF']) + # TODO: Is 'tt' a default we should put here? + scale = meta.get('TIMESYS', 'tt').lower() + # Note: we could call .copy('iso') or .replicate('iso') + # here if we prefer 'iso' over 'mjd' format in most places. + + return Time(mjd, format='mjd', scale=scale) + + +def time_relative_to_ref(time, meta): + """Convert a time using an existing reference. + + The time reference is built as MJDREFI + MJDREFF in units of MJD. + The time will be converted to seconds after the reference. + + Parameters + ---------- + time: `~astropy.time.Time` + time to be converted. + meta : dict + dictionary with the keywords `MJDREFI` and `MJDREFF` + + Returns + ------- + time_delta : `~astropy.time.TimeDelta` + time in seconds after the reference + """ + time_ref = time_ref_from_dict(meta) + delta_time = TimeDelta(time - time_ref, format='sec') + + return delta_time