Fixed TS map boundary handling

CHANGES.rst

@@ -24,13 +24,14 @@ For plans and progress see https://github.com/gammapy/gammapy/milestones/0.4
 Contributors
 ++++++++++++
 
-- ...
+- Axel Donath
 
 Pull requests
 +++++++++++++
 
 - Make background cube models [#319] (Manuel Paz Arribas)
--  Production of true/reco bg cube models should use the same model. [#335] (Manuel Paz Arribas)
+- Production of true/reco bg cube models should use the same model. [#335] (Manuel Paz Arribas)
+- Fix TS map boundary handling [#332] (Axel Donath) 
 
 .. _gammapy_0p3_release:
 

docs/development/howto.rst

@@ -704,7 +704,7 @@ Note that if you distribute Gammapy together with one of the GPL dependencies,
 the whole distribution then falls under the GPL license.
 
 Gammapy plotting style
-++++++++++++++++++++++
+----------------------
 
 Figures and plots in the Gammapy docs use the same consistent plotting style,
 that is defined in `gammapy.utils.mpl_style`.  The style is derived from the

gammapy/detect/_test_statistics_cython.pyx

@@ -34,8 +34,9 @@ def _f_cash_root_cython(np.float_t x, np.ndarray[np.float_t, ndim=2] counts,
     sum = 0
     for j in range(nj):
         for i in range(ni):
-            sum += (model[j, i] * (counts[j, i] / (x * FLUX_FACTOR * model[j, i]
-                                                   + background[j, i]) - 1))
+            if model[j, i] > 0:
+                sum += (model[j, i] * (counts[j, i] / (x * FLUX_FACTOR * model[j, i]
+                                                       + background[j, i]) - 1))
     return sum
 
 
@@ -57,8 +58,8 @@ def _amplitude_bounds_cython(np.ndarray[np.float_t, ndim=2] counts,
         Source template (multiplied with exposure).
     """
 
-    cdef float s_model = 0, s_counts = 0, sn, sn_min = 1E14, c_min
-    cdef float b_min, b_max
+    cdef np.float_t s_model = 0, s_counts = 0, sn, sn_min = 1E14, c_min = 1
+    cdef np.float_t b_min, b_max
     cdef unsigned int i, j, ni, nj
     ni = counts.shape[1]
     nj = counts.shape[0]
@@ -67,7 +68,7 @@ def _amplitude_bounds_cython(np.ndarray[np.float_t, ndim=2] counts,
             s_model += model[j, i]
             if counts[j, i] > 0:
                 s_counts += counts[j, i]
-                if model[j, i] != 0:
+                if model[j, i] > 0:
                     sn = background[j, i] / model[j, i]
                     if sn < sn_min:
                         sn_min = sn

gammapy/detect/test_statistics.py

@@ -166,7 +166,7 @@ def compute_ts_map_multiscale(maps, psf_parameters, scales=[0], downsample='auto
         if downsample == 'auto':
             factor = int(np.select([scale < 5 * BINSZ, scale < 10 * BINSZ,
                                     scale < 20 * BINSZ, scale < 40 * BINSZ],
-                                   [1, 2, 4, 8], 16))
+                                   [1, 2, 4, 4], 8))
         else:
             factor = int(downsample)
         if factor == 1:
@@ -323,13 +323,22 @@ def compute_ts_map(counts, background, exposure, kernel, mask=None, flux=None,
     ----------
     [Stewart2009]_
     """
-    from scipy.ndimage.morphology import binary_erosion
     from time import time
     t_0 = time()
 
     assert counts.shape == background.shape
     assert counts.shape == exposure.shape
 
+    # in some maps there are pixels, which have exposure, but zero
+    # background, which doesn't make sense and causes the TS computation
+    # to fail, this is a temporary fix
+    mask_ = np.logical_and(background == 0,  exposure > 0)
+    if mask_.any():
+        log.warn('There are pixels in the data, that have exposure, but zero '
+                 'background, which can cause the ts computation to fail. '
+                 'Setting exposure of this pixels to zero.')
+        exposure[mask_] = 0
+
     if (flux is None and method != 'root brentq') or threshold is not None:
         from scipy.ndimage import convolve
         radius = _flux_correlation_radius(kernel)
@@ -349,7 +358,7 @@ def compute_ts_map(counts, background, exposure, kernel, mask=None, flux=None,
 
     # Positions where exposure == 0 are not processed
     if mask is None:
-        mask = binary_erosion(exposure > 0, np.ones(np.array(kernel.shape) + 2))
+        mask = exposure > 0
     positions = [(j, i) for j, i in positions if mask[j][i]]
 
     wrap = partial(_ts_value, counts=counts, exposure=exposure,

gammapy/scripts/ts_image.py

@@ -5,7 +5,7 @@
 import json
 import logging
 log = logging.getLogger(__name__)
-from ..utils.scripts import get_parser
+from ..utils.scripts import get_parser, set_up_logging_from_args
 
 __all__ = ['ts_image']
 
@@ -42,7 +42,11 @@ def main(args=None):
                         " that the fit is done at all.")
     parser.add_argument('--overwrite', action='store_true',
                         help='Overwrite output files.')
+    parser.add_argument("-l", "--loglevel", default='info',
+                        choices=['debug', 'info', 'warning', 'error', 'critical'],
+                        help="Set the logging level")
     args = parser.parse_args()
+    set_up_logging_from_args(args)
     ts_image(**vars(args))