Add WCS map cutout method

gammapy/maps/tests/test_wcsnd.py

@@ -369,22 +369,25 @@ def test_wcsndmap_upsample(npix, binsz, coordsys, proj, skydir, axes):
     m2 = m.upsample(2, order=0, preserve_counts=True)
     assert_allclose(np.nansum(m.data), np.nansum(m2.data))
 
+
 def test_coadd_unit():
-    geom = WcsGeom.create(npix=(10,10), binsz=1,
+    geom = WcsGeom.create(npix=(10, 10), binsz=1,
                           proj='CAR', coordsys='GAL')
-    m1 = WcsNDMap(geom, data=np.ones((10,10)), unit='m2')
-    m2 = WcsNDMap(geom, data=np.ones((10,10)), unit='cm2')
+    m1 = WcsNDMap(geom, data=np.ones((10, 10)), unit='m2')
+    m2 = WcsNDMap(geom, data=np.ones((10, 10)), unit='cm2')
 
     m1.coadd(m2)
 
     assert_allclose(m1.data, 1.0001)
 
+
 def test_make_region_mask():
     from regions import CircleSkyRegion
-    geom = WcsGeom.create(npix=(3,3), binsz=2,
+    geom = WcsGeom.create(npix=(3, 3), binsz=2,
                           proj='CAR', coordsys='GAL')
     m = WcsNDMap(geom)
-    region = CircleSkyRegion(SkyCoord(0, 0, unit='deg', frame='galactic'), 1.0*u.deg)
+    region = CircleSkyRegion(SkyCoord(0, 0, unit='deg', frame='galactic'),
+                             1.0 * u.deg)
     maskmap = m.make_region_mask(region)
 
     assert maskmap.data.dtype == bool
@@ -405,3 +408,14 @@ def test_smooth(kernel):
     smoothed = m.smooth(0.2 * u.deg, kernel)
     actual = smoothed.data.sum()
     assert_allclose(actual, desired)
+
+
+def test_make_cutout():
+    pos = SkyCoord(0, 0, unit='deg', frame='galactic')
+    geom = WcsGeom.create(npix=(10, 10), binsz=1, skydir=pos,
+                          proj='CAR', coordsys='GAL', axes=axes2)
+    m = WcsNDMap(geom, data=np.ones((3, 2, 10, 10)), unit='m2')
+    cutout = m.make_cutout(position=pos, width=(2.0, 3.0) * u.deg)
+    actual = cutout.data.sum()
+    assert_allclose(actual, 36.0)
+    assert_allclose(cutout.geom.shape, m.geom.shape)

gammapy/maps/wcsnd.py

@@ -5,6 +5,8 @@
 import numpy as np
 from astropy.io import fits
 from astropy.units import Quantity
+from astropy.nddata import Cutout2D
+import astropy.units as u
 from astropy.convolution import Tophat2DKernel
 from .utils import unpack_seq
 from .geom import pix_tuple_to_idx, axes_pix_to_coord
@@ -599,3 +601,43 @@ def smooth(self, radius, kernel='gauss', **kwargs):
         image = copy.copy(self)
         image.data = data
         return image
+
+    def make_cutout(self, position, width, mode="strict", copy=True):
+        """
+        Create a cutout of a WcsNDMap around a given position.
+
+        Parameters
+        ----------
+        position : `~astropy.coordinates.SkyCoord`
+            Center position of the cutout region.
+        width : tuple of `~astropy.coordinates.Angle`
+            Angular sizes of the region in (lon, lat). If only one value is passed,
+            a square region is extracted. For more options see also `~astropy.nddata.utils.Cutout2D`.
+        mode : {'trim', 'partial', 'strict'}
+            Mode option for Cutout2D, for details see `~astropy.nddata.utils.Cutout2D`.
+
+        copy : bool, optional
+               If False (default), then the cutout data will be a view into the original data  array. 
+               If True, then the cutout data will hold a copy of the original data array.
+
+        Returns
+        -------
+        cutout : `~gammapy.maps.WcsNDMap`
+            The cutout map itself
+        """
+
+        idx = (0,) * len(self.geom.axes)
+
+        cutout2d = Cutout2D(data=self.data[idx], wcs=self.geom.wcs,
+                            position=position, size=width, mode=mode)
+
+        # Create the slices with the non-spatial axis
+        cutout_slices = Ellipsis, cutout2d.slices_original[0], cutout2d.slices_original[1]
+
+        geom = WcsGeom(cutout2d.wcs, cutout2d.shape[::-1], axes=self.geom.axes)
+        data = self.data[cutout_slices]
+
+        if copy:
+            data=data.copy()
+
+        return WcsNDMap(geom, data, meta=self.meta, unit=self.unit)