Merge pull request #538 from brechmos-stsci/i535-check-continuum-subt…

…racted

Simple utility continuum baseline checker

specutils/__init__.py

@@ -7,6 +7,7 @@
 # should keep this content at the top.
 # ----------------------------------------------------------------------------
 from ._astropy_init import *
+from astropy import config as _config
 # ----------------------------------------------------------------------------
 
 # Enforce Python version check during package import.
@@ -34,3 +35,16 @@ class UnsupportedPythonError(Exception):
     _load_user_io()
 
 __citation__ = 'https://doi.org/10.5281/zenodo.1421356'
+
+class Conf(_config.ConfigNamespace):
+    """
+    Configuration parameters for specutils.
+    """
+
+    do_continuum_function_check = _config.ConfigItem(
+        True,
+        'Whether to check the spectrum baseline value is close'
+        'to zero. If it is not within ``threshold`` then a warning is raised.'
+    )
+
+conf = Conf()

specutils/analysis/flux.py

@@ -2,12 +2,21 @@
 A module for analysis tools focused on determining fluxes of spectral features.
 """
 
+from functools import wraps
+import warnings
+
 import numpy as np
+from .. import conf
 from ..manipulation import extract_region
 from .utils import computation_wrapper
+import astropy.units as u
+from astropy.stats import sigma_clip
+from astropy.stats import mad_std
+from astropy.utils.exceptions import AstropyUserWarning
 
 
-__all__ = ['line_flux', 'equivalent_width']
+__all__ = ['line_flux', 'equivalent_width', 'is_continuum_below_threshold',
+           'warn_continuum_below_threshold']
 
 
 def line_flux(spectrum, regions=None):
@@ -116,3 +125,84 @@ def _compute_equivalent_width(spectrum, continuum=1, regions=None):
     ew =  dx - (line_flux / continuum)
 
     return ew.to(calc_spectrum.spectral_axis.unit)
+
+def is_continuum_below_threshold(spectrum, threshold=0.01):
+    """
+    Determine if the baseline of this spectrum is less than a threshold.
+    I.e., an estimate of whether or not the continuum has been subtracted.
+
+    If ``threshold`` is an `~astropy.units.Quantity` with flux units, this
+    directly compares the median of the spectrum to the threshold.
+    of the flux to the threshold.
+
+    If the threshold is a float or dimensionless quantity then the spectrum's uncertainty will be
+    used or an estimate of the uncertainty. If the uncertainty is present then the
+    threshold is compared to the median of the flux divided by the
+    uncertainty.  If the uncertainty is not present then the threshold
+    is compared to the median of the flux divided by the
+    `~astropy.stats.mad_std`.
+
+    Parameters
+    ----------
+    spectrum : `~specutils.spectra.spectrum1d.Spectrum1D`
+        The spectrum object over which the width will be calculated.
+
+    threshold: float or `~astropy.units.Quantity`
+        The tolerance on the quantification to confirm the continuum is
+        near zero.
+
+    Returns
+    -------
+    is_continuum_below_threshold: bool
+        Return True if the continuum of the spectrum is below the threshold, False otherwise.
+
+    """
+
+    flux = spectrum.flux
+    uncertainty = spectrum.uncertainty if hasattr(spectrum, 'uncertainty') else None
+
+    # Apply the mask if it exists.
+    if hasattr(spectrum, 'mask') and spectrum.mask is not None:
+        flux = flux[~spectrum.mask]
+        uncertainty = uncertainty[~spectrum.mask] if uncertainty else uncertainty
+
+    # If the threshold has units then the assumption is that we want
+    # to compare the median of the flux regardless of the
+    # existence of the uncertainty.
+    if hasattr(threshold, 'unit') and not threshold.unit == u.dimensionless_unscaled:
+        return np.median(flux) < threshold
+
+    # If threshold does not have a unit, ie it is not a quantity, then
+    # we are going to calculate based on the S/N if the uncertainty
+    # exists.
+    if uncertainty and uncertainty.uncertainty_type != 'std':
+        return np.median(flux / uncertainty.quantity) < threshold
+    else:
+        return np.median(flux) / mad_std(flux) < threshold
+
+def warn_continuum_below_threshold(threshold=0.01):
+    """
+    Decorator for methods that should warn if the baseline
+    of the spectrum does not appear to be below a threshold.
+
+    The ``check`` parameter is based on the
+    `astropy configuration system <http://docs.astropy.org/en/stable/config/#adding-new-configuration-items>`_.
+    Examples are on that page to show how to turn off this type of warning checking.
+    """
+    def actual_decorator(function):
+        @wraps(function)
+        def wrapper(*args, **kwargs):
+            if conf.do_continuum_function_check:
+                spectrum = args[0]
+                if not is_continuum_below_threshold(spectrum, threshold):
+                    message = "Spectrum is not below the threshold {}.\n\n".format(threshold)
+                    message += ("""If you want to suppress this warning either type """
+                                """'specutils.conf.do_continuum_function_check = False' or """
+                                """see http://docs.astropy.org/en/stable/config/#adding-new-configuration-items """
+                                """for other ways to configure the warning.""")
+
+                    warnings.warn(message, AstropyUserWarning)
+            result = function(*args, **kwargs)
+            return result
+        return wrapper
+    return actual_decorator

