Added in spectrum_from_model

specutils/manipulation/__init__.py

@@ -1,3 +1,4 @@
 from .smoothing import *  # noqa
 from .estimate_uncertainty import *  # noqa
 from .extract_spectral_region import *  # noqa
+from .utils import *  # noqa

specutils/manipulation/utils.py

@@ -2,13 +2,13 @@
 
 from ..spectra import Spectrum1D, SpectralRegion
 
-__all__ = ['excise_regions', 'linear_exciser']
+__all__ = ['excise_regions', 'linear_exciser', 'spectrum_from_model']
 
 
 def linear_exciser(spectrum, region):
     """
     Basic spectral excise method where the spectral region defined by the
-    2-tuple parameter `region` (start and end wavelengths) will result
+    2-tuple parameter ``region`` (start and end wavelengths) will result
     in the flux between those regions set to a linear ramp of the
     two points immediately before and after the start and end regions.
 
@@ -30,12 +30,12 @@ def linear_exciser(spectrum, region):
     Raises
     ------
     ValueError
-       In the case that ``spectrum`` and ``regions`` are not the correct types.
+       In the case that ``spectrum`` and ``region`` are not the correct types.
 
     """
 
     #
-    # Find the indices of the wavelengths in the range `range`
+    # Find the indices of the wavelengths in the range ``range``
     #
 
     wavelengths = spectrum.spectral_axis
@@ -77,7 +77,7 @@ def excise_regions(spectrum, regions, exciser=linear_exciser):
 
     regions : list of `~specutils.SpectralRegion`
         Each element of the list is a `~specutils.SpectralRegion`. The flux
-        between these wavelengths will be "cut out" using the `exciser`
+        between these wavelengths will be "cut out" using the ``exciser``
         method.
 
     exciser: method
@@ -151,3 +151,42 @@ def excise_region(spectrum, region, exciser=linear_exciser):
     #
 
     return exciser(spectrum, region)
+
+
+def spectrum_from_model(model_input, spectrum):
+    """
+    This method will create a `~specutils.Spectrum1D` object
+    with the flux defined by calling the input ``model``. All
+    other parameters for the output `~specutils.Spectrum1D` object
+    will be the same as the input `~specutils.Spectrum1D` object.
+
+    Parameters
+    ----------
+    model : `~astropy.modeling.Model`
+        The input model or compound model from which flux is calculated.
+
+    spectrum : `~specutils.Spectrum1D`
+        The `~specutils.Spectrum1D` object to use as the model template.
+
+    Returns
+    -------
+    spectrum : `~specutils.Spectrum1D`
+        Output `~specutils.Spectrum1D` which is copy of the one passed in with the updated flux.
+        The uncertainty will not be copied as it is not necessarily the same.
+
+    """
+
+    # If the input model has units then we will call it normally.
+    if getattr(model_input, model_input.param_names[0]).unit is not None:
+        flux = model_input(spectrum.spectral_axis)
+
+    # If the input model does not have units, then assume it is in
+    # the same units as the input spectrum.
+    else:
+        flux = model_input(spectrum.spectral_axis.value)*spectrum.flux.unit
+
+    return Spectrum1D(flux=flux,
+                      spectral_axis=spectrum.spectral_axis,
+                      wcs=spectrum.wcs,
+                      velocity_convention=spectrum.velocity_convention,
+                      rest_value=spectrum.rest_value)

specutils/tests/test_fitting.py

@@ -7,7 +7,7 @@
 from ..fitting import (fit_lines, find_lines_derivative,
                        find_lines_threshold, estimate_line_parameters)
 from ..analysis import fwhm, centroid
-from ..manipulation import noise_region_uncertainty
+from ..manipulation import noise_region_uncertainty, spectrum_from_model
 
 
 def single_peak():
@@ -512,3 +512,45 @@ def test_fitter_parameters():
                2.88900005e-01, 6.24602556e-02, 9.22061121e-03, 9.29427266e-04]) * u.Jy
 
     assert np.allclose(y_single_fit.value[::10], y_single_fit_expected.value, atol=1e-5)
+
+
+def test_spectrum_from_model():
+    """
+    This test fits the the first simulated spectrum from the fixture.  The
+    initial guesses are manually set here with bounds that essentially make
+    sense as the functionality of the test is to make sure the fit works and
+    we get a reasonable answer out **given** good initial guesses.
+    """
+
+    np.random.seed(0)
+    x = np.linspace(0., 10., 200)
+    y = 3 * np.exp(-0.5 * (x - 6.3)**2 / 0.1**2)
+    y += np.random.normal(0., 0.2, x.shape)
+
+    y_continuum = 3.2 * np.exp(-0.5 * (x - 5.6)**2 / 4.8**2)
+    y += y_continuum
+
+    spectrum = Spectrum1D(flux=y*u.Jy, spectral_axis=x*u.um)
+
+    # Unitless test
+    chebyshev = models.Chebyshev1D(3, c0=0.1, c1=4, c2=5)
+    spectrum_chebyshev = spectrum_from_model(chebyshev, spectrum)
+
+    flux_expected = np.array([-4.90000000e+00, -3.64760991e-01,  9.22085553e+00,  2.38568496e+01,
+                              4.35432211e+01,  6.82799702e+01,  9.80670968e+01,  1.32904601e+02,
+                              1.72792483e+02,  2.17730742e+02,  2.67719378e+02,  3.22758392e+02,
+                              3.82847784e+02,  4.47987553e+02,  5.18177700e+02,  5.93418224e+02,
+                              6.73709126e+02,  7.59050405e+02,  8.49442062e+02,  9.44884096e+02])
+
+    assert np.allclose(spectrum_chebyshev.flux.value[::10], flux_expected, atol=1e-5)
+
+    # Unitfull test
+    gaussian = models.Gaussian1D(amplitude=5*u.Jy, mean=4*u.um, stddev=2.3*u.um)
+    spectrum_gaussian = spectrum_from_model(gaussian, spectrum)
+
+    flux_expected = np.array([1.1020263, 1.57342489, 2.14175093, 2.77946243, 3.4389158,
+                              4.05649712, 4.56194132, 4.89121902, 4.99980906, 4.872576,
+                              4.52723165, 4.01028933, 3.3867847, 2.72689468, 2.09323522,
+                              1.5319218, 1.06886794, 0.71101768, 0.45092638, 0.27264641])
+
+    assert np.allclose(spectrum_gaussian.flux.value[::10], flux_expected, atol=1e-5)