New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Add function to compute counts maps #1317
Merged
Merged
Changes from all commits
Commits
Show all changes
14 commits
Select commit
Hold shift + click to select a range
6910834
First commit for new map/cube analysis scheme.
registerrier 819d319
Corrected small issues
registerrier 8d32b2a
Create fill_map_counts function and remove masking
registerrier 7c0c524
Missing line in doc string
registerrier e0fcb2e
Adding a type to MapAxis based on its unit
registerrier 2ae0cc3
Modified fill_map_counts to use the new maps.fill_by_coords
registerrier f0df042
now with the correct fill_by_coord. Also use the MapAxis.type
registerrier 6c5dab6
Some code implifications
registerrier b580d99
updated fill_map_counts
registerrier cffee8c
updated fill_map_counts with simpler treatment of column names
registerrier fb92cdd
added test for HpxMap and include function in __all__
registerrier 173d173
remove make_map_count for now
registerrier 1cc7d61
Minor cleanup in fill_map_counts
cdeil 995bab3
Update gammapy/cube/tests/test_basic_cube.py
cdeil File filter
Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -8,3 +8,4 @@ | |
from .utils import * | ||
from .models import * | ||
from .cube_pipe import * | ||
from .basic_cube import * |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,44 @@ | ||
# Licensed under a 3-clause BSD style license - see LICENSE.rst | ||
"""Functions to perform basic functions for map and cube analysis. | ||
""" | ||
from __future__ import absolute_import, division, print_function, unicode_literals | ||
|
||
__all__ = [ | ||
'fill_map_counts', | ||
] | ||
|
||
|
||
def fill_map_counts(count_map, event_list): | ||
"""Fill events into a counts map. | ||
|
||
The energy of the events is used for a non-spatial axis homogeneous to energy. | ||
The other non-spatial axis names should have an entry in the column names of the ``EventList`` | ||
|
||
Parameters | ||
---------- | ||
count_map : `~gammapy.maps.Map` | ||
Map object, will be filled by this function. | ||
event_list : `~gammapy.data.EventList` | ||
Event list | ||
""" | ||
geom = count_map.geom | ||
|
||
# Make a coordinate dictionary; skycoord is always added | ||
coord_dict = dict(skycoord=event_list.radec) | ||
|
||
# Now add one coordinate for each extra map axis | ||
for axis in geom.axes: | ||
if axis.type == 'energy': | ||
# This axis is the energy. We treat it differently because axis.name could be e.g. 'energy_reco' | ||
coord_dict[axis.name] = event_list.energy.to(axis.unit) | ||
# TODO: add proper extraction for time | ||
else: | ||
# We look for other axes name in the table column names (case insensitive) | ||
colnames = [_.upper() for _ in event_list.table.colnames] | ||
if axis.name.upper() in colnames: | ||
column_name = event_list.table.colnames[colnames.index(axis.name.upper())] | ||
coord_dict.update({axis.name: event_list.table[column_name].to(axis.unit)}) | ||
else: | ||
raise ValueError("Cannot find MapGeom axis {!r} in EventList".format(axis.name)) | ||
|
||
count_map.fill_by_coord(coord_dict) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,87 @@ | ||
# Licensed under a 3-clause BSD style license - see LICENSE.rst | ||
import pytest | ||
import numpy as np | ||
from astropy.coordinates import SkyCoord | ||
import astropy.units as u | ||
from ...utils.testing import requires_dependency, requires_data | ||
from ...maps import MapAxis, WcsGeom, HpxGeom, Map | ||
from ...data import DataStore, EventList | ||
from ...cube.basic_cube import fill_map_counts | ||
|
||
|
||
pytest.importorskip('scipy') | ||
|
||
# TODO: move these two cases to different test functions? | ||
# Would allow to change the asserts in `test_fill_map_counts` to something much more specific / useful, no? | ||
axis_energy_reco = MapAxis(np.logspace(-1., 1.5, 10), interp='log', node_type='edge', | ||
name='energy_reco', unit='TeV') | ||
|
||
# initial time of run 110380 from cta 1DC | ||
times = np.linspace(664504199, 664504199 + 1900., 10) | ||
axis_time = MapAxis(times, node_type='edge', name='time', unit='s') | ||
|
||
pos_cta = SkyCoord(0.0, 0.0, frame='galactic', unit='deg') | ||
|
||
geom_cta = {'binsz': 0.02, 'coordsys': 'GAL', 'width': 15 * u.deg, | ||
'skydir': pos_cta, 'axes': [axis_energy_reco]} | ||
geom_cta_time = {'binsz': 0.02, 'coordsys': 'GAL', 'width': 15 * u.deg, | ||
'skydir': pos_cta, 'axes': [axis_energy_reco, axis_time]} | ||
|
||
|
||
# TODO: change the test event list to something that's created from scratch, | ||
# using values so that it's possible to make simple assert statements on the | ||
# map data in the tests below, i.e. have pixels that should receive 0, 1 or 2 counts | ||
def make_test_event_list(): | ||
ds = DataStore.from_dir('$GAMMAPY_EXTRA/datasets/cta-1dc/index/gps/') | ||
return ds.obs(110380).events | ||
|
||
|
||
@pytest.fixture(scope='session') | ||
def evt_2fhl(): | ||
return EventList.read('$GAMMAPY_EXTRA/datasets/fermi_2fhl/2fhl_events.fits.gz') | ||
|
||
|
||
@requires_data('gammapy-extra') | ||
@pytest.mark.parametrize('geom_opts', [geom_cta, geom_cta_time]) | ||
def test_fill_map_counts(geom_opts): | ||
events = make_test_event_list() | ||
geom = WcsGeom.create(**geom_opts) | ||
cntmap = Map.from_geom(geom) | ||
|
||
fill_map_counts(cntmap, events) | ||
|
||
# Number of entries in the map | ||
nmap = cntmap.data.sum() | ||
# number of events | ||
valid = np.ones_like(events.energy.value) | ||
for axis in geom.axes: | ||
if axis.name == 'energy_reco': | ||
valid = np.logical_and(events.energy.value >= axis.edges[0], valid) | ||
valid = np.logical_and(events.energy.value <= axis.edges[-1], valid) | ||
else: | ||
valid = np.logical_and(events.table[axis.name.upper()] >= axis.edges[0], valid) | ||
valid = np.logical_and(events.table[axis.name.upper()] <= axis.edges[-1], valid) | ||
|
||
nevt = valid.sum() | ||
assert nmap == nevt | ||
|
||
|
||
@requires_data('gammapy-extra') | ||
@requires_dependency('healpy') | ||
def test_fill_map_counts_hpx(evt_2fhl): | ||
# This tests healpix maps fill with non standard non spatial axis | ||
|
||
axis_zen = MapAxis([0, 45, 180], node_type='edge', name='zenith_angle', unit=u.deg) | ||
geom = HpxGeom(256, coordsys='GAL', axes=[axis_zen]) | ||
m = Map.from_geom(geom) | ||
|
||
fill_map_counts(m, evt_2fhl) | ||
|
||
nmap_l = np.sum(m.data[0]) | ||
nmap_h = np.sum(m.data[1]) | ||
|
||
nevt_l = np.sum(evt_2fhl.table['ZENITH_ANGLE'] < 45) | ||
nevt_h = np.sum(evt_2fhl.table['ZENITH_ANGLE'] > 45) | ||
|
||
assert nmap_l == nevt_l | ||
assert nmap_h == nevt_h |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Please add a test for this type, or ideally even one that covers the three cases of "energy", "time" and "any".