diff --git a/AUTHORS.rst b/AUTHORS.rst index 000a7028..0d0377e0 100644 --- a/AUTHORS.rst +++ b/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) diff --git a/CHANGES.rst b/CHANGES.rst index f4cf7b17..4c0b6d68 100644 --- a/CHANGES.rst +++ b/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 ^^^^^^^^^ diff --git a/ccdproc/combiner.py b/ccdproc/combiner.py index 8f1eec7d..c8476677 100644 --- a/ccdproc/combiner.py +++ b/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: diff --git a/ccdproc/tests/test_combiner.py b/ccdproc/tests/test_combiner.py index 222051f6..4c4efb47 100644 --- a/ccdproc/tests/test_combiner.py +++ b/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() diff --git a/setup.cfg b/setup.cfg index 4f0e3907..9d1c4e75 100644 --- a/setup.cfg +++ b/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