Always use astropy's sigma clipping and default to string-based functions

AUTHORS.rst

@@ -49,6 +49,7 @@ Alphabetical list of code contributors
 * Adrian Price-Whelan (@adrn)
 * JVSN Reddy (@janga1997)
 * Luca Rizzi (@lucarizzi)
+* Thomas Robitaille (@astrofrog)
 * Evert Rol (@evertrol)
 * Jenna Ryon (@jryon)
 * William Schoenell (@wschoenell)

CHANGES.rst

@@ -1,4 +1,4 @@
-3.0.0 (unreleased)
+2.4.0 (unreleased)
 ------------------
 
 New Features
@@ -7,6 +7,11 @@ New Features
 Other Changes and Additions
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
+- The sigma clipping option in the image combiner now always uses the
+  astropy sigma clipping function, and supports specifying the
+  functions to use for estimating the center and deviation values
+  as strings for common cases (which significantly improves performance). [#794]
+
 Bug Fixes
 ^^^^^^^^^
 

ccdproc/combiner.py

@@ -16,6 +16,7 @@
 
 from astropy.nddata import CCDData, StdDevUncertainty
 from astropy.stats import sigma_clip
+from astropy.utils import deprecated_renamed_argument
 from astropy import log
 
 __all__ = ['Combiner', 'combine']
@@ -294,9 +295,13 @@ def minmax_clipping(self, min_clip=None, max_clip=None):
             self.data_arr.mask[mask] = True
 
     # set up sigma  clipping algorithms
+    @deprecated_renamed_argument('use_astropy', None, arg_in_kwargs=True,
+                                 since='2.4.0',
+                                 message='The use_astropy argument has been removed because '
+                                         'astropy sigma clipping is now always used.'
+                                 )
     def sigma_clipping(self, low_thresh=3, high_thresh=3,
-                       func=ma.mean, dev_func=ma.std, use_astropy=False,
-                       **kwd):
+                       func='mean', dev_func='std', **kwd):
         """
         Pixels will be rejected if they have deviations greater than those
         set by the threshold values. The algorithm will first calculated
@@ -307,6 +312,7 @@ def sigma_clipping(self, low_thresh=3, high_thresh=3,
 
         Parameters
         -----------
+
         low_thresh : positive float or None, optional
             Threshold for rejecting pixels that deviate below the baseline
             value. If negative value, then will be convert to a positive
@@ -318,54 +324,40 @@ def sigma_clipping(self, low_thresh=3, high_thresh=3,
             value. If None, no rejection will be done based on high_thresh.
             Default is 3.
 
-        func : function, optional
-            Function for calculating the baseline values (i.e. `numpy.ma.mean`
-            or `numpy.ma.median`). This should be a function that can handle
-            `numpy.ma.MaskedArray` objects. **Set to ``'median'`` and
-            set ``use_astropy=True`` for best performance if using a
-            median.**
-            Default is `numpy.ma.mean`.
-
-        dev_func : function, optional
-            Function for calculating the deviation from the baseline value
-            (i.e. `numpy.ma.std`). This should be a function that can handle
-            `numpy.ma.MaskedArray` objects.
-            Default is `numpy.ma.std`.
-
-        use_astropy : bool, optional
-            If ``True``, use astropy's `~astropy.stats.sigma_clip`, which is faster
-            and more flexible. The high/low sigma clip parameters are set
-            from ``low_thresh`` and ``high_thresh``. Any remaining keywords are passed
-            in to astropy's `~astropy.stats.sigma_clip`. By default, the
-            number of iterations and other settings will be made to reproduce
-            the behavior of ccdproc's ``sigma_clipping``.
+        func : {'median', 'mean'} or callable, optional
+            The statistic or callable function/object used to compute
+            the center value for the clipping. If using a callable
+            function/object and the ``axis`` keyword is used, then it must
+            be able to ignore NaNs (e.g., `numpy.nanmean`) and it must have
+            an ``axis`` keyword to return an array with axis dimension(s)
+            removed. The default is ``'median'``.
+
+        dev_func : {'std', 'mad_std'} or callable, optional
+            The statistic or callable function/object used to compute the
+            standard deviation about the center value. If using a callable
+            function/object and the ``axis`` keyword is used, then it must
+            be able to ignore NaNs (e.g., `numpy.nanstd`) and it must have
+            an ``axis`` keyword to return an array with axis dimension(s)
+            removed. The default is ``'std'``.
+
+        kwd
+            Any remaining keyword arguments are passed to astropy's
+            :func:`~astropy.stats.sigma_clip` function.
         """
