Source code for gammapy.data.observations

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

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

log = logging.getLogger(__name__)


[docs]class ObservationCTA(object): """Container class for an CTA observation Parameters follow loosely the "Required columns" here: http://gamma-astro-data-formats.readthedocs.io/en/latest/data_storage/obs_index/index.html#required-columns Parameters ---------- obs_id : int Observation ID gti : `~gammapy.data.GTI` Good Time Intervals events : `~gammapy.data.EventList` Event list aeff : `~gammapy.irf.EffectiveAreaTable2D` Effective area edisp : `~gammapy.irf.EnergyDispersion2D` Energy dispersion psf : `~gammapy.irf.PSF3D` or `~gammapy.irf.EnergyDependentMultiGaussPSF` or `~gammapy.irf.PSFKing` Tabled Point Spread Function bkg : `~gammapy.irf.Background2D` or `~gammapy.irf.Background3D` Background rate pointing_radec : `~astropy.coordinates.SkyCoord` Pointing RA / DEC sky coordinates observation_live_time_duration : `~astropy.units.Quantity` Live-time duration in seconds The dead-time-corrected observation time. Computed as ``t_live = t_observation * (1 - f_dead)`` where ``f_dead`` is the dead-time fraction. observation_dead_time_fraction : float Dead-time fraction Defined as dead-time over observation time. Dead-time is defined as the time during the observation where the detector did not record events: https://en.wikipedia.org/wiki/Dead_time https://adsabs.harvard.edu/abs/2004APh....22..285F The dead-time fraction is used in the live-time computation, which in turn is used in the exposure and flux computation. meta : `~collections.OrderedDict` Dictionary to store metadata Examples -------- A minimal working example of how to create an observation from CTA's 1DC is given in examples/example_observation_cta.py """ def __init__( self, obs_id=None, gti=None, events=None, aeff=None, edisp=None, psf=None, bkg=None, pointing_radec=None, observation_live_time_duration=None, observation_dead_time_fraction=None, meta=None, ): self.obs_id = obs_id self.gti = gti self.events = events self.aeff = aeff self.edisp = edisp self.psf = psf self.bkg = bkg self.pointing_radec = pointing_radec self.observation_live_time_duration = observation_live_time_duration self.observation_dead_time_fraction = observation_dead_time_fraction self.meta = meta or OrderedDict() def __str__(self): ss = "Info for OBS_ID = {}\n".format(self.obs_id) ss += "- Pointing pos: RA {:.2f} / Dec {:.2f}\n".format( self.pointing_radec.ra if self.pointing_radec else "None", self.pointing_radec.dec if self.pointing_radec else "None", ) tstart = np.atleast_1d(self.gti.time_start.fits)[0] if self.gti else "None" ss += "- Start time: {}\n".format(tstart) ss += "- Observation duration: {}\n".format( self.gti.time_sum if self.gti else "None" ) ss += "- Dead-time fraction: {:5.3f} %\n".format( 100 * self.observation_dead_time_fraction ) return ss
[docs]class DataStoreObservation(object): """IACT data store observation. Parameters ---------- obs_id : int Observation ID data_store : `~gammapy.data.DataStore` Data store obs_filter : `~gammapy.data.ObservationFilter`, optional Filter for the observation """ def __init__(self, obs_id, data_store, obs_filter=None): # Assert that `obs_id` is available if obs_id not in data_store.obs_table["OBS_ID"]: raise ValueError("OBS_ID = {} not in obs index table.".format(obs_id)) if obs_id not in data_store.hdu_table["OBS_ID"]: raise ValueError("OBS_ID = {} not in HDU index table.".format(obs_id)) self.obs_id = obs_id self.data_store = data_store self.obs_filter = obs_filter or ObservationFilter() def __str__(self): ss = "Info for OBS_ID = {}\n".format(self.obs_id) ss += "- Start time: {:.2f}\n".format(self.tstart.mjd) ss += "- Pointing pos: RA {:.2f} / Dec {:.2f}\n".format( self.pointing_radec.ra, self.pointing_radec.dec ) ss += "- Observation duration: {}\n".format(self.observation_time_duration) ss += "- Dead-time fraction: {:5.3f} %\n".format( 100 * self.observation_dead_time_fraction ) # TODO: Which target was observed? # TODO: print info about available HDUs for this observation ... return ss
[docs] def location(self, hdu_type=None, hdu_class=None): """HDU location object. Parameters ---------- hdu_type : str HDU type (see `~gammapy.data.HDUIndexTable.VALID_HDU_TYPE`) hdu_class : str HDU class (see `~gammapy.data.HDUIndexTable.VALID_HDU_CLASS`) Returns ------- location : `~gammapy.data.HDULocation` HDU location """ return self.data_store.hdu_table.hdu_location( obs_id=self.obs_id, hdu_type=hdu_type, hdu_class=hdu_class )
[docs] def load(self, hdu_type=None, hdu_class=None): """Load data file as appropriate object. Parameters ---------- hdu_type : str HDU type (see `~gammapy.data.HDUIndexTable.VALID_HDU_TYPE`) hdu_class : str HDU class (see `~gammapy.data.HDUIndexTable.VALID_HDU_CLASS`) Returns ------- object : object Object depends on type, e.g. for `events` it's a `~gammapy.data.EventList`. """ location = self.location(hdu_type=hdu_type, hdu_class=hdu_class) return location.load()
@property def events(self): """Load `gammapy.data.EventList` object and apply the filter.""" events = self.load(hdu_type="events") return self.obs_filter.filter_events(events) @property def gti(self): """Load `gammapy.data.GTI` object and apply the filter.""" try: gti = self.load(hdu_type="gti") except IndexError: # For now we support data without GTI HDUs # TODO: if GTI becomes required, we should drop this case # CTA discussion in https://github.com/open-gamma-ray-astro/gamma-astro-data-formats/issues/20 # Added in Gammapy in https://github.com/gammapy/gammapy/pull/1908 gti = self.data_store.obs_table.create_gti(obs_id=self.obs_id) return self.obs_filter.filter_gti(gti) @property def aeff(self): """Load effective area object.""" return self.load(hdu_type="aeff") @property def edisp(self): """Load energy dispersion object.""" return self.load(hdu_type="edisp") @property def psf(self): """Load point spread function object.""" return self.load(hdu_type="psf") @property def bkg(self): """Load background object.""" return self.load(hdu_type="bkg") @property def obs_info(self): """Observation information (`~collections.OrderedDict`).""" row = self.data_store.obs_table.select_obs_id(obs_id=self.obs_id)[0] return table_row_to_dict(row) @property def tstart(self): """Observation start time (`~astropy.time.Time`).""" met_ref = time_ref_from_dict(self.data_store.obs_table.meta) met = Quantity(self.obs_info["TSTART"].astype("float64"), "second") time = met_ref + met return time @property def tstop(self): """Observation stop time (`~astropy.time.Time`).""" met_ref = time_ref_from_dict(self.data_store.obs_table.meta) met = Quantity(self.obs_info["TSTOP"].astype("float64"), "second") time = met_ref + met return time @property def observation_time_duration(self): """Observation time duration in seconds (`~astropy.units.Quantity`). The wall time, including dead-time. """ return self.gti.time_sum @property def observation_live_time_duration(self): """Live-time duration in seconds (`~astropy.units.Quantity`). The dead-time-corrected observation time. Computed as ``t_live = t_observation * (1 - f_dead)`` where ``f_dead`` is the dead-time fraction. """ return self.gti.time_sum * (1 - self.observation_dead_time_fraction) @property def observation_dead_time_fraction(self): """Dead-time fraction (float). Defined as dead-time over observation time. Dead-time is defined as the time during the observation where the detector didn't record events: https://en.wikipedia.org/wiki/Dead_time https://adsabs.harvard.edu/abs/2004APh....22..285F The dead-time fraction is used in the live-time computation, which in turn is used in the exposure and flux computation. """ return 1 - self.obs_info["DEADC"] @property def pointing_radec(self): """Pointing RA / DEC sky coordinates (`~astropy.coordinates.SkyCoord`).""" lon, lat = self.obs_info["RA_PNT"], self.obs_info["DEC_PNT"] return SkyCoord(lon, lat, unit="deg", frame="icrs") @property def pointing_altaz(self): """Pointing ALT / AZ sky coordinates (`~astropy.coordinates.SkyCoord`).""" alt, az = self.obs_info["ALT_PNT"], self.obs_info["AZ_PNT"] return SkyCoord(az, alt, unit="deg", frame="altaz") @property def pointing_zen(self): """Pointing zenith angle sky (`~astropy.units.Quantity`).""" return Quantity(self.obs_info["ZEN_PNT"], unit="deg") @property def target_radec(self): """Target RA / DEC sky coordinates (`~astropy.coordinates.SkyCoord`).""" lon, lat = self.obs_info["RA_OBJ"], self.obs_info["DEC_OBJ"] return SkyCoord(lon, lat, unit="deg", frame="icrs") @property def observatory_earth_location(self): """Observatory location (`~astropy.coordinates.EarthLocation`).""" return earth_location_from_dict(self.obs_info) @property def muoneff(self): """Observation muon efficiency.""" return self.obs_info["MUONEFF"]
[docs] def peek(self): """Quick-look plots in a few panels.""" raise NotImplementedError
[docs] def to_observation_cta(self): """Convert to `~gammapy.data.ObservationCTA`. This loads all observation-related info from disk and stores it in the in-memory ``ObservationCTA``. Returns ------- observation : `~gammapy.data.ObservationCTA` Observation """ # maps the ObservationCTA class attributes to the DataStoreObservation properties props = { "obs_id": "obs_id", "gti": "gti", "events": "events", "aeff": "aeff", "edisp": "edisp", "psf": "psf", "bkg": "bkg", "pointing_radec": "pointing_radec", "observation_live_time_duration": "observation_live_time_duration", "observation_dead_time_fraction": "observation_dead_time_fraction", } for obs_cta_kwarg, ds_obs_prop in props.items(): try: props[obs_cta_kwarg] = getattr(self, ds_obs_prop) except Exception as exception: log.warning(exception) props[obs_cta_kwarg] = None return ObservationCTA(**props)
[docs] def check(self, checks="all"): """Run checks. This is a generator that yields a list of dicts. """ checker = ObservationChecker(self) return checker.run(checks=checks)
[docs]class Observations(object): """Container class that holds a list of observations. Parameters ---------- obs_list : list A list of `~gammapy.data.DataStoreObservation` """ def __init__(self, obs_list=None): self._obs_list = obs_list or [] def __getitem__(self, key): return self._obs_list[key] def __len__(self): return len(self._obs_list) def __str__(self): s = self.__class__.__name__ + "\n" s += "Number of observations: {}\n".format(len(self)) for obs in self: s += str(obs) return s
class ObservationChecker(Checker): """Check an observation. Checks data format and a bit about the content. """ CHECKS = { "events": "check_events", "gti": "check_gti", "aeff": "check_aeff", "edisp": "check_edisp", "psf": "check_psf", } def __init__(self, observation): self.observation = observation def _record(self, level="info", msg=None): return {"level": level, "obs_id": self.observation.obs_id, "msg": msg} def check_events(self): yield self._record(level="debug", msg="Starting events check") try: events = self.observation.load("events") except Exception: yield self._record(level="warning", msg="Loading events failed") return for record in EventListChecker(events).run(): yield record # TODO: split this out into a GTIChecker def check_gti(self): yield self._record(level="debug", msg="Starting gti check") try: gti = self.observation.load("gti") except Exception: yield self._record(level="warning", msg="Loading GTI failed") return if len(gti.table) == 0: yield self._record(level="error", msg="GTI table has zero rows") columns_required = ["START", "STOP"] for name in columns_required: if name not in gti.table.colnames: yield self._record( level="error", msg="Missing table column: {!r}".format(name) ) # TODO: Check that header keywords agree with table entries # TSTART, TSTOP, MJDREFI, MJDREFF # Check that START and STOP times are consecutive # times = np.ravel(self.table['START'], self.table['STOP']) # # TODO: not sure this is correct ... add test with a multi-gti table from Fermi. # if not np.all(np.diff(times) >= 0): # yield 'GTIs are not consecutive or sorted.' # TODO: add reference times for all instruments and check for this # Use TELESCOP header key to check which instrument it is. def _check_times(self): """Check if various times are consistent. The headers and tables of the FITS EVENTS and GTI extension contain various observation and event time information. """ # http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_Data/Time_in_ScienceTools.html # https://hess-confluence.desy.de/confluence/display/HESS/HESS+FITS+data+-+References+and+checks#HESSFITSdata-Referencesandchecks-Time telescope_met_refs = OrderedDict( FERMI=Time("2001-01-01T00:00:00"), HESS=Time("2001-01-01T00:00:00") ) meta = self.dset.event_list.table.meta telescope = meta["TELESCOP"] if telescope in telescope_met_refs.keys(): dt = self.time_ref - telescope_met_refs[telescope] if dt > self.accuracy["time"]: yield self._record( level="error", msg="Reference time incorrect for telescope" ) def check_aeff(self): yield self._record(level="debug", msg="Starting aeff check") try: aeff = self.observation.load("aeff") except Exception: yield self._record(level="warning", msg="Loading aeff failed") return # Check that thresholds are meaningful for aeff if ( "LO_THRES" in aeff.meta and "HI_THRES" in aeff.meta and aeff.meta["LO_THRES"] >= aeff.meta["HI_THRES"] ): yield self._record( level="error", msg="LO_THRES >= HI_THRES in effective area meta data" ) # Check that data isn't all null if np.max(aeff.data.data) <= 0: yield self._record( level="error", msg="maximum entry of effective area is <= 0" ) def check_edisp(self): yield self._record(level="debug", msg="Starting edisp check") try: edisp = self.observation.load("edisp") except Exception: yield self._record(level="warning", msg="Loading edisp failed") return # Check that data isn't all null if np.max(edisp.data.data) <= 0: yield self._record(level="error", msg="maximum entry of edisp is <= 0") def check_psf(self): yield self._record(level="debug", msg="Starting psf check") try: psf = self.observation.load("psf") except Exception: yield self._record(level="warning", msg="Loading psf failed") return # TODO: implement some basic check # The following doesn't work, because EnergyDependentMultiGaussPSF # has no attribute `data` # Check that data isn't all null # if np.max(psf.data.data) <= 0: # yield self._record( # level="error", msg="maximum entry of psf is <= 0" # )