Merge 16f072b into 4d935ae

gammapy/background/tests/test_reflected.py

@@ -2,15 +2,14 @@
 from __future__ import print_function, division
 import unittest
 from astropy.tests.helper import pytest
-from ...obs import RunList
 from ...background import Maps, ReflectedRegionMaker
 
 
 @pytest.mark.xfail
 class TestReflectedBgMaker(unittest.TestCase):
 
     def test_analysis(self):
-        runs = RunList('runs.lis')
+        runs = 'TODO'
         maps = Maps('maps.fits')
         reflected_bg_maker = ReflectedRegionMaker(runs, maps, psi=2, theta=0.1)
         total_maps = Maps('total_maps.fits')

gammapy/catalog/utils.py

@@ -246,17 +246,18 @@ def skycoord_from_table(table):
     """
 
     if set(['RAJ2000', 'DEJ2000']).issubset(table.colnames):
-        lon, lat, frame = 'RAJ2000', 'DEJ2000', 'fk5'
+        lon, lat, frame = 'RAJ2000', 'DEJ2000', 'icrs'
     elif set(['RA', 'DEC']).issubset(table.colnames):
-        lon, lat, frame = 'RA', 'DEC', 'galactic'
+        lon, lat, frame = 'RA', 'DEC', 'icrs'
     elif set(['GLON', 'GLAT']).issubset(table.colnames):
         lon, lat, frame = 'GLON', 'GLAT', 'galactic'
     else:
-        raise KeyError('No column GLON / GLAT or RAJ2000 / DEJ2000 found.')
+        raise KeyError('No column GLON / GLAT or RA / DEC or RAJ2000 / DEJ2000 found.')
 
     unit = table[lon].unit if table[lon].unit else 'deg'
 
     skycoord = SkyCoord(table[lon], table[lat], unit=unit, frame=frame)
+
     return skycoord
 
 

gammapy/obs/__init__.py

@@ -3,4 +3,5 @@
 Observation utility functions and classes
 """
 from .observers import *
-from .run import *
+from .observation import *
+from .datastore import *

gammapy/obs/datastore.py

