Repository: gammapy/gammapy
Pull request #5280 (2024-05-31 12:43:24)

Add min_pix option in WcsGeom.cutout and set it to 3 for IRFs 
Add min_pix option in WcsGeom.cutout and set it to 3 for IRFs in order to ensure that IRF interpolation is done correctly if the binning of the IRFs is larger than the width of the analysis region.

bkhelifi: Thanks @QRemy... What happens where there is not enough bin for the interpolation?
QRemy: > What happens where there is not enough bin for the interpolation?
It uses the single bin that is left in the psf cutout. So if we reduce the region size the psf used at a given position can be different to the one used before cutout.


Patch:
diff --git a/gammapy/irf/core.py b/gammapy/irf/core.py
index a0fdabb77..41e8f3168 100644
--- a/gammapy/irf/core.py
+++ b/gammapy/irf/core.py
@@ -939,7 +939,7 @@ class IRFMap:
         """Copy IRF map."""
         return deepcopy(self)
 
-    def cutout(self, position, width, mode="trim"):
+    def cutout(self, position, width, mode="trim", min_npix=3):
         """Cutout IRF map.
 
         Parameters
@@ -952,15 +952,22 @@ class IRFMap:
         mode : {'trim', 'partial', 'strict'}, optional
             Mode option for Cutout2D, for details see `~astropy.nddata.utils.Cutout2D`.
             Default is "trim".
+        min_npix : bool, optional
+            Force width to a minimmum number of pixels.
+            Default is 3. The default is 3 pixels so interpolation is done correctly
+            if the binning of the IRF is larger than the width of the analysis region.
 
         Returns
         -------
         cutout : `IRFMap`
             Cutout IRF map.
         """
-        irf_map = self._irf_map.cutout(position, width, mode)
+
+        irf_map = self._irf_map.cutout(position, width, mode, min_npix=min_npix)
         if self.exposure_map:
-            exposure_map = self.exposure_map.cutout(position, width, mode)
+            exposure_map = self.exposure_map.cutout(
+                position, width, mode, min_npix=min_npix
+            )
         else:
             exposure_map = None
         return self.__class__(irf_map, exposure_map=exposure_map)
diff --git a/gammapy/maps/wcs/geom.py b/gammapy/maps/wcs/geom.py
index f5b21796f..ee9d3ff50 100644
--- a/gammapy/maps/wcs/geom.py
+++ b/gammapy/maps/wcs/geom.py
@@ -885,7 +885,7 @@ class WcsGeom(Geom):
         coord = self.to_image().get_coord()
         return center.separation(coord.skycoord)
 
-    def cutout(self, position, width, mode="trim", odd_npix=False):
+    def cutout(self, position, width, mode="trim", odd_npix=False, min_npix=1):
         """
         Create a cutout around a given position.
 
@@ -902,6 +902,10 @@ class WcsGeom(Geom):
         odd_npix : bool, optional
             Force width to odd number of pixels.
             Default is False.
+        min_npix : bool, optional
+            Force width to a minimmum number of pixels.
+            Default is 1.
+
 
         Returns
         -------
@@ -911,7 +915,7 @@ class WcsGeom(Geom):
         width = _check_width(width) * u.deg
 
         binsz = self.pixel_scales
-        width_npix = np.clip((width / binsz).to_value(""), 1, None)
+        width_npix = np.clip((width / binsz).to_value(""), min_npix, None)
 
         if odd_npix:
             width_npix = round_up_to_odd(width_npix)
diff --git a/gammapy/maps/wcs/ndmap.py b/gammapy/maps/wcs/ndmap.py
index cf57345bb..97f34e237 100644
--- a/gammapy/maps/wcs/ndmap.py
+++ b/gammapy/maps/wcs/ndmap.py
@@ -975,7 +975,7 @@ class WcsNDMap(WcsMap):
 
         return self._init_copy(data=smoothed_data)
 
-    def cutout(self, position, width, mode="trim", odd_npix=False):
+    def cutout(self, position, width, mode="trim", odd_npix=False, min_npix=1):
         """
         Create a cutout around a given position.
 
