(Re)added functionality to process datacubes. Includes testing. Rebased.

photutils/aperture/core.py

@@ -393,21 +393,22 @@ def do_photometry(self, data, error=None, mask=None, method='exact',
 
         aperture_sums = []
         aperture_sum_errs = []
+        output_shape = (1,) if data.ndim==2 else data.shape[0]
         for mask in self.to_mask(method=method, subpixels=subpixels):
             data_cutout = mask.cutout(data)
 
             if data_cutout is None:
-                aperture_sums.append(np.nan)
+                aperture_sums.append(np.repeat(np.nan, output_shape))
             else:
-                aperture_sums.append(np.sum(data_cutout * mask.data))
+                aperture_sums.append(np.sum(data_cutout * mask.data, axis=(-2,-1)))
 
             if error is not None:
                 error_cutout = mask.cutout(error)
 
                 if error_cutout is None:
-                    aperture_sum_errs.append(np.nan)
+                    aperture_sum_errs.append(np.repeat(np.nan, output_shape))
                 else:
-                    aperture_var = np.sum(error_cutout ** 2 * mask.data)
+                    aperture_var = np.sum(error_cutout ** 2 * mask.data, axis=(-2,-1))
                     aperture_sum_errs.append(np.sqrt(aperture_var))
 
         # handle Quantity objects and input units
@@ -726,8 +727,8 @@ def _prepare_photometry_input(data, error, mask, wcs, unit):
             pass
 
     data = np.asanyarray(data)
-    if data.ndim != 2:
-        raise ValueError('data must be a 2D array.')
+    if (data.ndim != 2) and (data.ndim !=3):
+        raise ValueError('data must be a 2D or 3D array.')
 
     if unit is not None:
         unit = u.Unit(unit, parse_strict='warn')

photutils/aperture/mask.py

@@ -178,8 +178,8 @@ def cutout(self, data, fill_value=0., copy=False):
         """
 
         data = np.asanyarray(data)
-        if data.ndim != 2:
-            raise ValueError('data must be a 2D array.')
+        if (data.ndim != 2) and (data.ndim != 3):
+            raise ValueError('data must be a 2D or 3D array.')
 
         partial_overlap = False
         if self.bbox.ixmin < 0 or self.bbox.iymin < 0:
@@ -189,20 +189,21 @@ def cutout(self, data, fill_value=0., copy=False):
             # try this for speed -- the result may still be a partial
             # overlap, in which case the next block will be triggered
             if copy:
-                cutout = np.copy(data[self.bbox.slices])
+                cutout = np.copy(data[(Ellipsis,)+self.bbox.slices])
             else:
-                cutout = data[self.bbox.slices]
+                cutout = data[(Ellipsis,)+self.bbox.slices]
 
-        if partial_overlap or (cutout.shape != self.shape):
-            slices_large, slices_small = self._overlap_slices(data.shape)
+        if partial_overlap or (cutout.shape[-2:] != self.shape):
+            slices_large, slices_small = self._overlap_slices(data.shape[-2:])
 
             if slices_small is None:
                 return None    # no overlap
 
             # cutout is a copy
-            cutout = np.zeros(self.shape, dtype=data.dtype)
+            output_shape = self.shape if data.ndim==2 else (data.shape[0],)+self.shape
+            cutout = np.zeros(output_shape, dtype=data.dtype)
             cutout[:] = fill_value
-            cutout[slices_small] = data[slices_large]
+            cutout[(Ellipsis,)+slices_small] = data[(Ellipsis,)+slices_large]
 
             if isinstance(data, u.Quantity):
                 cutout = u.Quantity(cutout, unit=data.unit)