Merge b363055 into b5720f6

docs/manipulation.rst

@@ -85,6 +85,9 @@ that takes the spectrum and an astropy 1D kernel.  So, one could also do:
 In this case, the ``spec1_bsmooth2`` result should be equivalent to the ``spec1_bsmooth`` in
 the section above (assuming the flux data of the input ``spec`` is the same).
 
+The errors are propagated using a standrard "propagation of errors" method, if the uncertainty 
+is defined for the spectrum *and* it is one of StdDevUncertainty, VarianceUncertainty or InverseVariance.
+*Note:* this does *not* consider covariance between points.
 
 Median Smoothing
 ^^^^^^^^^^^^^^^^
@@ -93,7 +96,7 @@ The median based smoothing  is implemented using `scipy.signal.medfilt` and
 has a similar call structure to the convolution-based smoothing methods. This
 method applys the median filter across the flux.
 
-Note: This method is not flux conserving.
+Note: This method is not flux conserving and errors are not propagated.
 
 .. code-block:: python
 

specutils/manipulation/smoothing.py

@@ -1,6 +1,12 @@
-import astropy.units as u
+import copy
+import warnings
+
 from astropy import convolution
+from astropy.nddata import StdDevUncertainty, VarianceUncertainty, InverseVariance
+import astropy.units as u
+from astropy.utils.exceptions import AstropyUserWarning
 from scipy.signal import medfilt
+import numpy as np
 
 from ..spectra import Spectrum1D
 
@@ -15,6 +21,10 @@ def convolution_smooth(spectrum, kernel):
 
     This method can be used along but also is used by other specific methods below.
 
+    If the spectrum uncertainty exists and is StdDevUncertainty, VarianceUncertainty or InverseVariance
+    then the errors will be propagated through the convolution using a standard propagation of errors. The
+    covariance is not considered, currently.
+
     Parameters
     ----------
     spectrum : `~specutils.Spectrum1D`
@@ -43,11 +53,51 @@ def convolution_smooth(spectrum, kernel):
     # Smooth based on the input kernel
     smoothed_flux = convolution.convolve(flux, kernel)
 
+    # Propagate the uncertainty if it exists...
+    uncertainty = copy.deepcopy(spectrum.uncertainty)
+    if uncertainty is not None:
+        if isinstance(uncertainty, StdDevUncertainty):
+            # Convert
+            values = uncertainty.array
+            ivar_values = 1 / values**2
+
+            # Propagate
+            prop_ivar_values = convolution.convolve(ivar_values, kernel)
+
+            # Put back in
+            uncertainty.array = 1 / np.sqrt(prop_ivar_values)
+
+        elif isinstance(uncertainty, VarianceUncertainty):
+            # Convert
+            values = uncertainty.array
+            ivar_values = 1 / values
+
+            # Propagate
+            prop_ivar_values = convolution.convolve(ivar_values, kernel)
+
+            # Put back in
+            uncertainty.array = 1 / prop_ivar_values
+
+        elif isinstance(uncertainty, InverseVariance):
+            # Convert
+            ivar_values = uncertainty.array
+
+            # Propagate
+            prop_ivar_values = convolution.convolve(ivar_values, kernel)
+
+            # Put back in
+            uncertainty.array = prop_ivar_values
+        else:
+            warnings.warn("Uncertainty is {} but convolutional error propagation is not defined for that type.".format(type(uncertainty)),
+                          AstropyUserWarning)
+
+
     # Return a new object with the smoothed flux.
     return Spectrum1D(flux=u.Quantity(smoothed_flux, spectrum.unit),
                       spectral_axis=u.Quantity(spectrum.spectral_axis,
                                                spectrum.spectral_axis_unit),
                       wcs=spectrum.wcs,
+                      uncertainty=uncertainty,
                       velocity_convention=spectrum.velocity_convention,
                       rest_value=spectrum.rest_value)
 

specutils/tests/test_smoothing.py

@@ -3,6 +3,7 @@
 from astropy import convolution
 from scipy.signal import medfilt
 import astropy.units as u
+from astropy.nddata import StdDevUncertainty, VarianceUncertainty, InverseVariance
 from ..spectra.spectrum1d import Spectrum1D
 from ..tests.spectral_examples import simulated_spectra
 
@@ -217,3 +218,42 @@ def test_smooth_median_bad(simulated_spectra, width):
     # Test bad parameters
     with pytest.raises(ValueError):
         median_smooth(spec1, width)
+
+
+def test_smooth_custom_kernel_uncertainty(simulated_spectra):
+    """
+    Test CustomKernel smoothing with correct parmaeters.
+    """
+
+    np.random.seed(42)
+
+    # Create a custom kernel (some weird asymmetric-ness)
+    numpy_kernel = np.array([0.5, 1, 2, 0.5, 0.2])
+    numpy_kernel = numpy_kernel / np.sum(numpy_kernel)
+    custom_kernel = convolution.CustomKernel(numpy_kernel)
+
+    spec1 = simulated_spectra.s1_um_mJy_e1
+    uncertainty = np.abs(np.random.random(spec1.flux.shape))
+
+    # Test StdDevUncertainty
+    spec1.uncertainty = StdDevUncertainty(uncertainty)
+
+    spec1_smoothed = convolution_smooth(spec1, custom_kernel)
+    tt = convolution.convolve(1/(spec1.uncertainty.array**2), custom_kernel)
+    uncertainty_smoothed_astropy = 1/np.sqrt(tt)
+
+    assert np.allclose(spec1_smoothed.uncertainty.array, uncertainty_smoothed_astropy)
+
+    # Test VarianceUncertainty
+    spec1.uncertainty = VarianceUncertainty(uncertainty)
+
+    spec1_smoothed = convolution_smooth(spec1, custom_kernel)
+    uncertainty_smoothed_astropy = 1/convolution.convolve(1/spec1.uncertainty.array, custom_kernel)
+    assert np.allclose(spec1_smoothed.uncertainty.array, uncertainty_smoothed_astropy)
+
+    # Test InverseVariance
+    spec1.uncertainty = InverseVariance(uncertainty)
+
+    spec1_smoothed = convolution_smooth(spec1, custom_kernel)
+    uncertainty_smoothed_astropy = convolution.convolve(spec1.uncertainty.array, custom_kernel)
+    assert np.allclose(spec1_smoothed.uncertainty.array, uncertainty_smoothed_astropy)