Merge pull request #375 from keflavich/mad_std

add mad_std and fix up the docs for regular std

spectral_cube/spectral_cube.py

@@ -20,6 +20,7 @@
 from astropy import log
 from astropy import wcs
 from astropy import convolution
+from astropy import stats
 
 import numpy as np
 
@@ -290,7 +291,9 @@ def apply_numpy_function(self, function, fill=np.nan,
             spatial axes?
         unit : None or `astropy.units.Unit`
             The unit to include for the output array.  For example,
-            `SpectralCube.max` calls ``SpectralCube.apply_numpy_function(np.max, unit=self.unit)``, inheriting the unit from the original cube.
+            `SpectralCube.max` calls
+            ``SpectralCube.apply_numpy_function(np.max, unit=self.unit)``,
+            inheriting the unit from the original cube.
             However, for other numpy functions, e.g. `numpy.argmax`, the return
             is an index and therefore unitless.
         check_endian : bool
@@ -495,6 +498,19 @@ def _count_nonzero_slicewise(self, axis=None):
     def std(self, axis=None, how='cube', ddof=0):
         """
         Return the standard deviation of the cube, optionally over an axis.
+
+        Parameters
+        ----------
+        axis : int or None
+            The axis to average over
+        how : str
+            The method to use.  This method supports slice-wise standard
+            deviation calculation ('slice') and full-cube ('cube'), but not
+            'ray'.
+        ddof : int
+            Means Delta Degrees of Freedom.  The divisor used in calculations
+            is ``N - ddof``, where ``N`` represents the number of elements.  By
+            default `ddof` is zero.
         """
 
         projection = self._naxes_dropped(axis) in (1,2)
@@ -538,6 +554,19 @@ def std(self, axis=None, how='cube', ddof=0):
                                          axis=axis, unit=self.unit,
                                          projection=projection)
 
+    @aggregation_docstring
+    def mad_std(self, axis=None):
+        """
+        Use astropy's mad_std to computer the standard deviation
+        """
+        if int(astropy.__version__[0]) < 2:
+            raise NotImplementedError("mad_std requires astropy >= 2")
+        projection = self._naxes_dropped(axis) in (1,2)
+        return self.apply_numpy_function(stats.mad_std, fill=np.nan,
+                                         how='cube', axis=axis, unit=self.unit,
+                                         ignore_nan=True,
+                                         projection=projection)
+
 
     @aggregation_docstring
     def max(self, axis=None, how='auto'):

spectral_cube/tests/test_spectral_cube.py

@@ -8,6 +8,7 @@
 
 import pytest
 
+import astropy
 from astropy.io import fits
 from astropy import units as u
 from astropy.wcs import WCS
@@ -1322,3 +1323,26 @@ def test_mask_bad_beams():
 
     mean2 = masked_cube2.mean(axis=0)
     assert np.all(mean2 == (cube[2,:,:]+cube[1,:,:])/2)
+
+def test_mad_std():
+    cube, data = cube_and_raw('adv.fits')
+
+    if int(astropy.__version__[0]) < 2:
+        with pytest.raises(NotImplementedError) as exc:
+            cube.mad_std()
+
+    else:
+        # mad_std run manually on data
+        result = np.array([[0.15509701,  0.45763670],
+                           [0.55907956,  0.42932451],
+                           [0.48819454,  0.25499305]])
+
+        np.testing.assert_almost_equal(cube.mad_std(axis=0).value, result)
+
+        mcube = cube.with_mask(cube < 0.98*u.K)
+
+        result2 = np.array([[0.15509701,  0.45763670],
+                            [0.55907956,  0.23835865],
+                            [0.48819454,  0.25499305]])
+
+        np.testing.assert_almost_equal(mcube.mad_std(axis=0).value, result2)