From fd43fe8c5b06edbeff220f1da3160258e126c0f8 Mon Sep 17 00:00:00 2001
From: "Adam Ginsburg (keflavich)" <keflavich@gmail.com>
Date: Fri, 24 Mar 2017 15:51:17 -0600
Subject: [PATCH 1/2] add mad_std and fix up the docs for regular std

---
 spectral_cube/spectral_cube.py | 28 +++++++++++++++++++++++++++-
 1 file changed, 27 insertions(+), 1 deletion(-)

diff --git a/spectral_cube/spectral_cube.py b/spectral_cube/spectral_cube.py
index fb6887a65..6d8393e67 100644
--- a/spectral_cube/spectral_cube.py
+++ b/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
 
@@ -331,7 +332,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
@@ -535,6 +538,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)
@@ -578,6 +594,16 @@ 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
+        """
+        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,
+                                         projection=projection)
+
 
     @aggregation_docstring
     def max(self, axis=None, how='auto'):

From 05d537bd580b07f57ae5d412a4d72dcb49f1f793 Mon Sep 17 00:00:00 2001
From: "Adam Ginsburg (keflavich)" <keflavich@gmail.com>
Date: Sun, 21 May 2017 18:24:42 -0600
Subject: [PATCH 2/2] add mad_std test

---
 spectral_cube/spectral_cube.py            |  3 +++
 spectral_cube/tests/test_spectral_cube.py | 24 +++++++++++++++++++++++
 2 files changed, 27 insertions(+)

diff --git a/spectral_cube/spectral_cube.py b/spectral_cube/spectral_cube.py
index 6d8393e67..e92d8c624 100644
--- a/spectral_cube/spectral_cube.py
+++ b/spectral_cube/spectral_cube.py
@@ -599,9 +599,12 @@ 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)
 
 
diff --git a/spectral_cube/tests/test_spectral_cube.py b/spectral_cube/tests/test_spectral_cube.py
index 758630e8f..e60b6d727 100644
--- a/spectral_cube/tests/test_spectral_cube.py
+++ b/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
@@ -1297,3 +1298,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)