Merge pull request #231 from cdeil/obs-locations

Observatory and data classes

docs/data/index.rst

@@ -11,20 +11,65 @@ Data classes (`gammapy.data`)
 Introduction
 ============
 
-`gammapy.data` contains classes to represent gamma-ray data:
+`gammapy.data` contains classes to represent gamma-ray data.
 
-* An `~gammapy.data.EventList` contains unbinned data,
-  i.e. the longitude, latitude, energy and time for each event.
-* A `~gammapy.data.SpectralCube` contains binned dat as a 3-dimensional array
-  with longitude, latitude and log energy axes.
-* The event time information has been removed ... it is represented by
-  a `~gammapy.data.GoodTimeInterval` (a.k.a. ``GTI``) class that is needed
-  for exposure and flux computations.
+We follow the Fermi data model and FITS formats as much as possible.
+
+The data format used is FITS and we define one container class to represent
+each FITS extension. In addition we define two high-level dataset classes
+that group all info data and metadata that's usually given for a single
+observation together:
+
+* Unbinned data is represented by a `~gammapy.data.EventListDataset`, which contains:
+
+  - `~gammapy.data.EventList` - table with time, position and energy for each event
+  - `~gammapy.data.GoodTimeIntervals` - table of good time intervals.
+    Used for livetime and exposure computation.
+  - `~gammapy.data.TelescopeArray` - some info about the array that took the data.
+    Optional, not used at the moment.
+
+* For most analysis the first step is to bin the data, which turns the
+  `~gammapy.data.EventListDataset` into one of:
+
+  - `~gammapy.data.CountsCubeDataset` (lon, lat, energy)
+  - `~gammapy.data.CountsSpectrumDataset` (energy)
+  - `~gammapy.data.CountsImageDataset` (lon, lat)
+  - `~gammapy.data.CountsLightCurveDataset` (time)
+
+* TODO: add IRFs to the dataset classes?
+
+* TODO: do we need the ``*Dataset`` wrappers or can we only have
+  ``EventList``, ``CountsCube``, ``CountsSpectrum``, ``CountsImage``, ...?
+
+* We'll have to see if we want to copy the IRF and / or GTI and / or TelescopeArray info
+  over to the binned dataset classes ... at the moment it's not clear if we need that info.
+
+Energy binning
+--------------
+
+* `~gammapy.spectrum.EnergyBinCenters` and FITS "ENERGIES" extensions.
+* `~gammapy.spectrum.EnergyBinEdges` and FITS "EBOUNDS" extensions.
+
+Spatial binning
+---------------
+
+TODO: Should we define a "SpatialGrid" or "SpatialBinning" class that wraps
+the 2d image FITS WCS and add convenience methods like generating them from scratch
+e.g. centered on the target or aligned to a survey map grid?
 
 Getting Started
 ===============
 
-TODO
+
+.. code-block:: python
+
+    >>> from gammapy.data import EventListDataset
+    >>> events_ds = EventListDataset.read('events.fits')
+    >>> events_ds.info()
+    >>> from gammapy.data import CountsCubeDataset
+    # This is just an idea ... not implemented!
+    >>> counts_ds = CountsCubeDataset.from_events(events_ds, ...)
+    >>> counts_ds = events_ds.make_counts_cube(...)
 
 Reference/API
 =============

docs/development/index.rst

@@ -159,3 +159,119 @@ Make a Gammapy release
 For now, see https://github.com/astropy/package-template/issues/103
 
 * Check external HTML links (see :ref:`here <development-check_html_links>`).
