Repository: gammapy/gammapy
Pull request #5392 (2024-07-25 14:19:46)

Backport PR #5280: 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 a691e9d8d..c4704ff2a 100644
--- a/gammapy/irf/core.py
+++ b/gammapy/irf/core.py
@@ -907,7 +907,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
@@ -919,15 +919,23 @@ class IRFMap:
             If only one value is passed, a square region is extracted.
         mode : {'trim', 'partial', 'strict'}
             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 3b578f84d..7867f761b 100644
--- a/gammapy/maps/wcs/geom.py
+++ b/gammapy/maps/wcs/geom.py
@@ -858,7 +858,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.
 
@@ -873,6 +873,11 @@ class WcsGeom(Geom):
             Mode option for Cutout2D, for details see `~astropy.nddata.utils.Cutout2D`.
         odd_npix : bool
             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
         -------
@@ -882,7 +887,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 c59f1a758..ae327cf01 100644
--- a/gammapy/maps/wcs/ndmap.py
+++ b/gammapy/maps/wcs/ndmap.py
@@ -820,7 +820,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.
 
@@ -835,6 +835,10 @@ class WcsNDMap(WcsMap):
             Mode option for Cutout2D, for details see `~astropy.nddata.utils.Cutout2D`.
         odd_npix : bool
             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
         -------
@@ -842,7 +846,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 580206275..60ac23bb3 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 8a62e02ae..d63e18d45 100644
--- a/gammapy/maps/wcs/tests/test_geom.py
+++ b/gammapy/maps/wcs/tests/test_geom.py
@@ -300,9 +300,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 f98e4a14f..a8f69d77e 100644
--- a/gammapy/maps/wcs/tests/test_ndmap.py
+++ b/gammapy/maps/wcs/tests/test_ndmap.py
@@ -565,6 +565,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():