@@ -0,0 +1,227 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+import os
+from astropy.table import Table
+from astropy.coordinates import SkyCoord, Angle
+from ..catalog import select_sky_box, skycoord_from_table
+from ..obs import ObservationTable
+
+__all__ = ['DataStore', 'DataStoreIndexTable']
+
+
+def _make_filename_hess_scheme(obs_id, filetype='events'):
+    """Make filename string for the HESS storage scheme.
+
+    Parameters
+    ----------
+    obs_id : int
+        Observation ID
+    filetype : {'events', 'effective area', 'psf', 'background'}
+        Type of file
+
+    Examples
+    --------
+    >>> _make_filename_hess_scheme(obs_id=89565, filetype='effective area')
+    'run089400-089599/run089565/hess_aeff_2d_089565.fits.gz'
+    """
+    obs_id_min = obs_id - (obs_id % 200)
+    obs_id_max = obs_id_min + 199
+    group_folder = 'run{:06d}-{:06d}'.format(obs_id_min, obs_id_max)
+    obs_folder = 'run{:06d}'.format(obs_id)
+
+    if filetype == 'events':
+        label = 'events'
+    elif filetype == 'psf':
+        label = 'psf_king'
+    elif filetype == 'effective area':
+        label = 'aeff_2d'
+    elif filetype == 'background':
+        label = 'bkg_offruns'
+    else:
+        raise ValueError('Unknown filetype: {}'.format(filetype))
+
+    filename = 'hess_{}_{:06d}.fits.gz'.format(label, obs_id)
+
+    return os.path.join(group_folder, obs_folder, filename)
+
+
+class DataStoreIndexTable(Table):
+    """Data store index table.
+
+    The index table is a FITS file that stores which observations
+    are available and what their most important parameters are.
+
+    This makes it possible to select observation of interest and find out
+    what data is available without opening up thousands of FITS files
+    that contain the event list and IRFs and have similar information in the
+    FITS header.
+    """
+    # For now I've decided to not do the cleanup in `__init__`,
+    # but instead in `read`.
+    # See https://groups.google.com/d/msg/astropy-dev/0EaOw9peWSk/MSjH7q_7htoJ
+    # def __init__(self, *args, **kwargs):
+    #     super(DataStoreIndexTable, self).__init__(*args, **kwargs)
+    #     self._fixes()
+
+    @staticmethod
+    def read(*args, **kwargs):
+        """Read from FITS file. See `~astropy.table.Table.read`."""
+        table = Table.read(*args, **kwargs)
+        table = DataStoreIndexTable(table)
+        table._init_cleanup()
+        return table
+
+    def _init_cleanup(self):
+        # Rename to canonical column names
+        renames = [('RA_PNT', 'RA'),
+                   ('DEC_PNT', 'DEC')
+                   ]
+        for name, new_name in renames:
+            self.rename_column(name, new_name)
+
+        # Add useful extra columns
+        if not set(['GLON', 'GLAT']).issubset(self.colnames):
+            skycoord = skycoord_from_table(self).galactic
+            self['GLON'] = skycoord.l.degree
+            self['GLAT'] = skycoord.b.degree
+
+    def info(self):
+        ss = 'Data store index table summary:\n'
+        ss += 'Number of observations: {}\n'.format(len(self))
+        obs_id = self['OBS_ID']
+        ss += 'Observation IDs: {} to {}\n'.format(obs_id.min(), obs_id.max())
+        ss += 'Dates: {} to {}\n'.format('TODO', 'TODO')
+        return ss
+
+    @property
+    def radec(self):
+        """ICRS sky coordinates (`~astropy.coordinates.SkyCoord`)"""
+        return SkyCoord(self['RA'], self['DEC'], unit='deg', frame='icrs')
+
+    @property
+    def galactic(self):
+        """Galactic sky coordinates (`~astropy.coordinates.SkyCoord`)"""
+        return SkyCoord(self['GLON'], self['GLAT'], unit='deg', frame='galactic')
+
+
+class DataStore(object):
+    """Data store - convenient way to access and select data.
+
+    TODO: add methods to sync with remote datastore...
+
+    Parameters
+    ----------
+    dir : `str`
+        Data store directory on user machine
+    scheme : {'hess'}
+        Scheme of organising and naming the files.
+    """
+
+    def __init__(self, dir, scheme='hess'):
+        self.dir = dir
+        self.index_table_filename = 'runinfo.fits'
+        filename = os.path.join(dir, self.index_table_filename)
+        print('Reading {}'.format(filename))
+        self.index_table = DataStoreIndexTable.read(filename)
+        self.scheme = scheme
+
+    def info(self):
+        """Summary info string."""
+        ss = 'Data store summary info:\n'
+        ss += 'Directory: {}\n'.format(self.dir)
+        ss += 'Index table: {}\n'.format(self.index_table_filename)
+        ss += 'Scheme: {}\n'.format(self.scheme)
+        ss += self.index_table.info()
+        return ss
+
+    def make_table_of_files(self):
+        """Make list of all files in the datastore directory.
+
+        Returns
+        -------
+        table : `~astropy.table.Table`
+            Table summarising info about files.
+        """
+        raise NotImplementedError
+
+    def make_summary_plots(self):
+        """Make some plots summarising the available observations.
+
+        E.g. histograms of time, run duration, muon efficiency, zenith angle, ...
+        Aitoff plot showing run locations.
+        """
+        raise NotImplementedError
+
+    def filename(self, obs_id, filetype='events'):
+        """File name (relative to datastore `dir`).
+
+        Parameters
+        ----------
+        obs_id : int
+            Observation ID
+        filetype : {'events', 'effective area', 'psf', 'background'}
+            Type of file
+
+        Returns
+        -------
+        filename : str
+            Filename (including path relative to datastore `dir`)
+        """
+        scheme = self.scheme
+
+        if scheme == 'hess':
+            return _make_filename_hess_scheme(obs_id, filetype)
+        else:
+            raise ValueError('Invalid scheme: {}'.format(scheme))
+
+    def make_run_list(self, selection):
+        """Make a run list, applying some selection.
+
+        TODO: implement a more flexible scheme to make box cuts
+        on any fields (including e.g. OBSID or TIME
+        Not sure what a simple yet powerful method to implement this is!?
+
+        Example
+        -------
+        >>> selection = dict(shape='box', frame='galactic',
+        ...                  lon=(-100, 50), lat=(-5, 5), border=2)
+        >>> run_list = data_store.make_run_list(selection)
+        """
+        selection_region_shape = selection['shape']
+
+        if selection_region_shape == 'box':
+            lon = selection['lon']
+            lat = selection['lat']
+            border = selection['border']
+            lon = Angle([lon[0] - border, lon[1] + border], 'deg')
+            lat = Angle([lat[0] - border, lat[1] + border], 'deg')
+            # print(lon, lat)
+            table = select_sky_box(self.index_table,
+                                   lon_lim=lon, lat_lim=lat,
+                                   frame=selection['frame'])
+        else:
+            raise ValueError('Invalid selection type: {}'.format(selection_region_shape))
+
+        return table
+
+
+def make_count_image(header, data_store, run_list):
+    """Make count image from event lists (like gtbin).
+
+    TODO: what's a good API and location for this?
+
+    Parameters
+    ----------
+    header : `~astropy.io.fits.Header`
+        FITS header
+    data_store : `~gammapy.obs.DataStore`
+        Data store to access the event lists
+    run_list : `~astropy
+
+    """
+    shape = (header['NAXIS2'], header['NAXIS1'])
+    data = np.zeros(shape, dtype='int')
+    for run in run_list:
+        filename = data_store.filename(run['OBS_ID'], filetype='events')
+        event_list = Event
+
+    return fits.ImageHDU(data=data, header=header)