specutils/fitting/fitmodels.py

@@ -8,7 +8,7 @@
 from astropy.stats import sigma_clipped_stats
 
 from ..manipulation.utils import excise_regions
-from ..analysis import fwhm, centroid
+from ..analysis import fwhm, centroid, warn_continuum_below_threshold
 from ..utils import QuantityModel
 from ..manipulation import extract_region, noise_region_uncertainty
 from ..spectra.spectral_region import SpectralRegion
@@ -97,6 +97,7 @@ def _consecutive(data, stepsize=1):
     return np.split(data, np.where(np.diff(data) != stepsize)[0]+1)
 
 
+@warn_continuum_below_threshold(threshold=0.01)
 def find_lines_threshold(spectrum, noise_factor=1):
     """
     Find the emission and absorption lines in a spectrum. The method
@@ -170,6 +171,7 @@ def find_lines_threshold(spectrum, noise_factor=1):
     return qtable
 
 
+@warn_continuum_below_threshold(threshold=0.01)
 def find_lines_derivative(spectrum, flux_threshold=None):
     """
     Find the emission and absorption lines in a spectrum. The method

specutils/tests/test_analysis.py

@@ -6,11 +6,13 @@
 from astropy.nddata import StdDevUncertainty
 from astropy.stats.funcs import gaussian_sigma_to_fwhm
 from astropy.tests.helper import quantity_allclose
+from astropy.utils.exceptions import AstropyUserWarning
 
 from ..spectra import Spectrum1D, SpectralRegion
 from ..analysis import (line_flux, equivalent_width, snr, centroid,
                         gaussian_sigma_width, gaussian_fwhm, fwhm,
-                        snr_derived, fwzi)
+                        snr_derived, fwzi, is_continuum_below_threshold)
+from ..fitting import find_lines_threshold
 from ..tests.spectral_examples import simulated_spectra
 
 
@@ -578,3 +580,50 @@ def test_fwzi_multi_spectrum():
     expected = [113.51706001 * u.AA, 567.21252727 * u.AA, 499.5024546 * u.AA]
 
     assert quantity_allclose(fwzi(spec), expected)
+
+
+def test_is_continuum_below_threshold():
+
+    # No mask, no uncertainty
+    wavelengths = [300, 500, 1000] * u.nm
+    data = [0.001, -0.003, 0.003] * u.Jy
+    spectrum = Spectrum1D(spectral_axis=wavelengths, flux=data)
+    assert True == is_continuum_below_threshold(spectrum, threshold=0.1*u.Jy)
+
+#    # No mask, no uncertainty, threshold is float
+#    wavelengths = [300, 500, 1000] * u.nm
+#    data = [0.0081, 0.0043, 0.0072] * u.Jy
+#    spectrum = Spectrum1D(spectral_axis=wavelengths, flux=data)
+#    assert True == is_continuum_below_threshold(spectrum, threshold=0.1)
+
+    # No mask, with uncertainty
+    wavelengths = [300, 500, 1000] * u.nm
+    data = [0.03, 0.029, 0.031] * u.Jy
+    uncertainty = StdDevUncertainty([1.01, 1.03, 1.01] * u.Jy)
+    spectrum = Spectrum1D(spectral_axis=wavelengths, flux=data,
+                          uncertainty=uncertainty)
+
+    assert True == is_continuum_below_threshold(spectrum, threshold=0.1*u.Jy)
+
+    # With mask, with uncertainty
+    wavelengths = [300, 500, 1000] * u.nm
+    data = [0.01, 1.029, 0.013] * u.Jy
+    mask = np.array([False, True, False])
+    uncertainty = StdDevUncertainty([1.01, 1.13, 1.1] * u.Jy)
+    spectrum = Spectrum1D(spectral_axis=wavelengths, flux=data,
+                          uncertainty=uncertainty, mask=mask)
+
+    assert True == is_continuum_below_threshold(spectrum, threshold=0.1*u.Jy)
+
+    # With mask, with uncertainty -- should throw an exception
+    wavelengths = [300, 500, 1000] * u.nm
+    data = [10.03, 10.029, 10.033] * u.Jy
+    mask = np.array([False, False, False])
+    uncertainty = StdDevUncertainty([1.01, 1.13, 1.1] * u.Jy)
+    spectrum = Spectrum1D(spectral_axis=wavelengths, flux=data,
+                          uncertainty=uncertainty, mask=mask)
+
+    print('spectrum has flux {}'.format(spectrum.flux))
+    with pytest.warns(AstropyUserWarning) as e_info:
+        find_lines_threshold(spectrum, noise_factor=1)
+        assert len(e_info)==1 and 'if you want to suppress this warning' in e_info[0].message.args[0].lower()