Add .select_region() and .select_mask() methods to Models

gammapy/datasets/map.py

@@ -12,7 +12,12 @@
 from gammapy.data import GTI
 from gammapy.irf import EDispKernelMap, EDispMap, PSFKernel, PSFMap
 from gammapy.maps import Map, MapAxis, RegionGeom, WcsGeom
-from gammapy.modeling.models import BackgroundModel, DatasetModels, FoVBackgroundModel
+from gammapy.modeling.models import (
+    Models,
+    BackgroundModel,
+    DatasetModels,
+    FoVBackgroundModel,
+)
 from gammapy.stats import (
     CashCountsStatistic,
     WStatCountsStatistic,
@@ -47,6 +52,16 @@
 USE_NPRED_CACHE = True
 
 
+def get_cutout_width(model, psf=None, margin=CUTOUT_MARGIN):
+    """Cutout width for the model component"""
+    if psf is not None:
+        psf_width = np.max(psf.psf_kernel_map.geom.width)
+    else:
+        psf_width = 0 * u.deg
+
+    return psf_width + 2 * (model.evaluation_radius + margin)
+
+
 def create_map_dataset_geoms(
     geom, energy_axis_true=None, migra_axis=None, rad_axis=None, binsz_irf=None,
 ):
@@ -404,7 +419,12 @@ def npred_signal(self, model=None):
 
             if evaluator.needs_update:
                 evaluator.update(
-                    self.exposure, self.psf, self.edisp, self._geom, self.mask_safe_psf
+                    self.exposure,
+                    self.psf,
+                    self.edisp,
+                    self._geom,
+                    self.mask_fit,
+                    self.mask_safe_psf,
                 )
 
             if evaluator.contributes:
@@ -2399,6 +2419,8 @@ class MapEvaluator:
         PSF kernel
     edisp : `~gammapy.irf.EDispKernel`
         Energy dispersion
+    mask : `~gammapy.maps.Map`
+        Mask to apply to the likelihood for fitting.
     gti : `~gammapy.data.GTI`
         GTI of the observation or union of GTI if it is a stacked observation
     evaluation_mode : {"local", "global"}
@@ -2418,6 +2440,7 @@ def __init__(
         psf=None,
         edisp=None,
         gti=None,
+        mask=None,
         evaluation_mode="local",
         use_cache=True,
     ):
@@ -2426,6 +2449,7 @@ def __init__(
         self.exposure = exposure
         self.psf = psf
         self.edisp = edisp
+        self.mask = mask
         self.gti = gti
         self.contributes = True
         self.use_cache = use_cache
@@ -2492,14 +2516,9 @@ def needs_update(self):
     @property
     def cutout_width(self):
         """Cutout width for the model component"""
-        if self.psf is not None:
-            psf_width = np.max(self.psf.psf_kernel_map.geom.width)
-        else:
-            psf_width = 0 * u.deg
-
-        return psf_width + 2 * (self.model.evaluation_radius + CUTOUT_MARGIN)
+        return get_cutout_width(self.model, psf=self.psf)
 
-    def update(self, exposure, psf, edisp, geom, mask_safe_psf):
+    def update(self, exposure, psf, edisp, geom, mask_fit, mask_safe_psf):
         """Update MapEvaluator, based on the current position of the model component.
 
         Parameters
@@ -2512,6 +2531,8 @@ def update(self, exposure, psf, edisp, geom, mask_safe_psf):
             Edisp map.
         geom : `WcsGeom`
             Counts geom
+        mask_fit : `~gammapy.maps.Map`
+            Mask to apply to the likelihood for fitting.
         mask_safe_psf : `~gammapy.maps.Map`
             Mask safe map of boolean type.
         """
@@ -2543,6 +2564,9 @@ def update(self, exposure, psf, edisp, geom, mask_safe_psf):
                 self.irf_position, energy_axis=energy_axis
             )
 
+        if mask_fit is None:
+            mask_fit = Map.from_geom(geom, data=np.ones(geom.data_shape, dtype=bool))
+
         # lookup psf
         if psf:
             if self.apply_psf_after_edisp:
@@ -2557,13 +2581,15 @@ def update(self, exposure, psf, edisp, geom, mask_safe_psf):
 
         if self.evaluation_mode == "local" and self.model.evaluation_radius is not None:
             self._init_position = self.model.position
+            self.contributes = self.model.contributes(
+                mask_fit, margin=self.cutout_width, use_evaluation_region=True
+            )
             try:
                 self.exposure = exposure.cutout(
                     position=self.model.position, width=self.cutout_width
                 )
-                self.contributes = True
             except (NoOverlapError, ValueError):
-                self.contributes = False
+                pass
         else:
             self.exposure = exposure
 

gammapy/datasets/simulate.py

@@ -73,6 +73,7 @@ def sample_sources(self, dataset):
                     dataset.psf,
                     dataset.edisp,
                     dataset._geom,
+                    dataset.mask_fit,
                     dataset.mask_safe_psf,
                 )
 