@@ -992,6 +992,9 @@ class WcsNDMap(WcsMap):
         odd_npix : bool, optional
             Force width to odd number of pixels.
             Default is False.
+        min_npix : bool, optional
+            Force width to a minimmum number of pixels.
+            Default is 1.
 
         Returns
         -------
@@ -999,7 +1002,11 @@ class WcsNDMap(WcsMap):
             Cutout map.
         """
         geom_cutout = self.geom.cutout(
-            position=position, width=width, mode=mode, odd_npix=odd_npix
+            position=position,
+            width=width,
+            mode=mode,
+            odd_npix=odd_npix,
+            min_npix=min_npix,
         )
         cutout_info = geom_cutout.cutout_slices(self.geom, mode=mode)
 


Test patch:
diff --git a/gammapy/irf/psf/tests/test_map.py b/gammapy/irf/psf/tests/test_map.py
index 786221d1a..d2ef82787 100644
--- a/gammapy/irf/psf/tests/test_map.py
+++ b/gammapy/irf/psf/tests/test_map.py
@@ -129,6 +129,14 @@ def test_psf_map_containment():
     assert_allclose(psf_map.containment(rad=10 * u.deg, energy_true=[10] * u.TeV), 1)
 
 
+def test_psf_map_cutout():
+    psf_map = make_test_psfmap(0.15 * u.deg)
+    psf_map_cutout = psf_map.cutout(
+        position=SkyCoord(1, 1, unit="deg"), width=(1, 0.01)
+    )
+    assert_allclose(psf_map_cutout.psf_map.geom.data_shape, (4, 100, 3, 5))
+
+
 def test_psfmap_to_psf_kernel():
     psfmap = make_test_psfmap(0.15 * u.deg)
 
diff --git a/gammapy/maps/wcs/tests/test_geom.py b/gammapy/maps/wcs/tests/test_geom.py
index 368603afd..7bc7cf6a4 100644
--- a/gammapy/maps/wcs/tests/test_geom.py
+++ b/gammapy/maps/wcs/tests/test_geom.py
@@ -313,9 +313,19 @@ def test_cutout_min_size():
     geom = WcsGeom.create(skydir=(0, 0), npix=10, binsz=0.5)
     position = SkyCoord(0, 0, unit="deg")
     cutout_geom = geom.cutout(position=position, width=["2 deg", "0.1 deg"])
-
     assert cutout_geom.data_shape == (1, 4)
 
+    cutout_geom = geom.cutout(position=position, width=["2 deg", "0.1 deg"], min_npix=3)
+    assert cutout_geom.data_shape == (3, 4)
+
+    geom = WcsGeom.create(skydir=(0, 0), npix=(10, 1), binsz=0.5)
+    position = SkyCoord(0, 0, unit="deg")
+    cutout_geom = geom.cutout(position=position, width=["1 deg", "0.1 deg"])
+    assert cutout_geom.data_shape == (1, 2)
+
+    cutout_geom = geom.cutout(position=position, width=["1 deg", "0.1 deg"], min_npix=3)
+    assert cutout_geom.data_shape == (1, 3)
+
 
 def test_wcs_geom_get_coord():
     geom = WcsGeom.create(
diff --git a/gammapy/maps/wcs/tests/test_ndmap.py b/gammapy/maps/wcs/tests/test_ndmap.py
index 5ec433b59..022ba4a99 100644
--- a/gammapy/maps/wcs/tests/test_ndmap.py
+++ b/gammapy/maps/wcs/tests/test_ndmap.py
@@ -581,6 +581,10 @@ def test_make_cutout(mode):
     assert_allclose(actual, 36.0)
     assert_allclose(cutout.geom.shape_axes, m.geom.shape_axes)
     assert_allclose(cutout.geom.width.to_value("deg"), [[2.0], [3.0]])
+    assert_allclose(cutout.geom.data_shape, (3, 2, 3, 2))
+
+    cutout = m.cutout(position=pos, width=(2.0, 3.0) * u.deg, mode=mode, min_npix=3)
+    assert_allclose(cutout.geom.data_shape, (3, 2, 3, 3))
 
 
 def test_convolve_vs_smooth():