# Licensed under a 3-clause BSD style license - see LICENSE.rst
from __future__ import absolute_import, division, print_function, unicode_literals
from astropy.units import Quantity
__all__ = ["fill_map_counts"]
[docs]def fill_map_counts(counts_map, events):
"""Fill events into a counts map.
This method handles sky coordinates automatically.
For all other map axes, an event list table column with the same name
has to be present (case-insensitive compare), or a KeyError will be raised.
Note that this function is just a thin wrapper around the sky map
``fill_by_coord`` method. For more complex scenarios, use that directly.
Parameters
----------
counts_map : `~gammapy.maps.Map`
Map object, will be filled by this function.
events : `~gammapy.data.EventList`
Event list
Examples
--------
To make a counts map, create an empty map with a geometry of your choice
and then fill it using this function::
from gammapy.maps import Map
from gammapy.data import EventList
from gammapy.cube import fill_map_counts
events = EventList.read('$GAMMAPY_DATA/cta-1dc/data/baseline/gps/gps_baseline_110380.fits')
counts = Map.create(coordsys='GAL', skydir=(0, 0), binsz=0.1, npix=(120, 100))
fill_map_counts(counts, events)
counts.plot()
If you have a given map already, and want to make a counts image
with the same geometry (not using the pixel data from the original map), do this:
from gammapy.maps import Map
from gammapy.data import EventList
from gammapy.cube import fill_map_counts
events = EventList.read('$GAMMAPY_DATA/fermi_2fhl/2fhl_events.fits.gz')
reference_map = Map.read('$GAMMAPY_DATA/fermi_2fhl/fermi_2fhl_gc.fits.gz')
counts = Map.from_geom(reference_map.geom)
fill_map_counts(counts, events)
counts.smooth(3).plot()
It works for IACT and Fermi-LAT events, for WCS or HEALPix map geometries,
and also for extra axes. Especially energy axes are automatically handled correctly.
"""
coord = dict(skycoord=events.radec)
cols = {k.upper(): v for k, v in events.table.columns.items()}
for axis in counts_map.geom.axes:
try:
col = cols[axis.name.upper()]
coord[axis.name] = Quantity(col).to(axis.unit)
except KeyError:
raise KeyError("Column not found in event list: {!r}".format(axis.name))
counts_map.fill_by_coord(coord)