gammapy/estimators/ts_map.py

@@ -209,6 +209,7 @@ def estimate_kernel(self, dataset):
             dataset.psf,
             dataset.edisp,
             dataset.counts.geom,
+            dataset.mask_fit,
             dataset.mask_safe_psf,
         )
 

gammapy/maps/core.py

@@ -390,7 +390,7 @@ def pad(self, pad_width, mode="constant", cval=0, order=1):
         cval : float
             Padding value when mode='consant'.
         order : int
-            Order of interpolation when mode='constant' (0 =
+            Order of interpolation when mode='interp' (0 =
             nearest-neighbor, 1 = linear, 2 = quadratic, 3 = cubic).
 
         Returns

gammapy/maps/tests/test_wcsnd.py

@@ -295,7 +295,7 @@ def test_wcsndmap_interp_by_coord(npix, binsz, frame, proj, skydir, axes):
     assert_allclose(coords[1].value, m.interp_by_coord(coords, method="linear"))
     assert_allclose(coords[1].value, m.interp_by_coord(coords, method="linear"))
 
-    #if geom.is_regular and not geom.is_allsky:
+    # if geom.is_regular and not geom.is_allsky:
     #    assert_allclose(
     #        coords[1].to_value("deg"), m.interp_by_coord(coords, interp="cubic")
     #    )
@@ -817,5 +817,7 @@ def test_mask_fit_modifications():
     mask_fit = mask_fit & geom.boundary_mask(width=(0.3 * u.deg, 0.1 * u.deg))
     assert np.sum(mask_fit.data[0, :, :]) == 6300
     assert np.sum(mask_fit.data[1, :, :]) == 6300
-    mask_fit = mask_fit.binary_dilate(width=(0.3 * u.deg, 0.1 * u.deg))
+    mask_fit_fft = mask_fit.binary_dilate(width=(0.3 * u.deg, 0.1 * u.deg))
+    mask_fit = mask_fit.binary_dilate(width=(0.3 * u.deg, 0.1 * u.deg), use_fft=False)
+    assert np.sum(mask_fit_fft.data) == np.prod(mask_fit_fft.data.shape)
     assert np.sum(mask_fit.data) == np.prod(mask_fit.data.shape)

gammapy/maps/wcs.py

@@ -877,8 +877,8 @@ def region_mask(self, regions, inside=True):
 
         Parameters
         ----------
-        regions : list of  `~regions.Region`
-            Python list of regions (pixel or sky regions accepted)
+        regions : list of `~regions.Region`
+            region or list of regions (pixel or sky regions accepted)
         inside : bool
             For ``inside=True``, pixels in the region to True (the default).
             For ``inside=False``, pixels in the region are False.
@@ -915,6 +915,9 @@ def region_mask(self, regions, inside=True):
         if not self.is_regular:
             raise ValueError("Multi-resolution maps not supported yet")
 
+        if not isinstance(regions, list):
+            regions = [regions]
+
         idx = self.get_idx()
         pixcoord = PixCoord(idx[0], idx[1])
 

gammapy/maps/wcsnd.py

@@ -10,7 +10,7 @@
 from astropy.convolution import Tophat2DKernel
 from astropy.coordinates import SkyCoord
 from astropy.io import fits
-from regions import PointSkyRegion, RectangleSkyRegion
+from regions import PointSkyRegion, RectangleSkyRegion, SkyRegion, PixCoord
 from gammapy.extern.skimage import block_reduce
 from gammapy.utils.interpolation import ScaledRegularGridInterpolator
 from gammapy.utils.random import InverseCDFSampler, get_random_state
@@ -25,7 +25,7 @@
 log = logging.getLogger(__name__)
 
 
-def _apply_binary_operations(map_, width, func):
+def _apply_binary_operations(map_, width, func, output=None, mode="same"):
     """ Apply ndi.binary_dilation or ndi.binary_erosion to a boolean-mask map"""
 
     if not map_.is_mask:
@@ -34,17 +34,23 @@ def _apply_binary_operations(map_, width, func):
     if not isinstance(width, tuple):
         width = (width, width)
 
+    if output is None:
+        output = map_
+
     shape = tuple(
         [
             int(np.ceil(x / scale).value * 2 + 1)
             for x, scale in zip(width, map_.geom.pixel_scales)
         ]
     )
 
-    mask_data = np.empty(map_.data.shape, dtype=bool)
-
+    mask_data = np.empty(output.data.shape, dtype=bool)
+    structure = np.ones(shape)
     for img, idx in map_.iter_by_image():
-        mask_data[idx] = func(img, structure=np.ones(shape))
+        if func == scipy.signal.fftconvolve:
+            mask_data[idx] = func(img, structure, mode=mode)
+        else:
+            mask_data[idx] = func(img, structure)
 
     return mask_data
 
@@ -176,13 +182,16 @@ def _interp_by_coord_griddata(self, coords, method="linear"):
 
         data = self.data[np.isfinite(self.data)]
         vals = scipy.interpolate.griddata(
-           tuple(grid_coords.flat), data, tuple(coords), method=method
+            tuple(grid_coords.flat), data, tuple(coords), method=method
         )
 
         m = ~np.isfinite(vals)
         if np.any(m):
             vals_fill = scipy.interpolate.griddata(
-                tuple(grid_coords.flat), data, tuple([c[m] for c in coords]), method="nearest"
+                tuple(grid_coords.flat),
+                data,
+                tuple([c[m] for c in coords]),
+                method="nearest",
             )
             vals[m] = vals_fill
 
@@ -569,6 +578,34 @@ def get_spectrum(self, region=None, func=np.nansum, weights=None):
 
         return self.to_region_nd_map(region=region, func=func, weights=weights)
 
+    def mask_contains_region(self, regions):
+        """Check if input regions are overlaping with a boolean mask map.
+
+        Parameters
+        ----------
+        regions: `~regions.Region`
+             Region or list of Regions (pixel or sky regions accepted).
+
+        Returns
+        -------
+        overlap : `numpy.array`
+            Boolean array of same length than regions
+        """
+
+        if not self.is_mask:
+            raise ValueError("mask_contains_region is only supported for boolean masks")
+
+        pixcoords = self.geom.get_idx()
+        pixcoords = PixCoord(pixcoords[0][self.data], pixcoords[1][self.data])
+        if not isinstance(regions, list):
+            regions = [regions]
+        overlap = np.zeros(len(regions), dtype=bool)
+        for k, region in enumerate(regions):
+            if isinstance(region, SkyRegion):
+                region = region.to_pixel(self.geom.wcs)
+            overlap[k] = np.any(region.contains(pixcoords))
+        return overlap
+
     def binary_erode(self, width):
         """Binary erosion of boolean mask removing a given margin
 
@@ -584,26 +621,49 @@ def binary_erode(self, width):
             Eroded mask map
 
         """
+
         mask_data = _apply_binary_operations(self, width, ndi.binary_erosion)
         return self._init_copy(data=mask_data)
 
-    def binary_dilate(self, width):
+    def binary_dilate(self, width, use_fft=True, mode="same"):
         """Binary dilation of boolean mask addding a given margin
 
         Parameters
         ----------
         width : tuple of `~astropy.units.Quantity`
             Angular sizes of the margin in (lon, lat) in that specific order.
             If only one value is passed, the same margin is applied in (lon, lat).
+        use_fft : bool
+            Use `scipy.signal.fftconvolve` if True (default)
+            and `ndi.binary_dilation` otherwise.
+
+        mode : str {'same', 'full'}
+            A string indicating the size of the output.
+            - 'same': The output is the same size as input (Default).
+            - 'full': The output size is extended by the width
 
         Returns
         -------
         map : `WcsNDMap`
             Dilated mask map
-
         """
-        mask_data = _apply_binary_operations(self, width, ndi.binary_dilation)
-        return self._init_copy(data=mask_data)
+
+        if use_fft:
+            func = scipy.signal.fftconvolve
+        else:
+            func = ndi.binary_dilation
+
+        if mode == "full":
+            pad_width = u.Quantity(width) / (self.geom.pixel_scales)
+            pad_width = list(np.ceil(pad_width.value).astype(int))
+            mask = self.pad(pad_width)
+        else:
+            mask = self
+
+        mask_data = _apply_binary_operations(self, width, func, output=mask, mode=mode)
+        if use_fft:
+            mask_data = np.rint(mask_data).astype(bool)
+        return self._init_copy(geom=mask.geom, data=mask_data)
 
     def convolve(self, kernel, use_fft=True, **kwargs):
         """
@@ -618,7 +678,8 @@ def convolve(self, kernel, use_fft=True, **kwargs):
         kernel : `~gammapy.irf.PSFKernel` or `numpy.ndarray`
             Convolution kernel.
         use_fft : bool
-            Use `scipy.signal.fftconvolve` or `ndi.convolve`.
+            Use `scipy.signal.fftconvolve` if True (default)
+            and `ndi.convolve` otherwise.
         kwargs : dict
             Keyword arguments passed to `scipy.signal.fftconvolve` or
             `ndi.convolve`.