From 483fc1866c026529f3399dccaa3e35765f176206 Mon Sep 17 00:00:00 2001 From: mapaz Date: Tue, 2 Jun 2015 16:13:49 -0700 Subject: [PATCH 01/15] Adding ObservationTable generator. --- gammapy/datasets/make.py | 95 ++++++++++++++++++++++++++++- gammapy/datasets/tests/test_make.py | 9 ++- gammapy/obs/observation.py | 6 +- 3 files changed, 107 insertions(+), 3 deletions(-) diff --git a/gammapy/datasets/make.py b/gammapy/datasets/make.py index 536303323c..220bebce6a 100644 --- a/gammapy/datasets/make.py +++ b/gammapy/datasets/make.py @@ -5,9 +5,15 @@ unicode_literals) import numpy as np from astropy.units import Quantity +from astropy.time import Time +from astropy.table import Table, Column +from astropy.coordinates import SkyCoord, AltAz, Angle from ..irf import EnergyDependentMultiGaussPSF +from ..obs import ObservationTable, observatory_locations -__all__ = ['make_test_psf'] +__all__ = ['make_test_psf', + 'generate_observation_table', + ] def make_test_psf(energy_bins=15, theta_bins=12): @@ -64,3 +70,90 @@ def sigma_energy_theta(energy, theta, sigma): zenith=zenith_hi) return psf + + +def generate_observation_table(observatory, n_obs): + """Generate an observation table. + + For the moment, only random observation tables are created. + + Parameters + ---------- + observatory : string + name of the observatory + n_obs : integer + number of observations for the obs table + + Returns + ------- + obs_table : `~gammapy.obs.ObservationTable` + observation table + """ + n_obs_start = 1 + + astro_table = Table() + + # obs id + obs_id = np.arange(n_obs_start, n_obs_start + n_obs) + col_obs_id = Column(name='OBS_ID', data=obs_id) + astro_table.add_column(col_obs_id) + + # on time + ontime = Quantity(30.*np.ones_like(col_obs_id.data), 'minute') + col_ontime = Column(name='ONTIME', data=ontime) + astro_table.add_column(col_ontime) + + # livetime + livetime = Quantity(25.*np.ones_like(col_obs_id.data), 'minute') + col_livetime = Column(name='LIVETIME', data=livetime) + astro_table.add_column(col_livetime) + + # TODO: add columns for coordinates: + # pointing observation -> done + # and date -> done + # then convert to alt az, considering the observatory location + # TODO: add column(s) for offset, and take it into account in the coord. transformation!!! + + # RA, Dec + # random points on a sphere ref: http://mathworld.wolfram.com/SpherePointPicking.html + + ra = Angle(360.*np.random.random(len(col_obs_id)), 'degree') + col_ra = Column(name='RA', data=ra) + astro_table.add_column(col_ra) + + dec = Angle(np.arccos(2.*np.random.random(len(col_obs_id)) - 1), 'radian').to('degree') + # translate angles from [0, 180) deg to [-90, 90) deg + dec = dec - Angle(90., 'degree') + col_dec = Column(name='DEC', data=dec) + astro_table.add_column(col_dec) + + # date + # random points between the start of 2010 and the end of 2014 + #TODO: should this represent the time at the beginning of the run? + # this has consequences for the ra/dec -> alt/az conversion + datestart = Time('2010-01-01T00:00:00', format='fits', scale='utc') + dateend = Time('2015-01-01T00:00:00', format='fits', scale='utc') + date = Time((dateend.mjd - datestart.mjd)*np.random.random(len(col_obs_id)) + datestart.mjd, format='mjd', scale='utc').fits + col_date = Column(name='DATE', data=date) + astro_table.add_column(col_date) + + + # alt, az + # TODO: should they be in the ObservationTable? (they are derived quantities, like dead time) + # TODO: since I randomized RA/DEC without taking into account the observatory, I'm getting negative altitudes!!! + # maybe I should randomize alt az, then transform to RA/DEC!!! + observatory_location = observatory_locations[observatory] + + # ref: http://astropy.readthedocs.org/en/latest/coordinates/observing-example.html + sky_coord = SkyCoord(astro_table['RA'], astro_table['DEC'], frame='icrs') + alt_az_coord = sky_coord.transform_to(AltAz(obstime=astro_table['DATE'], location=observatory_location)) + #print(alt_az_coord) + #col_alt = alt_az_coord... + + + # TODO: operate row by row (using Observation class) instead of by columns? + # TODO: general methods for filling obs tables run by run; this function should call it. + + obs_table = ObservationTable(astro_table) + + return obs_table diff --git a/gammapy/datasets/tests/test_make.py b/gammapy/datasets/tests/test_make.py index 41f116d522..39af3d1386 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, generate_observation_table +from ...obs import ObservationTable def test_make_test_psf_fits_table(): @@ -15,3 +16,9 @@ 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_generate_observation_table(): + observatory='HESS' + n_obs = 10 + obs_table = generate_observation_table(observatory, n_obs) + assert len(obs_table) == n_obs diff --git a/gammapy/obs/observation.py b/gammapy/obs/observation.py index 23d102f462..c37dc42f14 100644 --- a/gammapy/obs/observation.py +++ b/gammapy/obs/observation.py @@ -11,7 +11,7 @@ class Observation(object): - """Observation (a.k.a. Run). + """Observation (a.k.a. run). TODO: not clear if this class is useful. @@ -43,6 +43,9 @@ class ObservationTable(Table): * ``OBS_ID`` * ``ONTIME`` * ``LIVETIME`` + * ``RA`` + * ``DEC`` + * ``DATE`` * ... """ @@ -56,6 +59,7 @@ def info(self): ss += 'Total live time: {}\n'.format(livetime) dtf = 100. * (1 - livetime / ontime) ss += 'Average dead time fraction: {:5.2f}%'.format(dtf) + #TODO: units are not shown!!! return ss def select_linspace_subset(self, num): From 8568fa71b3adce579899a1554a9bcb2636371df8 Mon Sep 17 00:00:00 2001 From: mapaz Date: Wed, 3 Jun 2015 16:49:43 -0700 Subject: [PATCH 02/15] Fixed astropy.Time format incompatibility with current astropy stable version. --- gammapy/datasets/make.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/gammapy/datasets/make.py b/gammapy/datasets/make.py index 220bebce6a..c9944250ac 100644 --- a/gammapy/datasets/make.py +++ b/gammapy/datasets/make.py @@ -131,13 +131,14 @@ def generate_observation_table(observatory, n_obs): # random points between the start of 2010 and the end of 2014 #TODO: should this represent the time at the beginning of the run? # this has consequences for the ra/dec -> alt/az conversion - datestart = Time('2010-01-01T00:00:00', format='fits', scale='utc') - dateend = Time('2015-01-01T00:00:00', format='fits', scale='utc') - date = Time((dateend.mjd - datestart.mjd)*np.random.random(len(col_obs_id)) + datestart.mjd, format='mjd', scale='utc').fits + datestart = Time('2010-01-01 00:00:00', format='iso', scale='utc') + dateend = Time('2015-01-01 00:00:00', format='iso', scale='utc') + date = Time((dateend.mjd - datestart.mjd)*np.random.random(len(col_obs_id)) + datestart.mjd, format='mjd', scale='utc').iso + # TODO: using format "iso" for ~astropy.Time for now; eventually change it + # to "fits" after the next astropy stable release (>1.0) is out. col_date = Column(name='DATE', data=date) astro_table.add_column(col_date) - # alt, az # TODO: should they be in the ObservationTable? (they are derived quantities, like dead time) # TODO: since I randomized RA/DEC without taking into account the observatory, I'm getting negative altitudes!!! From 90f7a08636fa9464da13f7eac533ab3a8570acb5 Mon Sep 17 00:00:00 2001 From: mapaz Date: Thu, 4 Jun 2015 13:45:04 -0700 Subject: [PATCH 03/15] In random obs table generator: calculating first ALT, AZ, then converting to RA, Dec. --- gammapy/datasets/make.py | 105 +++++++++++++++++++++++-------------- gammapy/obs/observation.py | 5 +- gammapy/utils/random.py | 18 ++++--- 3 files changed, 83 insertions(+), 45 deletions(-) diff --git a/gammapy/datasets/make.py b/gammapy/datasets/make.py index c9944250ac..a3c9ffe988 100644 --- a/gammapy/datasets/make.py +++ b/gammapy/datasets/make.py @@ -7,9 +7,10 @@ from astropy.units import Quantity from astropy.time import Time from astropy.table import Table, Column -from astropy.coordinates import SkyCoord, AltAz, Angle +from astropy.coordinates import SkyCoord, AltAz, FK5, Angle from ..irf import EnergyDependentMultiGaussPSF from ..obs import ObservationTable, observatory_locations +from ..utils.random import sample_sphere __all__ = ['make_test_psf', 'generate_observation_table', @@ -109,52 +110,80 @@ def generate_observation_table(observatory, n_obs): astro_table.add_column(col_livetime) # TODO: add columns for coordinates: - # pointing observation -> done - # and date -> done + # alt az + # date -> done + # calculate pointing observation (ra dec), considering the observatory location, and alt az is at the middle of the observation # then convert to alt az, considering the observatory location - # TODO: add column(s) for offset, and take it into account in the coord. transformation!!! - # RA, Dec - # random points on a sphere ref: http://mathworld.wolfram.com/SpherePointPicking.html + # TODO: is there a way to comment on the column names?!!!!! + # TODO: store obs name as a column or as a header value?!!! + # TODO: format times!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + # TODO: alt az at mean time!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + # TODO: restrict to night time? (or dark time?)!!!!!!!!!!!! - ra = Angle(360.*np.random.random(len(col_obs_id)), 'degree') - col_ra = Column(name='RA', data=ra) - astro_table.add_column(col_ra) - - dec = Angle(np.arccos(2.*np.random.random(len(col_obs_id)) - 1), 'radian').to('degree') - # translate angles from [0, 180) deg to [-90, 90) deg - dec = dec - Angle(90., 'degree') - col_dec = Column(name='DEC', data=dec) - astro_table.add_column(col_dec) - - # date + # start time # random points between the start of 2010 and the end of 2014 - #TODO: should this represent the time at the beginning of the run? - # this has consequences for the ra/dec -> alt/az conversion + # TODO: restrict to night time? (or dark time?)!!! + # if so: take into acount that enough time has to be permitted for the observation to finish (~ 30 min) datestart = Time('2010-01-01 00:00:00', format='iso', scale='utc') dateend = Time('2015-01-01 00:00:00', format='iso', scale='utc') - date = Time((dateend.mjd - datestart.mjd)*np.random.random(len(col_obs_id)) + datestart.mjd, format='mjd', scale='utc').iso + time_start = Time((dateend.mjd - datestart.mjd)*np.random.random(len(col_obs_id)) + datestart.mjd, format='mjd', scale='utc').iso # TODO: using format "iso" for ~astropy.Time for now; eventually change it # to "fits" after the next astropy stable release (>1.0) is out. - col_date = Column(name='DATE', data=date) - astro_table.add_column(col_date) - - # alt, az - # TODO: should they be in the ObservationTable? (they are derived quantities, like dead time) - # TODO: since I randomized RA/DEC without taking into account the observatory, I'm getting negative altitudes!!! - # maybe I should randomize alt az, then transform to RA/DEC!!! - observatory_location = observatory_locations[observatory] - - # ref: http://astropy.readthedocs.org/en/latest/coordinates/observing-example.html - sky_coord = SkyCoord(astro_table['RA'], astro_table['DEC'], frame='icrs') - alt_az_coord = sky_coord.transform_to(AltAz(obstime=astro_table['DATE'], location=observatory_location)) - #print(alt_az_coord) - #col_alt = alt_az_coord... - - - # TODO: operate row by row (using Observation class) instead of by columns? - # TODO: general methods for filling obs tables run by run; this function should call it. + col_time_start = Column(name='TSTART', data=time_start) + astro_table.add_column(col_time_start) + + # stop time + # calculated as TSTART + ONTIME + time_stop = Time(astro_table['TSTART']) + astro_table['ONTIME'] + col_time_stop = Column(name='TSTOP', data=time_stop) + astro_table.add_column(col_time_stop) + + # az, alt + # random points in a sphere above 45 deg altitude + az, alt = sample_sphere(len(col_obs_id), (0, 360), (45, 90), 'degree') + az = Angle(az, 'degree') + alt = Angle(alt, 'degree') + col_az = Column(name='AZ', data=az) + astro_table.add_column(col_az) + col_alt = Column(name='ALT', data=alt) + astro_table.add_column(col_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 (TSTART + TSTOP)/2 (or TSTART + (ONTIME/2)) + az = Angle(astro_table['AZ']) + alt = Angle(astro_table['ALT']) + obstime = astro_table['TSTART'] + ##obstime = astro_table['TSTART'] + astro_table['TSTOP'] + ##obstime = Time(astro_table['TSTART']) + Time(astro_table['TSTOP']) + location = observatory_locations[observatory] + alt_az_coord = AltAz(az = az, alt = alt, obstime = obstime, location = location) + # TODO: make it depend on other pars: temperature, pressure, humidity,... + sky_coord = alt_az_coord.transform_to(FK5) + ra = sky_coord.ra + col_ra = Column(name='RA', data=ra) + astro_table.add_column(col_ra) + dec = sky_coord.dec + col_dec = Column(name='DEC', data=dec) + astro_table.add_column(col_dec) obs_table = ObservationTable(astro_table) +# #t1 = Time('2010-01-01 00:00:00') +# #t2 = Time('2010-02-01 00:00:00') +# #dt = t2 - t1 +# #dt +# #print(dt) +# #print(dt.iso) +# +# +# #from astropy.time import TimeDelta +### from astropy.utils.data import download_file +### from astropy.utils import iers +### from astropy.coordinates import builtin_frames +### iers.IERS.iers_table = iers.IERS_A.open(download_file(iers.IERS_A_URL, cache=True)) [builtin_frames.utils] +# +# #dtss = astro_table['TSTART'] + astro_table['TSTOP'] +# #dtss = Time(astro_table['TSTART']) + Time(astro_table['TSTOP']) + return obs_table diff --git a/gammapy/obs/observation.py b/gammapy/obs/observation.py index c37dc42f14..c4f787ab71 100644 --- a/gammapy/obs/observation.py +++ b/gammapy/obs/observation.py @@ -43,9 +43,12 @@ class ObservationTable(Table): * ``OBS_ID`` * ``ONTIME`` * ``LIVETIME`` + * ``TSTART`` + * ``TSTOP`` + * ``AZ`` + * ``ALT`` * ``RA`` * ``DEC`` - * ``DATE`` * ... """ 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): From 7b1bb3db0f5088529cfd7529156b0bad453b4ddd Mon Sep 17 00:00:00 2001 From: mapaz Date: Thu, 4 Jun 2015 15:25:51 -0700 Subject: [PATCH 04/15] Moving time_ref_from_dict from data/utils to utils/time. --- docs/utils/index.rst | 4 +++- gammapy/data/event_list.py | 5 +++-- gammapy/data/gti.py | 5 +++-- gammapy/data/utils.py | 11 ----------- gammapy/utils/tests/test_time.py | 21 +++++++++++++++++++++ gammapy/utils/time.py | 30 ++++++++++++++++++++++++++++++ 6 files changed, 60 insertions(+), 16 deletions(-) create mode 100644 gammapy/utils/tests/test_time.py create mode 100644 gammapy/utils/time.py 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/utils/tests/test_time.py b/gammapy/utils/tests/test_time.py new file mode 100644 index 0000000000..83b3a9c7b6 --- /dev/null +++ b/gammapy/utils/tests/test_time.py @@ -0,0 +1,21 @@ +# 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 ...utils import time_ref_from_dict +#from ...utils.time import time_ref_from_dict +#from ..utils.time import time_ref_from_dict +#from .utils.time import time_ref_from_dict +#from ..time import time_ref_from_dict +from gammapy.utils.time import time_ref_from_dict #OK! +#from . import time_ref_from_dict +#from .. import time_ref_from_dict + +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) + decimal = 4 + s_error = "time reference not compatible with defined values" + np.testing.assert_almost_equal(time_ref.mjd, mjd_int + mjd_frac, decimal, s_error) diff --git a/gammapy/utils/time.py b/gammapy/utils/time.py new file mode 100644 index 0000000000..e4837d2396 --- /dev/null +++ b/gammapy/utils/time.py @@ -0,0 +1,30 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +"""Time related utility functions.""" +from astropy.time import Time + +__all__ = ['time_ref_from_dict', + ] + + +def time_ref_from_dict(meta): + """ Calculte 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 + dictionnary with the keywords 'MJDREFI' and 'MJDREFF' + + Returns + ------- + `~astropy.Time` object with the reference time in units of 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) From 1388d8d9809d47928faf82589beb522204e9b364 Mon Sep 17 00:00:00 2001 From: mapaz Date: Fri, 5 Jun 2015 16:56:01 -0700 Subject: [PATCH 05/15] Renamed generate_observation_table to make_test_observation_table. --- gammapy/datasets/make.py | 48 +++++++++++++++-------------- gammapy/datasets/tests/test_make.py | 6 ++-- gammapy/utils/tests/test_time.py | 9 +----- 3 files changed, 29 insertions(+), 34 deletions(-) diff --git a/gammapy/datasets/make.py b/gammapy/datasets/make.py index a3c9ffe988..482c97dbac 100644 --- a/gammapy/datasets/make.py +++ b/gammapy/datasets/make.py @@ -13,7 +13,7 @@ from ..utils.random import sample_sphere __all__ = ['make_test_psf', - 'generate_observation_table', + 'make_test_observation_table', ] @@ -73,7 +73,7 @@ def sigma_energy_theta(energy, theta, sigma): return psf -def generate_observation_table(observatory, n_obs): +def make_test_observation_table(observatory, n_obs): """Generate an observation table. For the moment, only random observation tables are created. @@ -92,22 +92,22 @@ def generate_observation_table(observatory, n_obs): """ n_obs_start = 1 - astro_table = Table() + obs_table = ObservationTable() # obs id obs_id = np.arange(n_obs_start, n_obs_start + n_obs) col_obs_id = Column(name='OBS_ID', data=obs_id) - astro_table.add_column(col_obs_id) + obs_table.add_column(col_obs_id) # on time ontime = Quantity(30.*np.ones_like(col_obs_id.data), 'minute') col_ontime = Column(name='ONTIME', data=ontime) - astro_table.add_column(col_ontime) + obs_table.add_column(col_ontime) # livetime livetime = Quantity(25.*np.ones_like(col_obs_id.data), 'minute') col_livetime = Column(name='LIVETIME', data=livetime) - astro_table.add_column(col_livetime) + obs_table.add_column(col_livetime) # TODO: add columns for coordinates: # alt az @@ -120,24 +120,28 @@ def generate_observation_table(observatory, n_obs): # TODO: format times!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! # TODO: alt az at mean time!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! # TODO: restrict to night time? (or dark time?)!!!!!!!!!!!! + # TODO: in utils.time + # - a function to convert a certain time to time ref (int + fraction) + seconds after ref + # - a function to convert a certain time to seconds after a given time ref (int + fraction) # 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 # TODO: restrict to night time? (or dark time?)!!! - # if so: take into acount that enough time has to be permitted for the observation to finish (~ 30 min) + # if so: take into account that enough time has to be granted for the last observation of the night to finish (~ 30 min) 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(col_obs_id)) + datestart.mjd, format='mjd', scale='utc').iso # TODO: using format "iso" for ~astropy.Time for now; eventually change it # to "fits" after the next astropy stable release (>1.0) is out. col_time_start = Column(name='TSTART', data=time_start) - astro_table.add_column(col_time_start) + obs_table.add_column(col_time_start) # stop time # calculated as TSTART + ONTIME - time_stop = Time(astro_table['TSTART']) + astro_table['ONTIME'] + time_stop = Time(obs_table['TSTART']) + obs_table['ONTIME'] col_time_stop = Column(name='TSTOP', data=time_stop) - astro_table.add_column(col_time_stop) + obs_table.add_column(col_time_stop) # az, alt # random points in a sphere above 45 deg altitude @@ -145,29 +149,27 @@ def generate_observation_table(observatory, n_obs): az = Angle(az, 'degree') alt = Angle(alt, 'degree') col_az = Column(name='AZ', data=az) - astro_table.add_column(col_az) + obs_table.add_column(col_az) col_alt = Column(name='ALT', data=alt) - astro_table.add_column(col_alt) + obs_table.add_column(col_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 (TSTART + TSTOP)/2 (or TSTART + (ONTIME/2)) - az = Angle(astro_table['AZ']) - alt = Angle(astro_table['ALT']) - obstime = astro_table['TSTART'] - ##obstime = astro_table['TSTART'] + astro_table['TSTOP'] - ##obstime = Time(astro_table['TSTART']) + Time(astro_table['TSTOP']) + az = Angle(obs_table['AZ']) + alt = Angle(obs_table['ALT']) + obstime = obs_table['TSTART'] + ##obstime = obs_table['TSTART'] + obs_table['TSTOP'] + ##obstime = Time(obs_table['TSTART']) + Time(obs_table['TSTOP']) location = observatory_locations[observatory] alt_az_coord = AltAz(az = az, alt = alt, obstime = obstime, location = location) # TODO: make it depend on other pars: temperature, pressure, humidity,... sky_coord = alt_az_coord.transform_to(FK5) ra = sky_coord.ra col_ra = Column(name='RA', data=ra) - astro_table.add_column(col_ra) + obs_table.add_column(col_ra) dec = sky_coord.dec col_dec = Column(name='DEC', data=dec) - astro_table.add_column(col_dec) - - obs_table = ObservationTable(astro_table) + obs_table.add_column(col_dec) # #t1 = Time('2010-01-01 00:00:00') # #t2 = Time('2010-02-01 00:00:00') @@ -183,7 +185,7 @@ def generate_observation_table(observatory, n_obs): ### from astropy.coordinates import builtin_frames ### iers.IERS.iers_table = iers.IERS_A.open(download_file(iers.IERS_A_URL, cache=True)) [builtin_frames.utils] # -# #dtss = astro_table['TSTART'] + astro_table['TSTOP'] -# #dtss = Time(astro_table['TSTART']) + Time(astro_table['TSTOP']) +# #dtss = obs_table['TSTART'] + obs_table['TSTOP'] +# #dtss = Time(obs_table['TSTART']) + Time(obs_table['TSTOP']) return obs_table diff --git a/gammapy/datasets/tests/test_make.py b/gammapy/datasets/tests/test_make.py index 39af3d1386..004273c83e 100644 --- a/gammapy/datasets/tests/test_make.py +++ b/gammapy/datasets/tests/test_make.py @@ -4,7 +4,7 @@ from astropy.coordinates import Angle from astropy.units import Quantity from ...utils.testing import assert_quantity -from ...datasets import make_test_psf, generate_observation_table +from ...datasets import make_test_psf, make_test_observation_table from ...obs import ObservationTable @@ -17,8 +17,8 @@ def test_make_test_psf_fits_table(): fraction = psf2.containment_fraction(0.1) assert fraction == 0.44485638998490795 -def test_generate_observation_table(): +def test_make_test_observation_table(): observatory='HESS' n_obs = 10 - obs_table = generate_observation_table(observatory, n_obs) + obs_table = make_test_observation_table(observatory, n_obs) assert len(obs_table) == n_obs diff --git a/gammapy/utils/tests/test_time.py b/gammapy/utils/tests/test_time.py index 83b3a9c7b6..e1d8226669 100644 --- a/gammapy/utils/tests/test_time.py +++ b/gammapy/utils/tests/test_time.py @@ -2,14 +2,7 @@ from __future__ import (absolute_import, division, print_function, unicode_literals) import numpy as np -#from ...utils import time_ref_from_dict -#from ...utils.time import time_ref_from_dict -#from ..utils.time import time_ref_from_dict -#from .utils.time import time_ref_from_dict -#from ..time import time_ref_from_dict -from gammapy.utils.time import time_ref_from_dict #OK! -#from . import time_ref_from_dict -#from .. import time_ref_from_dict +from ..time import time_ref_from_dict def test_time_ref_from_dict(): mjd_int = 500 From dba3aeaa30ebc1658f0fc927ade345899d408388 Mon Sep 17 00:00:00 2001 From: mapaz Date: Fri, 5 Jun 2015 16:56:44 -0700 Subject: [PATCH 06/15] Docs: defining obs lists format for Gammapy. --- docs/dataformats/observation_lists.rst | 47 ++++++++++++++++++++------ 1 file changed, 36 insertions(+), 11 deletions(-) diff --git a/docs/dataformats/observation_lists.rst b/docs/dataformats/observation_lists.rst index a5cb0b3239..ac4a4ca8fc 100644 --- a/docs/dataformats/observation_lists.rst +++ b/docs/dataformats/observation_lists.rst @@ -22,23 +22,48 @@ This has the advantage that everyone can work with whatever tool they like. Here * `csvkit `_ * Your own script using e.g the `Python standard library CSV module `_ or `pandas `_ -A run list must have at least a column called ``Run``:: +A run list must have at least a column called ``OBS_ID``:: - Run + OBS_ID 42 43 -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 +Usually it has many more columns with information about each observation. A list of Gammapy supported columns is:: + + OBS_ID: observation ID as an integer (starting at 0 or 1? does it matter?) + RA: pointing position right ascension in equatorial (ICRS) coordinates + DEC: pointing position declination in equatorial (ICRS) coordinates + AZ: average azimuth angle during the observarion + ALT: average altitude angle during the observarion + MUON_EFFICIENCY: average muon efficiency of the telescopes + TIME_START: start time of the observation stored as number of seconds after the time reference in the header + TIME_STOP: end time of the observation in the same forma as TIME_START + TIME_OBSERVATION: duration of the observation + TIME_LIVE: duration of the observation without dead time + TRIGGER_RATE: average trigger rate of the system + MEAN_TEMPERATURE: average temperature of the atmosphere + N_TELS: number of telescopes participating in the observation + TEL_LIST: string with a CSV list of telescope IDs participating in the observation + QUALITY: data quality; recommended: "spectral" or "detection" (not used by Gammapy at the moment) + +Extra user defined columns are allowed; Gammapy will ignore them. + +In order for the extra columns to have full meaning a header is needed with at least the following keywords:: + + OBSERVATORY_NAME: name of the observatory where the observations were taken. This is important for instance for coordinate transformations between celestial (i.e. RA/dec) and terrestrial (i.e. az/alt) coordinate systems. + TIME_REF_MJD_INT: reference time for other times in the list (i.e. TIME_START/TIME_STOP). Integer value in mean julian days. + TIME_REF_MJD_FRA: fraction of integer value defined in TIME_REF_MJD_INT. + +TODO: change names already defined in event lists (and in `~gammapy.utils.time.time_ref_from_dict`) ``MJDREFI`` and ``MJDREFF``? Or adopt ``TIME_REF_MJD_INT`` and ``TIME_REF_MJD_FRA``? + +TODO: should the observation list have already a header? Or should the header values be defined as columns, then interpreted correspondingly when stored into an ObservationTable? (I might be mixing up 2 things: obs lists and obs tables here? or should they have the same format?) + +TODO: should the observation list already define the times in this ``TIME_REF_MJD_INT + TIME_REF_MJD_FRA`` format? Or should it be converted once its read into an ObservationTable? + +TODO: ``TEL_LIST``: incompatibility wit CSV format, since it's a CSV list of tels???!!! -Special column names that the Gammapy analysis tools understand: +TODO: how to biuld sphinx only for a specific module or file? -* ``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 XML format ---------- From 9d2f1d1f4c5a73f0254265f79ba0b95eaac9337b Mon Sep 17 00:00:00 2001 From: mapaz Date: Mon, 8 Jun 2015 12:40:09 -0700 Subject: [PATCH 07/15] Added utility function to calculate time differences w.r.t. a time reference. --- gammapy/utils/tests/test_time.py | 17 ++++++++++++++++- gammapy/utils/time.py | 27 +++++++++++++++++++++++++-- 2 files changed, 41 insertions(+), 3 deletions(-) diff --git a/gammapy/utils/tests/test_time.py b/gammapy/utils/tests/test_time.py index e1d8226669..bac1ab0e6d 100644 --- a/gammapy/utils/tests/test_time.py +++ b/gammapy/utils/tests/test_time.py @@ -2,7 +2,9 @@ from __future__ import (absolute_import, division, print_function, unicode_literals) import numpy as np -from ..time import time_ref_from_dict +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 @@ -12,3 +14,16 @@ def test_time_ref_from_dict(): decimal = 4 s_error = "time reference not compatible with defined values" np.testing.assert_almost_equal(time_ref.mjd, mjd_int + mjd_frac, decimal, s_error) + + +def test_time_relative_to_ref(): + 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) + delta_time_1sec = TimeDelta(1., format='sec') + time = time_ref + delta_time_1sec + delta_time = time_relative_to_ref(time, time_ref_dict) + decimal = 4 + s_error = "delta time not compatible with defined value" + np.testing.assert_almost_equal(delta_time.sec, delta_time_1sec.sec, decimal, s_error) diff --git a/gammapy/utils/time.py b/gammapy/utils/time.py index e4837d2396..6313b711eb 100644 --- a/gammapy/utils/time.py +++ b/gammapy/utils/time.py @@ -1,13 +1,14 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst """Time related utility functions.""" -from astropy.time import Time +from astropy.time import Time, TimeDelta __all__ = ['time_ref_from_dict', + 'time_relative_to_ref', ] def time_ref_from_dict(meta): - """ Calculte the time reference from metadata. + """Calculte 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. @@ -28,3 +29,25 @@ def time_ref_from_dict(meta): # 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 to be converted. + meta : dict + dictionnary with the keywords 'MJDREFI' and 'MJDREFF' + + Returns + ------- + `~astropy.TimeDelta` object with the time in seconds after the reference. + """ + time_ref = time_ref_from_dict(meta) + delta_time = TimeDelta(time - time_ref, format='sec') + + return delta_time From 8b49c31ef10cb4686b2b37be697b0f980bac69af Mon Sep 17 00:00:00 2001 From: mapaz Date: Tue, 9 Jun 2015 12:36:05 -0700 Subject: [PATCH 08/15] Generation of observation tables uses Quantity objects for time interval columns, because TimeDelta seems faulty. --- gammapy/datasets/make.py | 123 +++++++++++++++++----------- gammapy/datasets/tests/test_make.py | 5 +- 2 files changed, 77 insertions(+), 51 deletions(-) diff --git a/gammapy/datasets/make.py b/gammapy/datasets/make.py index 482c97dbac..9668505959 100644 --- a/gammapy/datasets/make.py +++ b/gammapy/datasets/make.py @@ -5,12 +5,13 @@ unicode_literals) import numpy as np from astropy.units import Quantity -from astropy.time import Time +from astropy.time import Time, TimeDelta from astropy.table import Table, Column 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', 'make_test_observation_table', @@ -73,17 +74,19 @@ def sigma_energy_theta(energy, theta, sigma): return psf -def make_test_observation_table(observatory, n_obs): +def make_test_observation_table(observatory_name, n_obs, debug=False): """Generate an observation table. For the moment, only random observation tables are created. Parameters ---------- - observatory : string + observatory_name : string name of the observatory n_obs : integer number of observations for the obs table + debug : bool + show UTC times instead of seconds after the reference Returns ------- @@ -94,52 +97,84 @@ def make_test_observation_table(observatory, n_obs): obs_table = ObservationTable() + # build a time reference as the start of 2010 + dateref = Time('2010-01-01 00:00:00', format='iso', scale='utc') + # TODO: using format "iso" for `~astropy.Time` for now; eventually change it + # to "fits" after the next astropy stable release (>1.0) is out. + 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) col_obs_id = Column(name='OBS_ID', data=obs_id) obs_table.add_column(col_obs_id) - # on time - ontime = Quantity(30.*np.ones_like(col_obs_id.data), 'minute') + # TODO: `~astropy.table.Table` doesn't handle `~astropy.time.TimeDelta` columns + # properly see: + # https://github.com/astropy/astropy/issues/3832 + # until this issue is solved (and included into a stable release), quantity + # objects are used. + + # on time: 30 min + ontime = TimeDelta(30.*60.*np.ones_like(col_obs_id.data), format='sec') + ontime = Quantity(ontime.sec, 'second') # converting to quantity col_ontime = Column(name='ONTIME', data=ontime) obs_table.add_column(col_ontime) - # livetime - livetime = Quantity(25.*np.ones_like(col_obs_id.data), 'minute') + # livetime: 25 min + livetime = TimeDelta(25.*60.*np.ones_like(col_obs_id.data), format='sec') + livetime = Quantity(livetime.sec, 'second') # converting to quantity col_livetime = Column(name='LIVETIME', data=livetime) obs_table.add_column(col_livetime) - # TODO: add columns for coordinates: - # alt az - # date -> done - # calculate pointing observation (ra dec), considering the observatory location, and alt az is at the middle of the observation - # then convert to alt az, considering the observatory location - - # TODO: is there a way to comment on the column names?!!!!! - # TODO: store obs name as a column or as a header value?!!! - # TODO: format times!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - # TODO: alt az at mean time!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - # TODO: restrict to night time? (or dark time?)!!!!!!!!!!!! - # TODO: in utils.time - # - a function to convert a certain time to time ref (int + fraction) + seconds after ref - # - a function to convert a certain time to seconds after a given time ref (int + fraction) + # TODO: adopt new name scheme defined in sphinx doc!!!!!! + # TODO: is there a way to comment on the column names?!!! # 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 - # TODO: restrict to night time? (or dark time?)!!! - # if so: take into account that enough time has to be granted for the last observation of the night to finish (~ 30 min) + # 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(col_obs_id)) + datestart.mjd, format='mjd', scale='utc').iso - # TODO: using format "iso" for ~astropy.Time for now; eventually change it - # to "fits" after the next astropy stable release (>1.0) is out. + time_start = Time((dateend.mjd - datestart.mjd)*np.random.random(len(col_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 = TimeDelta(22.*60.*60., format='sec') + night_duration = TimeDelta(5.5*60.*60., format='sec') + hour_start = night_start + TimeDelta(night_duration.sec*np.random.random(len(col_obs_id)), format='sec') + + # 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) + time_start = Quantity(time_start.sec, 'second') # converting to quantity + col_time_start = Column(name='TSTART', data=time_start) obs_table.add_column(col_time_start) # stop time # calculated as TSTART + ONTIME - time_stop = Time(obs_table['TSTART']) + obs_table['ONTIME'] + if debug : + time_stop = Time(obs_table['TSTART']) + TimeDelta(obs_table['ONTIME']) + else : + time_stop = TimeDelta(obs_table['TSTART']) + TimeDelta(obs_table['ONTIME']) + time_stop = Quantity(time_stop.sec, 'second') # converting to quantity + col_time_stop = Column(name='TSTOP', data=time_stop) obs_table.add_column(col_time_stop) @@ -154,15 +189,20 @@ def make_test_observation_table(observatory, n_obs): obs_table.add_column(col_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 (TSTART + TSTOP)/2 (or TSTART + (ONTIME/2)) + # derive from az, alt taking into account that alt, az represent the values + # at the middle of the observation, i.e. at time_ref + (TSTART + TSTOP)/2 + # (or better: time_ref + TSTART + (ONTIME/2)) + # in debug modus, the time_ref should not be added, since it's already included + # in TSTART and TSTOP az = Angle(obs_table['AZ']) alt = Angle(obs_table['ALT']) - obstime = obs_table['TSTART'] - ##obstime = obs_table['TSTART'] + obs_table['TSTOP'] - ##obstime = Time(obs_table['TSTART']) + Time(obs_table['TSTOP']) - location = observatory_locations[observatory] + if debug : + obstime = Time(obs_table['TSTART']) + TimeDelta(obs_table['ONTIME'])/2. + else : + obstime = time_ref_from_dict(obs_table.meta) + TimeDelta(obs_table['TSTART']) + TimeDelta(obs_table['ONTIME'])/2. + location = observatory_locations[observatory_name] alt_az_coord = AltAz(az = az, alt = alt, obstime = obstime, location = location) - # TODO: make it depend on other pars: temperature, pressure, humidity,... + # optional: make it depend on other pars: temperature, pressure, humidity,... sky_coord = alt_az_coord.transform_to(FK5) ra = sky_coord.ra col_ra = Column(name='RA', data=ra) @@ -171,21 +211,6 @@ def make_test_observation_table(observatory, n_obs): col_dec = Column(name='DEC', data=dec) obs_table.add_column(col_dec) -# #t1 = Time('2010-01-01 00:00:00') -# #t2 = Time('2010-02-01 00:00:00') -# #dt = t2 - t1 -# #dt -# #print(dt) -# #print(dt.iso) -# -# -# #from astropy.time import TimeDelta -### from astropy.utils.data import download_file -### from astropy.utils import iers -### from astropy.coordinates import builtin_frames -### iers.IERS.iers_table = iers.IERS_A.open(download_file(iers.IERS_A_URL, cache=True)) [builtin_frames.utils] -# -# #dtss = obs_table['TSTART'] + obs_table['TSTOP'] -# #dtss = Time(obs_table['TSTART']) + Time(obs_table['TSTOP']) + # optional: it would be nice to plot a skymap with the simulated RA/dec positions return obs_table diff --git a/gammapy/datasets/tests/test_make.py b/gammapy/datasets/tests/test_make.py index 004273c83e..ec769d6bac 100644 --- a/gammapy/datasets/tests/test_make.py +++ b/gammapy/datasets/tests/test_make.py @@ -17,8 +17,9 @@ def test_make_test_psf_fits_table(): fraction = psf2.containment_fraction(0.1) assert fraction == 0.44485638998490795 + def test_make_test_observation_table(): - observatory='HESS' + observatory_name='HESS' n_obs = 10 - obs_table = make_test_observation_table(observatory, n_obs) + obs_table = make_test_observation_table(observatory_name, n_obs) assert len(obs_table) == n_obs From a8bf41d222cec31bae9ecad8b417b149c63926ab Mon Sep 17 00:00:00 2001 From: mapaz Date: Tue, 9 Jun 2015 12:56:00 -0700 Subject: [PATCH 09/15] Adopting new naming scheme from the doc for the obs table columns. --- gammapy/datasets/make.py | 41 ++++++++++++++++++-------------------- gammapy/obs/observation.py | 12 +++++------ 2 files changed, 25 insertions(+), 28 deletions(-) diff --git a/gammapy/datasets/make.py b/gammapy/datasets/make.py index 9668505959..d1ff64120a 100644 --- a/gammapy/datasets/make.py +++ b/gammapy/datasets/make.py @@ -119,20 +119,17 @@ def make_test_observation_table(observatory_name, n_obs, debug=False): # until this issue is solved (and included into a stable release), quantity # objects are used. - # on time: 30 min - ontime = TimeDelta(30.*60.*np.ones_like(col_obs_id.data), format='sec') - ontime = Quantity(ontime.sec, 'second') # converting to quantity - col_ontime = Column(name='ONTIME', data=ontime) - obs_table.add_column(col_ontime) + # obs time: 30 min + time_observation = TimeDelta(30.*60.*np.ones_like(col_obs_id.data), format='sec') + time_observation = Quantity(time_observation.sec, 'second') # converting to quantity + col_time_observation = Column(name='TIME_OBSERVATION', data=time_observation) + obs_table.add_column(col_time_observation) # livetime: 25 min - livetime = TimeDelta(25.*60.*np.ones_like(col_obs_id.data), format='sec') - livetime = Quantity(livetime.sec, 'second') # converting to quantity - col_livetime = Column(name='LIVETIME', data=livetime) - obs_table.add_column(col_livetime) - - # TODO: adopt new name scheme defined in sphinx doc!!!!!! - # TODO: is there a way to comment on the column names?!!! + time_live = TimeDelta(25.*60.*np.ones_like(col_obs_id.data), format='sec') + time_live = Quantity(time_live.sec, 'second') # converting to quantity + col_time_live = Column(name='TIME_LIVE', data=time_live) + obs_table.add_column(col_time_live) # start time # random points between the start of 2010 and the end of 2014 @@ -164,18 +161,18 @@ def make_test_observation_table(observatory_name, n_obs, debug=False): time_start = time_relative_to_ref(time_start, header) time_start = Quantity(time_start.sec, 'second') # converting to quantity - col_time_start = Column(name='TSTART', data=time_start) + col_time_start = Column(name='TIME_START', data=time_start) obs_table.add_column(col_time_start) # stop time - # calculated as TSTART + ONTIME + # calculated as TIME_START + TIME_OBSERVATION if debug : - time_stop = Time(obs_table['TSTART']) + TimeDelta(obs_table['ONTIME']) + time_stop = Time(obs_table['TIME_START']) + TimeDelta(obs_table['TIME_OBSERVATION']) else : - time_stop = TimeDelta(obs_table['TSTART']) + TimeDelta(obs_table['ONTIME']) + time_stop = TimeDelta(obs_table['TIME_START']) + TimeDelta(obs_table['TIME_OBSERVATION']) time_stop = Quantity(time_stop.sec, 'second') # converting to quantity - col_time_stop = Column(name='TSTOP', data=time_stop) + col_time_stop = Column(name='TIME_STOP', data=time_stop) obs_table.add_column(col_time_stop) # az, alt @@ -190,16 +187,16 @@ def make_test_observation_table(observatory_name, n_obs, debug=False): # 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 + (TSTART + TSTOP)/2 - # (or better: time_ref + TSTART + (ONTIME/2)) + # 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 TSTART and TSTOP + # in TIME_START and TIME_STOP az = Angle(obs_table['AZ']) alt = Angle(obs_table['ALT']) if debug : - obstime = Time(obs_table['TSTART']) + TimeDelta(obs_table['ONTIME'])/2. + obstime = Time(obs_table['TIME_START']) + TimeDelta(obs_table['TIME_OBSERVATION'])/2. else : - obstime = time_ref_from_dict(obs_table.meta) + TimeDelta(obs_table['TSTART']) + TimeDelta(obs_table['ONTIME'])/2. + 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) # optional: make it depend on other pars: temperature, pressure, humidity,... diff --git a/gammapy/obs/observation.py b/gammapy/obs/observation.py index c4f787ab71..0f4b913311 100644 --- a/gammapy/obs/observation.py +++ b/gammapy/obs/observation.py @@ -41,10 +41,10 @@ class ObservationTable(Table): convenience methods and the following columns: * ``OBS_ID`` - * ``ONTIME`` - * ``LIVETIME`` - * ``TSTART`` - * ``TSTOP`` + * ``TIME_OBSERVATION`` + * ``TIME_LIVE`` + * ``TIME_START`` + * ``TIME_STOP`` * ``AZ`` * ``ALT`` * ``RA`` @@ -56,9 +56,9 @@ class ObservationTable(Table): def info(self): ss = 'Observation table:\n' ss += 'Number of observations: {}\n'.format(len(self)) - ontime = self['ONTIME'].sum() + ontime = self['TIME_OBSERVATION'].sum() ss += 'Total observation time: {}\n'.format(ontime) - livetime = self['LIVETIME'].sum() + livetime = self['TIME_LIVE'].sum() ss += 'Total live time: {}\n'.format(livetime) dtf = 100. * (1 - livetime / ontime) ss += 'Average dead time fraction: {:5.2f}%'.format(dtf) From c27dc7c827c4ffa050f2e8eb19e4a536d8f863f9 Mon Sep 17 00:00:00 2001 From: mapaz Date: Tue, 9 Jun 2015 15:01:21 -0700 Subject: [PATCH 10/15] Added columns for n_tels and muon_efficiency to table generator. --- gammapy/datasets/make.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/gammapy/datasets/make.py b/gammapy/datasets/make.py index d1ff64120a..e80fb7bbea 100644 --- a/gammapy/datasets/make.py +++ b/gammapy/datasets/make.py @@ -210,4 +210,20 @@ def make_test_observation_table(observatory_name, n_obs, debug=False): # optional: it would be nice to plot a skymap with the simulated RA/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(col_obs_id)) + col_n_tels = Column(name='N_TELS', data=n_tels) + obs_table.add_column(col_n_tels) + + # muon efficiency + # random between 0.6 and 1.0 + muon_efficiency_min = 0.6 + muon_efficiency_max = 1.0 + muon_efficiency = (muon_efficiency_max - muon_efficiency_min)*np.random.random(len(col_obs_id)) + muon_efficiency_min + col_muon_efficiency = Column(name='MUON_EFFICIENCY', data=muon_efficiency) + obs_table.add_column(col_muon_efficiency) + return obs_table From 76c67efc053780a0e3412034dcbf8f809f76a615 Mon Sep 17 00:00:00 2001 From: mapaz Date: Tue, 9 Jun 2015 15:07:30 -0700 Subject: [PATCH 11/15] Revised structure of observation lists doc and moved the file format stuff to its rightful place. --- docs/dataformats/file_formats.rst | 3 + docs/dataformats/observation_lists.rst | 76 ++++++++++++-------------- gammapy/obs/observation.py | 24 ++++---- 3 files changed, 49 insertions(+), 54 deletions(-) diff --git a/docs/dataformats/file_formats.rst b/docs/dataformats/file_formats.rst index 58ff09e7b4..1f253f132f 100644 --- a/docs/dataformats/file_formats.rst +++ b/docs/dataformats/file_formats.rst @@ -194,6 +194,9 @@ XML **XML** files (see `Wikipedia `__) +GammaLib / ctools uses an "observation definition" XML format described +`here `__. + TODO: describe FITS diff --git a/docs/dataformats/observation_lists.rst b/docs/dataformats/observation_lists.rst index ac4a4ca8fc..e3be601512 100644 --- a/docs/dataformats/observation_lists.rst +++ b/docs/dataformats/observation_lists.rst @@ -10,63 +10,59 @@ 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``. - -CSV format ----------- - -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: - -* `LibreOffice Calc `_ -* `TOPCAT `_ or `STILTS `_ -* `csvkit `_ -* Your own script using e.g the `Python standard library CSV module `_ or `pandas `_ - A run list must have at least a column called ``OBS_ID``:: OBS_ID 42 43 -Usually it has many more columns with information about each observation. A list of Gammapy supported columns is:: - - OBS_ID: observation ID as an integer (starting at 0 or 1? does it matter?) - RA: pointing position right ascension in equatorial (ICRS) coordinates - DEC: pointing position declination in equatorial (ICRS) coordinates - AZ: average azimuth angle during the observarion - ALT: average altitude angle during the observarion - MUON_EFFICIENCY: average muon efficiency of the telescopes - TIME_START: start time of the observation stored as number of seconds after the time reference in the header - TIME_STOP: end time of the observation in the same forma as TIME_START - TIME_OBSERVATION: duration of the observation - TIME_LIVE: duration of the observation without dead time - TRIGGER_RATE: average trigger rate of the system - MEAN_TEMPERATURE: average temperature of the atmosphere - N_TELS: number of telescopes participating in the observation - TEL_LIST: string with a CSV list of telescope IDs participating in the observation - QUALITY: data quality; recommended: "spectral" or "detection" (not used by Gammapy at the moment) +Usually it has many more columns with information about each observation. A list of Gammapy supported columns is: + +================ ================================================================================================ ========= +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 yes? +ALT average altitude angle during the observation yes? +MUON_EFFICIENCY average muon efficiency of the telescopes yes? +TIME_START start time of the observation stored as number of seconds after the time reference in the header yes? +TIME_STOP end time of the observation in the same forma as TIME_START yes? +TIME_OBSERVATION duration of the observation yes? +TIME_LIVE duration of the observation without dead time yes? +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" (not used by Gammapy at the moment) no +================ ================================================================================================ ========= Extra user defined columns are allowed; Gammapy will ignore them. -In order for the extra columns to have full meaning a header is needed with at least the following keywords:: +In order for the extra columns to have full meaning the following is needed: + + * 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: - OBSERVATORY_NAME: name of the observatory where the observations were taken. This is important for instance for coordinate transformations between celestial (i.e. RA/dec) and terrestrial (i.e. az/alt) coordinate systems. - TIME_REF_MJD_INT: reference time for other times in the list (i.e. TIME_START/TIME_STOP). Integer value in mean julian days. - TIME_REF_MJD_FRA: fraction of integer value defined in TIME_REF_MJD_INT. +================ ========================================================================================================================================================================================================== +keyword description +================ ========================================================================================================================================================================================================== +OBSERVATORY_NAME name of the observatory where the observations were taken. This is important for instance for coordinate transformations between celestial (i.e. RA/dec) and terrestrial (i.e. az/alt) coordinate systems. +MJDREFI reference time for other times in the list (i.e. TIME_START/TIME_STOP). Integer value in mean julian days. +MJDREFF fraction of integer value defined in MJDREFI. +================ ========================================================================================================================================================================================================== -TODO: change names already defined in event lists (and in `~gammapy.utils.time.time_ref_from_dict`) ``MJDREFI`` and ``MJDREFF``? Or adopt ``TIME_REF_MJD_INT`` and ``TIME_REF_MJD_FRA``? +Extra user defined header entries are allowed; Gammapy will ignore them. TODO: should the observation list have already a header? Or should the header values be defined as columns, then interpreted correspondingly when stored into an ObservationTable? (I might be mixing up 2 things: obs lists and obs tables here? or should they have the same format?) -TODO: should the observation list already define the times in this ``TIME_REF_MJD_INT + TIME_REF_MJD_FRA`` format? Or should it be converted once its read into an ObservationTable? +TODO: should the observation list already define the times in this ``MJDREFI + MJDREFF`` format? Or should it be converted once its read into an ObservationTable? TODO: ``TEL_LIST``: incompatibility wit CSV format, since it's a CSV list of tels???!!! TODO: how to biuld sphinx only for a specific module or file? +TODO: add an example table?!!! (and header? how?) -XML format ----------- - -GammaLib / ctools uses an "observation definition" XML format described -`here `__. +TODO: restructure (as requested in PR)!!! diff --git a/gammapy/obs/observation.py b/gammapy/obs/observation.py index 0f4b913311..f77613c875 100644 --- a/gammapy/obs/observation.py +++ b/gammapy/obs/observation.py @@ -3,6 +3,7 @@ unicode_literals) import numpy as np from astropy.table import Table +from ..utils.time import time_ref_from_dict __all__ = [ # 'Observation', @@ -38,30 +39,26 @@ 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`` - * ``TIME_OBSERVATION`` - * ``TIME_LIVE`` - * ``TIME_START`` - * ``TIME_STOP`` - * ``AZ`` - * ``ALT`` - * ``RA`` - * ``DEC`` - * ... + convenience methods. The format of the observation table + is described in: + http://gammapy.readthedocs.org/en/latest/dataformats/observation_lists.html + TODO: is there a better way to refer to the gammapy doc?!!! """ 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['TIME_OBSERVATION'].sum() ss += 'Total observation time: {}\n'.format(ontime) livetime = self['TIME_LIVE'].sum() 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) + ss += 'Time reference: {}'.format(time_ref) #TODO: units are not shown!!! return ss @@ -89,4 +86,3 @@ def select_linspace_subset(self, num): # Round down to nearest integer indices = indices.astype('int') return self[indices] - From 86008f55d55a302b5339e525233f0f9517b0d5c6 Mon Sep 17 00:00:00 2001 From: mapaz Date: Tue, 9 Jun 2015 15:50:10 -0700 Subject: [PATCH 12/15] Clean up doc and added an example table. --- docs/dataformats/observation_lists.rst | 45 +++++++++++++++++++++----- gammapy/datasets/make.py | 4 +-- gammapy/obs/observation.py | 19 ++++++----- gammapy/utils/tests/test_time.py | 8 ++--- gammapy/utils/time.py | 14 ++++---- 5 files changed, 61 insertions(+), 29 deletions(-) diff --git a/docs/dataformats/observation_lists.rst b/docs/dataformats/observation_lists.rst index e3be601512..1fb4128a3e 100644 --- a/docs/dataformats/observation_lists.rst +++ b/docs/dataformats/observation_lists.rst @@ -55,14 +55,43 @@ MJDREFF fraction of integer value defined in MJDREFI. Extra user defined header entries are allowed; Gammapy will ignore them. -TODO: should the observation list have already a header? Or should the header values be defined as columns, then interpreted correspondingly when stored into an ObservationTable? (I might be mixing up 2 things: obs lists and obs tables here? or should they have the same format?) - -TODO: should the observation list already define the times in this ``MJDREFI + MJDREFF`` format? Or should it be converted once its read into an ObservationTable? - -TODO: ``TEL_LIST``: incompatibility wit CSV format, since it's a CSV list of tels???!!! - -TODO: how to biuld sphinx only for a specific module or file? +TODO: ``TEL_LIST``: since the final format is not clear, please mind having a look to this issue and its outcome: +https://github.com/gammapy/gammapy/issues/282 TODO: add an example table?!!! (and header? how?) -TODO: restructure (as requested in PR)!!! +Example +------- +The tool `~gammapy.datasets.make.make_test_observation_table` can generate a `~gammapy.obs.ObservationTable` with dummy values. + +Header (metadata):: + + {u'MJDREFI': 55197.0, u'OBSERVATORY_NAME': u'HESS', u'MJDREFF': 0.0} + +Table: + ++------+----------------+---------+-------------+-------------+-------------+-------------+-------------+--------------+------+---------------+ +|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/gammapy/datasets/make.py b/gammapy/datasets/make.py index e80fb7bbea..95fc447edc 100644 --- a/gammapy/datasets/make.py +++ b/gammapy/datasets/make.py @@ -75,14 +75,14 @@ def sigma_energy_theta(energy, theta, sigma): def make_test_observation_table(observatory_name, n_obs, debug=False): - """Generate an observation table. + """Make a test observation table. For the moment, only random observation tables are created. Parameters ---------- observatory_name : string - name of the observatory + 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 diff --git a/gammapy/obs/observation.py b/gammapy/obs/observation.py index f77613c875..6f74805952 100644 --- a/gammapy/obs/observation.py +++ b/gammapy/obs/observation.py @@ -3,6 +3,8 @@ 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__ = [ @@ -40,10 +42,7 @@ class ObservationTable(Table): This is an `~astropy.table.Table` sub-class, with a few convenience methods. The format of the observation table - is described in: - http://gammapy.readthedocs.org/en/latest/dataformats/observation_lists.html - TODO: is there a better way to refer to the gammapy doc?!!! - + is described in :ref:`dataformats_observation_lists`. """ def info(self): @@ -51,15 +50,19 @@ def info(self): obs_name = self.meta['OBSERVATORY_NAME'] ss += 'Observatory name: {}\n'.format(obs_name) ss += 'Number of observations: {}\n'.format(len(self)) - ontime = self['TIME_OBSERVATION'].sum() + ontime = Quantity(self['TIME_OBSERVATION'].sum(), self['TIME_OBSERVATION'].unit) ss += 'Total observation time: {}\n'.format(ontime) - livetime = self['TIME_LIVE'].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}%\n'.format(dtf) time_ref = time_ref_from_dict(self.meta) - ss += 'Time reference: {}'.format(time_ref) - #TODO: units are not shown!!! + time_ref_unit = time_ref_from_dict(self.meta).format + ss += 'Time reference: {} {}'.format(time_ref, time_ref_unit) + # TODO: for an unknown reason, I can't enable astropy.units.cds + # as explained in the following link + # http://astropy.readthedocs.org/en/latest/units/#module-astropy.units.cds + # for using a Quantity object with MJD units. return ss def select_linspace_subset(self, num): diff --git a/gammapy/utils/tests/test_time.py b/gammapy/utils/tests/test_time.py index bac1ab0e6d..985845d690 100644 --- a/gammapy/utils/tests/test_time.py +++ b/gammapy/utils/tests/test_time.py @@ -1,7 +1,7 @@ # 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 numpy.testing import assert_almost_equal from astropy.time import Time, TimeDelta from ..time import time_ref_from_dict, time_relative_to_ref @@ -12,8 +12,7 @@ def test_time_ref_from_dict(): time_ref_dict = dict(MJDREFI=mjd_int, MJDREFF=mjd_frac) time_ref = time_ref_from_dict(time_ref_dict) decimal = 4 - s_error = "time reference not compatible with defined values" - np.testing.assert_almost_equal(time_ref.mjd, mjd_int + mjd_frac, decimal, s_error) + assert_almost_equal(time_ref.mjd, mjd_int + mjd_frac, decimal) def test_time_relative_to_ref(): @@ -25,5 +24,4 @@ def test_time_relative_to_ref(): time = time_ref + delta_time_1sec delta_time = time_relative_to_ref(time, time_ref_dict) decimal = 4 - s_error = "delta time not compatible with defined value" - np.testing.assert_almost_equal(delta_time.sec, delta_time_1sec.sec, decimal, s_error) + assert_almost_equal(delta_time.sec, delta_time_1sec.sec, decimal) diff --git a/gammapy/utils/time.py b/gammapy/utils/time.py index 6313b711eb..d8ea0a5f04 100644 --- a/gammapy/utils/time.py +++ b/gammapy/utils/time.py @@ -1,5 +1,7 @@ # 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', @@ -8,18 +10,18 @@ def time_ref_from_dict(meta): - """Calculte the time reference from metadata. + """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 - dictionnary with the keywords 'MJDREFI' and 'MJDREFF' + dictionary with the keywords `MJDREFI` and `MJDREFF` Returns ------- - `~astropy.Time` object with the reference time in units of MJD. + `~astropy.time.Time` object with the reference time in format of `MJD`. """ # Note: the `float` call here is to make sure we use 64-bit mjd = float(meta['MJDREFI']) + float(meta['MJDREFF']) @@ -38,14 +40,14 @@ def time_relative_to_ref(time, meta): Parameters ---------- - time: `~astropy.Time` + time: `~astropy.time.Time` time to be converted. meta : dict - dictionnary with the keywords 'MJDREFI' and 'MJDREFF' + dictionary with the keywords `MJDREFI` and `MJDREFF` Returns ------- - `~astropy.TimeDelta` object with the time in seconds after the reference. + `~astropy.time.TimeDelta` object with the time in seconds after the reference. """ time_ref = time_ref_from_dict(meta) delta_time = TimeDelta(time - time_ref, format='sec') From 4e699e867693d395a3e5f6c35f26779f6dbac4e1 Mon Sep 17 00:00:00 2001 From: mapaz Date: Wed, 10 Jun 2015 16:25:03 -0700 Subject: [PATCH 13/15] Using Quantity instead of TimeDelta for many columns in the dummy obs table generator. --- docs/dataformats/observation_lists.rst | 5 +- gammapy/datasets/make.py | 83 +++++++++++--------------- gammapy/utils/time.py | 7 ++- 3 files changed, 43 insertions(+), 52 deletions(-) diff --git a/docs/dataformats/observation_lists.rst b/docs/dataformats/observation_lists.rst index 1fb4128a3e..e00238f7d7 100644 --- a/docs/dataformats/observation_lists.rst +++ b/docs/dataformats/observation_lists.rst @@ -19,7 +19,7 @@ A run list must have at least a column called ``OBS_ID``:: Usually it has many more columns with information about each observation. A list of Gammapy supported columns is: ================ ================================================================================================ ========= -column name description required? +column name description required? ================ ================================================================================================ ========= OBS_ID observation ID as an integer yes RA pointing position right ascension in equatorial (ICRS) coordinates yes? @@ -58,11 +58,10 @@ Extra user defined header entries are allowed; Gammapy will ignore them. TODO: ``TEL_LIST``: since the final format is not clear, please mind having a look to this issue and its outcome: https://github.com/gammapy/gammapy/issues/282 -TODO: add an example table?!!! (and header? how?) Example ------- -The tool `~gammapy.datasets.make.make_test_observation_table` can generate a `~gammapy.obs.ObservationTable` with dummy values. +The tool `~gammapy.datasets.make_test_observation_table` can generate a `~gammapy.obs.ObservationTable` with dummy values. Header (metadata):: diff --git a/gammapy/datasets/make.py b/gammapy/datasets/make.py index 95fc447edc..33948dceb0 100644 --- a/gammapy/datasets/make.py +++ b/gammapy/datasets/make.py @@ -6,7 +6,7 @@ import numpy as np from astropy.units import Quantity from astropy.time import Time, TimeDelta -from astropy.table import Table, Column +from astropy.table import Table from astropy.coordinates import SkyCoord, AltAz, FK5, Angle from ..irf import EnergyDependentMultiGaussPSF from ..obs import ObservationTable, observatory_locations @@ -110,26 +110,15 @@ def make_test_observation_table(observatory_name, n_obs, debug=False): # obs id obs_id = np.arange(n_obs_start, n_obs_start + n_obs) - col_obs_id = Column(name='OBS_ID', data=obs_id) - obs_table.add_column(col_obs_id) - - # TODO: `~astropy.table.Table` doesn't handle `~astropy.time.TimeDelta` columns - # properly see: - # https://github.com/astropy/astropy/issues/3832 - # until this issue is solved (and included into a stable release), quantity - # objects are used. + obs_table['OBS_ID'] = obs_id # obs time: 30 min - time_observation = TimeDelta(30.*60.*np.ones_like(col_obs_id.data), format='sec') - time_observation = Quantity(time_observation.sec, 'second') # converting to quantity - col_time_observation = Column(name='TIME_OBSERVATION', data=time_observation) - obs_table.add_column(col_time_observation) + time_observation = Quantity(30.*np.ones_like(obs_id), 'minute').to('second') + obs_table['TIME_OBSERVATION'] = time_observation # livetime: 25 min - time_live = TimeDelta(25.*60.*np.ones_like(col_obs_id.data), format='sec') - time_live = Quantity(time_live.sec, 'second') # converting to quantity - col_time_live = Column(name='TIME_LIVE', data=time_live) - obs_table.add_column(col_time_live) + 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 @@ -139,16 +128,19 @@ def make_test_observation_table(observatory_name, n_obs, debug=False): # 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(col_obs_id)) + datestart.mjd, format='mjd', 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 = TimeDelta(22.*60.*60., format='sec') - night_duration = TimeDelta(5.5*60.*60., format='sec') - hour_start = night_start + TimeDelta(night_duration.sec*np.random.random(len(col_obs_id)), format='sec') + # 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 @@ -159,10 +151,10 @@ def make_test_observation_table(observatory_name, n_obs, debug=False): else : # show the observation times in seconds after the reference time_start = time_relative_to_ref(time_start, header) - time_start = Quantity(time_start.sec, 'second') # converting to quantity + # converting to quantity (beter treatment of units) + time_start = Quantity(time_start.sec, 'second') - col_time_start = Column(name='TIME_START', data=time_start) - obs_table.add_column(col_time_start) + obs_table['TIME_START'] = time_start # stop time # calculated as TIME_START + TIME_OBSERVATION @@ -170,20 +162,18 @@ def make_test_observation_table(observatory_name, n_obs, debug=False): 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']) - time_stop = Quantity(time_stop.sec, 'second') # converting to quantity + # converting to quantity (beter treatment of units) + time_stop = Quantity(time_stop.sec, 'second') - col_time_stop = Column(name='TIME_STOP', data=time_stop) - obs_table.add_column(col_time_stop) + obs_table['TIME_STOP'] = time_stop # az, alt # random points in a sphere above 45 deg altitude - az, alt = sample_sphere(len(col_obs_id), (0, 360), (45, 90), 'degree') + az, alt = sample_sphere(len(obs_id), (0, 360), (45, 90), 'degree') az = Angle(az, 'degree') alt = Angle(alt, 'degree') - col_az = Column(name='AZ', data=az) - obs_table.add_column(col_az) - col_alt = Column(name='ALT', data=alt) - obs_table.add_column(col_alt) + obs_table['AZ'] = az + obs_table['ALT'] = alt # RA, dec # derive from az, alt taking into account that alt, az represent the values @@ -194,19 +184,18 @@ def make_test_observation_table(observatory_name, n_obs, debug=False): az = Angle(obs_table['AZ']) alt = Angle(obs_table['ALT']) if debug : - obstime = Time(obs_table['TIME_START']) + TimeDelta(obs_table['TIME_OBSERVATION'])/2. + obstime = Time(obs_table['TIME_START']) + obstime += 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. + obstime = time_ref_from_dict(obs_table.meta) + obstime += TimeDelta(obs_table['TIME_START']) + obstime += TimeDelta(obs_table['TIME_OBSERVATION'])/2. location = observatory_locations[observatory_name] alt_az_coord = AltAz(az = az, alt = alt, obstime = obstime, location = location) # optional: make it depend on other pars: temperature, pressure, humidity,... sky_coord = alt_az_coord.transform_to(FK5) - ra = sky_coord.ra - col_ra = Column(name='RA', data=ra) - obs_table.add_column(col_ra) - dec = sky_coord.dec - col_dec = Column(name='DEC', data=dec) - obs_table.add_column(col_dec) + obs_table['RA'] = sky_coord.ra + obs_table['DEC'] = sky_coord.dec # optional: it would be nice to plot a skymap with the simulated RA/dec positions @@ -214,16 +203,16 @@ def make_test_observation_table(observatory_name, n_obs, debug=False): # 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(col_obs_id)) - col_n_tels = Column(name='N_TELS', data=n_tels) - obs_table.add_column(col_n_tels) + 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 = (muon_efficiency_max - muon_efficiency_min)*np.random.random(len(col_obs_id)) + muon_efficiency_min - col_muon_efficiency = Column(name='MUON_EFFICIENCY', data=muon_efficiency) - obs_table.add_column(col_muon_efficiency) + muon_efficiency = np.random.random(len(obs_id)) + muon_efficiency *= (muon_efficiency_max - muon_efficiency_min) + muon_efficiency += muon_efficiency_min + obs_table['MUON_EFFICIENCY'] = muon_efficiency return obs_table diff --git a/gammapy/utils/time.py b/gammapy/utils/time.py index d8ea0a5f04..68c412250f 100644 --- a/gammapy/utils/time.py +++ b/gammapy/utils/time.py @@ -11,17 +11,19 @@ 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 + meta : `~dict` dictionary with the keywords `MJDREFI` and `MJDREFF` Returns ------- - `~astropy.time.Time` object with the reference time in format of `MJD`. + 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']) @@ -35,6 +37,7 @@ def time_ref_from_dict(meta): 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. From d1da484674a37c1b986b94942d7897d03930e815 Mon Sep 17 00:00:00 2001 From: mapaz Date: Thu, 11 Jun 2015 15:04:48 -0700 Subject: [PATCH 14/15] Added a few more tests to the dummy obs table generator. Added doc about times in Gammapy and a glossary. Fixed many small formatting isues. --- docs/dataformats/file_formats.rst | 4 +- docs/dataformats/index.rst | 10 ++++ docs/dataformats/observation_lists.rst | 67 +++++++++++++------------- docs/references.rst | 14 ++++++ gammapy/datasets/make.py | 61 ++++++++++++----------- gammapy/datasets/tests/test_make.py | 20 ++++++++ gammapy/obs/observation.py | 7 +-- gammapy/utils/tests/test_time.py | 10 ++-- 8 files changed, 116 insertions(+), 77 deletions(-) diff --git a/docs/dataformats/file_formats.rst b/docs/dataformats/file_formats.rst index 1f253f132f..72c63af002 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: + CSV +++ @@ -195,7 +197,7 @@ XML **XML** files (see `Wikipedia `__) GammaLib / ctools uses an "observation definition" XML format described -`here `__. +`here `__. TODO: describe diff --git a/docs/dataformats/index.rst b/docs/dataformats/index.rst index 57c9b15426..55eb16df61 100644 --- a/docs/dataformats/index.rst +++ b/docs/dataformats/index.rst @@ -19,6 +19,16 @@ 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 e00238f7d7..5d90330ed9 100644 --- a/docs/dataformats/observation_lists.rst +++ b/docs/dataformats/observation_lists.rst @@ -18,45 +18,46 @@ A run list must have at least a column called ``OBS_ID``:: Usually it has many more columns with information about each observation. A list of Gammapy supported columns is: -================ ================================================================================================ ========= -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 yes? -ALT average altitude angle during the observation yes? -MUON_EFFICIENCY average muon efficiency of the telescopes yes? -TIME_START start time of the observation stored as number of seconds after the time reference in the header yes? -TIME_STOP end time of the observation in the same forma as TIME_START yes? -TIME_OBSERVATION duration of the observation yes? -TIME_LIVE duration of the observation without dead time yes? -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" (not used by Gammapy at the moment) no -================ ================================================================================================ ========= - -Extra user defined columns are allowed; Gammapy will ignore them. +================ =========================================================================== ========= +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 +================ =========================================================================== ========= + +Extra user-defined columns are allowed; Gammapy will ignore them. In order for the extra columns to have full meaning the following is needed: * 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 -================ ========================================================================================================================================================================================================== -OBSERVATORY_NAME name of the observatory where the observations were taken. This is important for instance for coordinate transformations between celestial (i.e. RA/dec) and terrestrial (i.e. az/alt) coordinate systems. -MJDREFI reference time for other times in the list (i.e. TIME_START/TIME_STOP). Integer value in mean julian days. -MJDREFF fraction of integer value defined in MJDREFI. -================ ========================================================================================================================================================================================================== - -Extra user defined header entries are allowed; Gammapy will ignore them. - -TODO: ``TEL_LIST``: since the final format is not clear, please mind having a look to this issue and its outcome: -https://github.com/gammapy/gammapy/issues/282 +================ =========================================================================== ========= +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 diff --git a/docs/references.rst b/docs/references.rst index 64cbf273d9..e9ba1b55c4 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` +.. [MET] mission elapsed time; see also :ref:`time_gammapy`. diff --git a/gammapy/datasets/make.py b/gammapy/datasets/make.py index 33948dceb0..d42dd0f4c6 100644 --- a/gammapy/datasets/make.py +++ b/gammapy/datasets/make.py @@ -82,16 +82,16 @@ def make_test_observation_table(observatory_name, n_obs, debug=False): Parameters ---------- observatory_name : string - name of the observatory; a list of choices is given in `~gammapy.obs.observatory_locations` + 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 + number of observations for the obs table debug : bool - show UTC times instead of seconds after the reference + show UTC times instead of seconds after the reference Returns ------- obs_table : `~gammapy.obs.ObservationTable` - observation table + observation table """ n_obs_start = 1 @@ -99,8 +99,6 @@ def make_test_observation_table(observatory_name, n_obs, debug=False): # build a time reference as the start of 2010 dateref = Time('2010-01-01 00:00:00', format='iso', scale='utc') - # TODO: using format "iso" for `~astropy.Time` for now; eventually change it - # to "fits" after the next astropy stable release (>1.0) is out. dateref_mjd_fra, dateref_mjd_int = np.modf(dateref.mjd) # header @@ -113,11 +111,12 @@ def make_test_observation_table(observatory_name, n_obs, debug=False): obs_table['OBS_ID'] = obs_id # obs time: 30 min - time_observation = Quantity(30.*np.ones_like(obs_id), 'minute').to('second') + 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') + time_live = Quantity(25. * np.ones_like(obs_id), 'minute').to('second') obs_table['TIME_LIVE'] = time_live # start time @@ -128,9 +127,8 @@ def make_test_observation_table(observatory_name, n_obs, debug=False): # 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') + 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) @@ -140,15 +138,15 @@ def make_test_observation_table(observatory_name, n_obs, debug=False): # 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)) + 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 : + if debug: # show the observation times in UTC time_start = time_start.iso - else : + 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) @@ -158,10 +156,12 @@ def make_test_observation_table(observatory_name, n_obs, debug=False): # 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']) + 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') @@ -183,21 +183,21 @@ def make_test_observation_table(observatory_name, n_obs, debug=False): # in TIME_START and TIME_STOP az = Angle(obs_table['AZ']) alt = Angle(obs_table['ALT']) - if debug : - obstime = Time(obs_table['TIME_START']) - obstime += TimeDelta(obs_table['TIME_OBSERVATION'])/2. - else : - obstime = time_ref_from_dict(obs_table.meta) - obstime += TimeDelta(obs_table['TIME_START']) - obstime += TimeDelta(obs_table['TIME_OBSERVATION'])/2. + 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) - # optional: make it depend on other pars: temperature, pressure, humidity,... + alt_az_coord = AltAz(az=az, alt=alt, obstime=obstime, location=location) + # optional: make it depend on other pars: temperature, pressure, + # humidity,... sky_coord = alt_az_coord.transform_to(FK5) obs_table['RA'] = sky_coord.ra obs_table['DEC'] = sky_coord.dec - # optional: it would be nice to plot a skymap with the simulated RA/dec positions + # positions # number of telescopes # random integers between 3 and 4 @@ -210,9 +210,8 @@ def make_test_observation_table(observatory_name, n_obs, debug=False): # 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 *= (muon_efficiency_max - muon_efficiency_min) - muon_efficiency += muon_efficiency_min + 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 ec769d6bac..447124138c 100644 --- a/gammapy/datasets/tests/test_make.py +++ b/gammapy/datasets/tests/test_make.py @@ -22,4 +22,24 @@ 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 6f74805952..97ac78d961 100644 --- a/gammapy/obs/observation.py +++ b/gammapy/obs/observation.py @@ -14,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 @@ -59,10 +60,6 @@ def info(self): 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) - # TODO: for an unknown reason, I can't enable astropy.units.cds - # as explained in the following link - # http://astropy.readthedocs.org/en/latest/units/#module-astropy.units.cds - # for using a Quantity object with MJD units. return ss def select_linspace_subset(self, num): diff --git a/gammapy/utils/tests/test_time.py b/gammapy/utils/tests/test_time.py index 985845d690..bc3a9d4cd0 100644 --- a/gammapy/utils/tests/test_time.py +++ b/gammapy/utils/tests/test_time.py @@ -11,17 +11,13 @@ def test_time_ref_from_dict(): mjd_frac = 0.5 time_ref_dict = dict(MJDREFI=mjd_int, MJDREFF=mjd_frac) time_ref = time_ref_from_dict(time_ref_dict) - decimal = 4 - assert_almost_equal(time_ref.mjd, mjd_int + mjd_frac, decimal) + assert_almost_equal(time_ref.mjd, mjd_int + mjd_frac, decimal=4) def test_time_relative_to_ref(): - mjd_int = 500 - mjd_frac = 0.5 - time_ref_dict = dict(MJDREFI=mjd_int, MJDREFF=mjd_frac) + 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) - decimal = 4 - assert_almost_equal(delta_time.sec, delta_time_1sec.sec, decimal) + assert_almost_equal(delta_time.sec, delta_time_1sec.sec, decimal=4) From bccbc6b16092060448812351b4e936738c21cef0 Mon Sep 17 00:00:00 2001 From: mapaz Date: Thu, 11 Jun 2015 23:07:38 -0700 Subject: [PATCH 15/15] Fixed double reference warning from sphinx build. Improved format of sphinx and string docs. --- docs/dataformats/file_formats.rst | 2 +- docs/dataformats/index.rst | 9 +++++++-- docs/dataformats/observation_lists.rst | 6 ++++-- docs/references.rst | 2 +- gammapy/datasets/make.py | 2 -- gammapy/utils/time.py | 11 ++++++----- 6 files changed, 19 insertions(+), 13 deletions(-) diff --git a/docs/dataformats/file_formats.rst b/docs/dataformats/file_formats.rst index 72c63af002..57d21f9ba2 100644 --- a/docs/dataformats/file_formats.rst +++ b/docs/dataformats/file_formats.rst @@ -145,7 +145,7 @@ INI files ... to use it in your code Unfortunately INI files are not standardised, so there's only conventions and tons of variants. -.. _CSV: +.. _CSV_files: CSV +++ diff --git a/docs/dataformats/index.rst b/docs/dataformats/index.rst index 55eb16df61..5b03bd490f 100644 --- a/docs/dataformats/index.rst +++ b/docs/dataformats/index.rst @@ -25,8 +25,13 @@ in the :ref:`dataformats_file_formats` section. 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. +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:: diff --git a/docs/dataformats/observation_lists.rst b/docs/dataformats/observation_lists.rst index 5d90330ed9..02d593f96e 100644 --- a/docs/dataformats/observation_lists.rst +++ b/docs/dataformats/observation_lists.rst @@ -16,7 +16,8 @@ A run list must have at least a column called ``OBS_ID``:: 42 43 -Usually it has many more columns with information about each observation. A list of Gammapy supported columns is: +Usually it has many more columns with information about each observation. A list of +Gammapy supported columns is: ================ =========================================================================== ========= column name description required? @@ -62,7 +63,8 @@ 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. +The tool `~gammapy.datasets.make_test_observation_table` can generate a `~gammapy.obs.ObservationTable` +with dummy values. Header (metadata):: diff --git a/docs/references.rst b/docs/references.rst index e9ba1b55c4..1999b7cc3f 100644 --- a/docs/references.rst +++ b/docs/references.rst @@ -70,5 +70,5 @@ Software references: Glossary -------- -.. [CSV] comma-separated values; see also :ref:`CSV` +.. [CSV] comma-separated values; see also :ref:`CSV_files`. .. [MET] mission elapsed time; see also :ref:`time_gammapy`. diff --git a/gammapy/datasets/make.py b/gammapy/datasets/make.py index d42dd0f4c6..0293e0900c 100644 --- a/gammapy/datasets/make.py +++ b/gammapy/datasets/make.py @@ -191,8 +191,6 @@ def make_test_observation_table(observatory_name, n_obs, debug=False): 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) - # optional: make it depend on other pars: temperature, pressure, - # humidity,... sky_coord = alt_az_coord.transform_to(FK5) obs_table['RA'] = sky_coord.ra obs_table['DEC'] = sky_coord.dec diff --git a/gammapy/utils/time.py b/gammapy/utils/time.py index 68c412250f..0b0c12ee31 100644 --- a/gammapy/utils/time.py +++ b/gammapy/utils/time.py @@ -18,12 +18,12 @@ def time_ref_from_dict(meta): Parameters ---------- meta : `~dict` - dictionary with the keywords `MJDREFI` and `MJDREFF` + dictionary with the keywords `MJDREFI` and `MJDREFF` Returns ------- time : `~astropy.time.Time` - reference time with ``format='MJD'`` + 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']) @@ -44,13 +44,14 @@ def time_relative_to_ref(time, meta): Parameters ---------- time: `~astropy.time.Time` - time to be converted. + time to be converted. meta : dict - dictionary with the keywords `MJDREFI` and `MJDREFF` + dictionary with the keywords `MJDREFI` and `MJDREFF` Returns ------- - `~astropy.time.TimeDelta` object with the time in seconds after the reference. + 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')