+
+Other codes
+-----------
+
+These projects are on Github, which is great because
+it has full-text search and git history view:
+
+* https://github.com/gammapy/gammapy
+* https://github.com/gammapy/gammapy-extra
+* https://github.com/astropy/astropy
+* https://github.com/astropy/photutils
+* https://github.com/gammalib/gammalib
+* https://github.com/ctools/ctools
+* https://github.com/zblz/naima
+* https://github.com/woodmd/gammatools
+* https://github.com/kialio/VHEObserverTools
+
+These are unofficial, unmaintained copies on open codes on Github:
+
+* https://github.com/brefsdal/sherpa
+* https://github.com/Xarthisius/yt-drone
+* https://github.com/cdeil/Fermi-ScienceTools-mirror
+* https://github.com/cdeil/kapteyn-mirror
+
+What checks and conversions should I do for inputs?
+---------------------------------------------------
+
+In Gammapy we assume that
+`"we're all consenting adults" <https://mail.python.org/pipermail/tutor/2003-October/025932.html>`_,
+which means that when you write a function you should write it like this:
+
+.. code-block:: python
+
+    def do_something(data, option):
+        """Do something.
+
+        Parameters
+        ----------
+        data : `numpy.ndarray`
+            Data
+        option : {'this', 'that'}
+            Option
+        """
+        if option == 'this':
+            out = 3 * data
+        elif option == 'that':
+            out = data ** 5
+        else:
+            ValueError('Invalid option: {}'.format(option))
+
+        return out
+
+* **Don't always add `isinstance` checks for everything** ... assume the caller passes valid inputs,
+  ... in the example above this is not needed::
+
+        assert isinstance(option, str)
+
+* **Don't always add `numpy.asanyarray` calls for all array-like inputs** ... the caller can do this if
+  it's really needed ... in the example above document ``data`` as type `~numpy.ndarray`
+  instead of array-like and don't put this line::
+
+        data = np.asanyarray(data)
+
+* **Do always add an `else` clause to your `if`-`elif` clauses** ... this is boilerplate code,
+  but not adding it would mean users get this error if they pass an invalid option::
+
+      UnboundLocalError: local variable 'out' referenced before assignment
+
+
+Now if you really want, you can add the `numpy.asanyarray` and `isinstance` checks
+for functions that end-users might often use for interactive work to provide them with
+better exception messages, but doing it everywhere would mean 1000s of lines of boilerplate
+code and take the fun out of Python programming.
+
+Float data type: 32 bit or 64 bit?
+----------------------------------
+
+Most of the time what we want is to use 32 bit to store data on disk and 64 bit to do
+computations in memory.
+
+Using 64 bit to store data and results (e.g. large images or cubes) on disk would mean
+a factor ~2 increase in file sizes and slower I/O, but I'm not aware of any case
+where we need that precision.
+
+On the other hand, doing computations with millions and billions of pixels very frequently
+results in inaccurate results ... e.g. the likelihood is the sum over per-pixel likelihoods
+and using 32-bit will usually result in erratic and hard-to-debug optimizer behaviour
+and even if the fit works incorrect results.
+
+Now you shouldn't put this line at the top of every function ... assume the caller
+passes 64-bit data::
+
+        data = np.asanyarray(data, dtype='float64')
+
+But you should add explicit type conversions to 64 bit when reading float data from files
+and explicit type conversions to 32 bit before writing to file.
+
+Clobber or overwrite?
+---------------------
+
+In Gammapy we use on ``overwrite`` bool option for `gammapy.scripts` and functions that
+write to files.
+
+Why not use ``clobber`` instead?
+After all the
+`FTOOLS <http://heasarc.gsfc.nasa.gov/ftools/ftools_menu.html>`__
+always use ``clobber``.
+
+The reason is that ``overwrite`` is clear to everyone, but ``clobber`` is defined by the dictionary
+(e.g. see `here <http://dictionary.reference.com/browse/clobber>`__)
+as "to batter severely; strike heavily. to defeat decisively. to denounce or criticize vigorously."
+and isn't intuitively clear to new users.
+
+Astropy uses both ``clobber`` and ``overwrite`` in various places at the moment.
+For Gammapy we can re-visit this decision before the 1.0 release, but for now,
+please be consistent and use ``overwrite``.

docs/obs/findruns.rst

@@ -0,0 +1,50 @@
+.. _obs_findruns:
+
+findruns
+========
+
+TODO: reference run list format
+
+The ``gammapy-findruns`` command line tool can be used to select a subset of runs from a given run list.
+
+Examples:
+
+* Find all runs within 5 deg of the Galactic center:
+
+  .. code-block:: bash
+
+      $ gammapy-findruns cone --x 0 --y 0 --r 5 --system galactic \
+                         --in all_runs.lis --out galactic_center_runs.lis
+
+* Select all runs in a box along the Galactic plane (GLON = -20 .. +20 deg, GLAT = -3 .. +3 deg):
+
+  .. code-block:: bash
+
+      $ gammapy-findruns box --x 0 --y 0 --dx 20 --dy 3 --system galactic \
+                         --in all_runs.lis  --out galactic_plane_runs.lis
+
+* Select all runs in a given date range (can of course be combined with other selections shown above):
+
+  .. code-block:: bash
+
+      $ gammapy-findruns --date_min 2010-04-26 --date_max "2010-04-29 12:42" \
+                         --in all_runs.lis --out run_042_to_100.lis
+
+* Select all runs in a given run range (can of course be combined with other selections shown above):
+
+  .. code-block:: bash
+
+      $ gammapy-findruns --run_min 42 --run_max 100 \
+                         --in all_runs.lis --out run_042_to_100.lis
+
+Using ``gammapy-findruns`` is easy, you don't have to remember all the options, just type:
+
+  .. code-block:: bash
+
+      $ findruns --help
+
+at the command line or read the usage help (TODO: add link here).
+
+.. note:: Noel Dawe's ``goodruns`` tool for `ATLAS <http://atlas.ch>`__ run selection
+      (`docs <http://ndawe.github.io/goodruns/>`__, `code <https://github.com/ndawe/goodruns>`__)
+      is a nice example for a run selection tool.