-        if use_astropy:
-            copy = kwd.get('copy', False)
-            axis = kwd.get('axis', 0)
-            maxiters = kwd.get('maxiters', 1)
-            self.data_arr.mask = \
-                sigma_clip(self.data_arr.data, sigma_lower=low_thresh,
-                           sigma_upper=high_thresh, axis=axis, copy=copy,
-                           maxiters=maxiters,
-                           cenfunc=func, stdfunc=dev_func,
-                           masked=True,
-                           **kwd).mask
-            return
-
-        # setup baseline values
-        baseline = func(self.data_arr, axis=0)
-        dev = dev_func(self.data_arr, axis=0)
-        # reject values
-        if low_thresh is not None:
-            # check for negative numbers in low_thresh
-            if low_thresh < 0:
-                low_thresh = abs(low_thresh)
-            mask = (self.data_arr - baseline < -low_thresh * dev)
-            self.data_arr.mask[mask] = True
-        if high_thresh is not None:
-            mask = (self.data_arr - baseline > high_thresh * dev)
-            self.data_arr.mask[mask] = True
+
+        # Remove in 3.0
+        _ = kwd.pop('use_astropy', True)
+
+        self.data_arr.mask = sigma_clip(self.data_arr.data,
+                                        sigma_lower=low_thresh,
+                                        sigma_upper=high_thresh,
+                                        axis=kwd.get('axis', 0),
+                                        copy=kwd.get('copy', False),
+                                        maxiters=kwd.get('maxiters', 1),
+                                        cenfunc=func,
+                                        stdfunc=dev_func,
+                                        masked=True,
+                                        **kwd).mask
 
     def _get_scaled_data(self, scale_arg):
         if scale_arg is not None:

ccdproc/tests/test_combiner.py

@@ -227,21 +227,6 @@ def test_combiner_sigmaclip_low():
     assert c.data_arr[5].mask.all()
 
 
-@pytest.mark.parametrize('threshold', [1, 10])
-def test_combiner_sigma_clip_use_astropy_same_result(threshold):
-    # If we turn on use_astropy and make no other changes we should get exactly
-    # the same result as if we use ccdproc sigma_clipping
-    ccd_list = [ccd_data_func(rng_seed=seed + 1) for seed in range(10)]
-    c_ccdp = Combiner(ccd_list)
-    c_apy = Combiner(ccd_list)
-
-    c_ccdp.sigma_clipping(low_thresh=threshold, high_thresh=threshold)
-    c_apy.sigma_clipping(low_thresh=threshold, high_thresh=threshold,
-                         use_astropy=True)
-
-    np.testing.assert_allclose(c_ccdp.data_arr.mask, c_apy.data_arr.mask)
-
-
 # test that the median combination works and returns a ccddata object
 def test_combiner_median():
     ccd_data = ccd_data_func()

setup.cfg

@@ -30,7 +30,7 @@ packages = find:
 zip_safe = False
 setup_requires = setuptools_scm
 install_requires = numpy>=1.18
-                   astropy>=4.0.6  # Support LTS, but only with bug fixes
+                   astropy>=4.3
                    scipy
                    astroscrappy>=1.0.8
                    reproject>=0.7