gammapy/obs/run.py → gammapy/obs/observation.py

@@ -1,13 +1,15 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
-"""Run and RunList class"""
 from __future__ import (absolute_import, division, print_function,
                         unicode_literals)
+from astropy.table import Table
 
-__all__ = ['Run', 'RunList']
+__all__ = ['Observation', 'ObservationTable']
 
 
-class Run(object):
-    """Run parameters container.
+class Observation(object):
+    """Observation (a.k.a. Run).
+
+    TODO: not clear if this class is useful.
 
     Parameters
     ----------
@@ -28,6 +30,6 @@ def wcs_header(self, system='FOV'):
         raise NotImplementedError
 
 
-class RunList(list):
+class ObservationTable(Table):
     """Run list container.
     """

gammapy/obs/tests/test_datastore.py

@@ -0,0 +1,9 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+from __future__ import print_function, division
+from ...obs import DataStore
+
+
+def test_DataStore():
+    # data_store = DataStore(dir='TODO')
+    # TODO: add tests
+    assert True

gammapy/obs/tests/test_observation.py

@@ -0,0 +1,11 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+from __future__ import print_function, division
+from ...obs import Observation, ObservationTable
+
+
+def test_Observation():
+    Observation(GLON=42, GLAT=43)
+
+
+def test_ObservationTable():
+    ObservationTable()

gammapy/scripts/__init__.py

@@ -13,7 +13,7 @@
 from .cwt import *
 from .derived_images import *
 from .detect import *
-from .find_runs import *
+from .find_obs import *
 from .image_decompose_a_trous import *
 from .info import *
 from .irf_info import *

gammapy/scripts/find_runs.py → gammapy/scripts/find_obs.py

@@ -5,11 +5,11 @@
 import argparse
 from ..utils.scripts import get_parser
 
-__all__ = ['find_runs']
+__all__ = ['find_obs']
 
 
 def main(args=None):
-    parser = get_parser(find_runs)
+    parser = get_parser(find_obs)
     parser.add_argument('infile', type=argparse.FileType('r'),
                         help='Input run list file name')
     parser.add_argument('outfile', nargs='?', type=argparse.FileType('w'),
@@ -25,16 +25,16 @@ def main(args=None):
     parser.add_argument('--overwrite', action='store_true',
                         help='Overwrite existing output file?')
     args = parser.parse_args(args)
-    find_runs(**vars(args))
+    find_obs(**vars(args))
 
 
-def find_runs(infile,
-              outfile,
-              x,
-              y,
-              pix,
-              overwrite):
-    """Select a subset of runs from a given run list.
+def find_obs(infile,
+             outfile,
+             x,
+             y,
+             pix,
+             overwrite):
+    """Select a subset of observations from a given observation list.
 
     TODO: explain.
     """

setup.py

@@ -86,7 +86,7 @@
     'gammapy-cwt = gammapy.scripts.cwt:main',
     'gammapy-derived-images = gammapy.scripts.derived_images:main',
     'gammapy-detect = gammapy.scripts.detect:main',
-    'gammapy-find-runs = gammapy.scripts.find_runs:main',
+    'gammapy-find-obs = gammapy.scripts.find_obs:main',
     'gammapy-image-decompose-a-trous = gammapy.scripts.image_decompose_a_trous:main',
     'gammapy-info = gammapy.scripts.info:main',
     'gammapy-irf-info = gammapy.scripts.irf_info:main',