Merge bfd8fc8 into dadab3b

specutils/analysis/flux.py

@@ -5,9 +5,11 @@
 import numpy as np
 from ..manipulation import extract_region
 from .utils import computation_wrapper
+from astropy.stats import sigma_clip
+from astropy.stats import median_absolute_deviation
 
 
-__all__ = ['line_flux', 'equivalent_width']
+__all__ = ['line_flux', 'equivalent_width', 'is_continuum_near_zero']
 
 
 def line_flux(spectrum, regions=None):
@@ -116,3 +118,46 @@ 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_near_zero(spectrum, eps=0.01):
+    """
+    Determine if the baseline continuum of the spectrum is near zero.
+
+    The value is scaled to match the sigma/standard deviation parameter of a
+    standard Gaussian profile. This will be calculated over the regions, if
+    they are specified.
+
+    Parameters
+    ----------
+    spectrum : `~specutils.spectra.spectrum1d.Spectrum1D`
+        The spectrum object over which the width will be calculated.
+
+    eps: float
+        The tolerance on the quantification to confirm the continuum is
+        near zero.
+
+    Returns
+    -------
+    is_near_zero: bool
+        A True if the continuum of the spectrum is with eps, False otherwise.
+
+    Notes
+    -----
+    The spectrum should be continuum subtracted before being passed to this
+    function.
+    """
+
+    # If the eps 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(eps, 'unit'):
+        return np.median(spectrum.flux) < eps
+
+    # If eps 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 hasattr(spectrum, 'uncertainty') and spectrum.uncertainty:
+        return np.median(spectrum.flux / spectrum.uncertainty.quantity) < eps
+    else:
+        raise Exception('Spectrum flux has units, eps does not, either include uncertainty or use units on eps')
+        #return np.median(spectrum.flux) / median_absolute_deviation(spectrum.flux) < eps

specutils/tests/test_analysis.py

@@ -10,7 +10,7 @@
 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_near_zero)
 from ..tests.spectral_examples import simulated_spectra
 
 
@@ -578,3 +578,27 @@ 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_near_zero():
+
+    # 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_near_zero(spectrum, eps=0.1*u.Jy)
+
+#    # No mask, no uncertainty, eps 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_near_zero(spectrum, eps=0.1)
+
+    # No mask, with uncertainty
+    wavelengths = [300, 500, 1000] * u.nm
+    data = [0.03, 0.029, 0.033] * u.Jy
+    uncertainty = StdDevUncertainty([1.01, 1.13, 1.1] * u.Jy)
+    spectrum = Spectrum1D(spectral_axis=wavelengths, flux=data,
+                          uncertainty=uncertainty)
+
+    assert True == is_continuum_near_zero(spectrum, eps=0.1)