docs/obs/index.rst

@@ -18,58 +18,24 @@ for a given amount of time (e.g. half an hour) and switching the central trigger
 The total dataset for a given target will usually consist of a few to a few 100 runs
 and some book-keeping is required when running the analysis.
 
-
 Getting Started
 ===============
 
-TODO: reference run list format
-
-findruns
---------
-
-The ``gammapy-findruns`` command line tool can be used to select a subset of runs from a given run list.
-
-Examples:
-
-* Find all runs within 5 deg of the Galactic center:
-
-  .. code-block:: bash
-
-      $ gammapy-findruns cone --x 0 --y 0 --r 5 --system galactic \
-                         --in all_runs.lis --out galactic_center_runs.lis
-   
-* Select all runs in a box along the Galactic plane (GLON = -20 .. +20 deg, GLAT = -3 .. +3 deg):
-
-  .. code-block:: bash
-
-      $ gammapy-findruns box --x 0 --y 0 --dx 20 --dy 3 --system galactic \
-                         --in all_runs.lis  --out galactic_plane_runs.lis
-
-* Select all runs in a given date range (can of course be combined with other selections shown above):
-
-  .. code-block:: bash
-
-      $ gammapy-findruns --date_min 2010-04-26 --date_max "2010-04-29 12:42" \
-                         --in all_runs.lis --out run_042_to_100.lis
-
-* Select all runs in a given run range (can of course be combined with other selections shown above):
-
-  .. code-block:: bash
-
-      $ gammapy-findruns --run_min 42 --run_max 100 \
-                         --in all_runs.lis --out run_042_to_100.lis
-   
-Using ``gammapy-findruns`` is easy, you don't have to remember all the options, just type:
+.. code-block:: python
 
-  .. code-block:: bash
+   >>> from gammapy.obs import observatory_locations
+   >>> observatory_locations.HESS
+   <EarthLocation (7237.152530011689, 2143.7727767623487, -3229.3927009565496) km>
+   >>> print(observatory_locations.HESS.geodetic)
+   (<Longitude 16.500222222222224 deg>, <Latitude -23.271777777772456 deg>, <Quantity 1835.0 km>)
 
-      $ findruns --help
+Using `gammapy.obs`
+=====================
 
-at the command line or read the usage help (TODO: add link here).
+.. toctree::
+   :maxdepth: 1
 
-.. note:: Noel Dawe's ``goodruns`` tool for `ATLAS <http://atlas.ch>`__ run selection
-      (`docs <http://ndawe.github.io/goodruns/>`__, `code <https://github.com/ndawe/goodruns>`__)
-      is a nice example for a run selection tool.
+   findruns
 
 Reference/API
 =============

examples/add_gnomonic_projection.py

@@ -16,6 +16,8 @@
 import numpy as np
 
 
+
+
 def tan_world_to_pix(lon, lat, lon_center, lat_center):
     """Hand-coded TAN (a.k.a. tangential or Gnomonic) projection.
 

examples/eventlist_coordinate_check.py

@@ -0,0 +1,15 @@
+"""Eventlist coordinate check.
+"""
+from gammapy.data import EventListDataset
+from gammapy.datasets import get_path
+
+
+filename = get_path('hess/run_0023037_hard_eventlist.fits.gz')
+event_list = EventListDataset.read(filename)
+print(event_list.info)
+event_list.check()
+
+"""
+TODO: figure out the origin of this offset:
+ALT / AZ not consistent with RA / DEC. Max separation: 823.3407076612169 arcsec
+"""