Repository: gammapy/gammapy
Pull request #5363 (2024-07-05 13:30:17)

Backport PR #5280 on branch v1.2.x (Add min_pix option in WcsGeom.cutout and set it to 3 for IRFs )
Backport PR #5280: Add min_pix option in WcsGeom.cutout and set it to 3 for IRFs 


Patch:
diff --git a/gammapy/irf/core.py b/gammapy/irf/core.py
index d2622df2b..42dda97d5 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():