Manipulating Spectra¶
While there are myriad ways you might want to alter a spectrum, specutils provides some specific functionality that is commonly used in astronomy. These tools are detailed here, but it is important to bear in mind that this is not intended to be exhaustive - the point of specutils is to provide a framework you can use to do your data analysis. Hence the functionality described here is best thought of as pieces you might string together with your own functionality to build a tailor-made spectral analysis environment.
In general, however, specutils is designed around the idea that spectral manipulations generally yield new spectrum objects, rather than in-place operations. This is not a true restriction, but is a guideline that is recommended primarily to keep you from accidentally modifying a spectrum you didn’t mean to change.
Smoothing¶
Specutils provides smoothing for spectra in two forms: 1) convolution based
using smoothing astropy.convolution and 2) median filtering
using the scipy.signal.medfilt(). Each of these act on the flux
of the Spectrum1D object.
Note
Specutils smoothing kernel widths and standard deviations are
in units of pixels and not Quantity.
Convolution Based Smoothing¶
While any kernel supported by astropy.convolution will work (using the
convolution_smooth function), several
commonly-used kernels have convenience functions wrapping them to simplify
the smoothing process into a simple one-line operation. Currently
implemented are: box_smooth()
(Box1DKernel),
gaussian_smooth()
(Gaussian1DKernel), and
trapezoid_smooth()
(Trapezoid1DKernel).
>>> from specutils import Spectrum1D
>>> import astropy.units as u
>>> import numpy as np
>>> from specutils.manipulation import (box_smooth, gaussian_smooth, trapezoid_smooth)
>>> spec1 = Spectrum1D(spectral_axis=np.arange(1, 50) * u.nm, flux=np.random.sample(49)*u.Jy)
>>> spec1_bsmooth = box_smooth(spec1, width=3)
>>> spec1_gsmooth = gaussian_smooth(spec1, stddev=3)
>>> spec1_tsmooth = trapezoid_smooth(spec1, width=3)
>>> gaussian_smooth(spec1, stddev=3)
Spectrum1D([0.22830748, 0.2783204 , 0.32007408, 0.35270403, 0.37899655,
0.40347983, 0.42974259, 0.45873436, 0.48875214, 0.51675647,
0.54007149, 0.55764758, 0.57052796, 0.58157173, 0.59448669,
0.61237409, 0.63635755, 0.66494062, 0.69436655, 0.7199299 ,
0.73754271, 0.74463192, 0.74067744, 0.72689092, 0.70569365,
0.6800534 , 0.65262146, 0.62504013, 0.59778884, 0.57072578,
0.54416776, 0.51984003, 0.50066938, 0.48944714, 0.48702192,
0.49126444, 0.49789092, 0.50276877, 0.50438924, 0.50458914,
0.50684731, 0.51321106, 0.52197328, 0.52782086, 0.52392599,
0.50453064, 0.46677128, 0.41125485, 0.34213489])
Each of the specific smoothing methods create the appropriate astropy.convolution.convolve
kernel and then call a helper function convolution_smooth()
that takes the spectrum and an astropy 1D kernel. So, one could also do:
>>> from astropy.convolution import Box1DKernel
>>> from specutils.manipulation import convolution_smooth
>>> box1d_kernel = Box1DKernel(width=3)
>>> spec1 = Spectrum1D(spectral_axis=np.arange(1, 50) * u.nm, flux=np.random.sample(49) * u.Jy)
>>> spec1_bsmooth2 = convolution_smooth(spec1, box1d_kernel)
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 uncertainties are propagated using a standard “propagation of errors” method, if the uncertainty is defined for the spectrum and it is one of StdDevUncertainty, VarianceUncertainty or InverseVariance. But note that this does not consider covariance between points.
Median Smoothing¶
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 and errors are not propagated.
>>> from specutils.manipulation import median_smooth
>>> spec1 = Spectrum1D(spectral_axis=np.arange(1, 50) * u.nm, flux=np.random.sample(49) * u.Jy)
>>> spec1_msmooth = median_smooth(spec1, width=3)
Resampling¶
specutils contains several classes for resampling the flux
in a Spectrum1D object. Currently supported methods of
resampling are integrated flux conserving with FluxConservingResampler,
linear interpolation with LinearInterpolatedResampler,
and cubic spline with SplineInterpolatedResampler.
Each of these classes takes in a Spectrum1D and a user
defined output dispersion grid, and returns a new Spectrum1D
with the resampled flux. Currently the resampling classes expect the new
dispersion grid unit to be the same as the input spectrum’s dispersion grid unit.
If the input Spectrum1D contains an uncertainty,
FluxConservingResampler will propogate the
uncertainty to the final output Spectrum1D. However, the
other two implemented resampling classes (LinearInterpolatedResampler
and SplineInterpolatedResampler) will ignore
any input uncertainty.
Here’s a set of simple examples showing each of the three types of resampling:
Splicing/Combining Multiple Spectra¶
The resampling functionality detailed above is also the default way specutils supports splicing multiple spectra together into a single spectrum. This can be achieved as follows:
Uncertainty Estimation¶
Some of the machinery in specutils (e.g.
snr) requires an uncertainty to be present.
While some data reduction pipelines generate this as part of the
reduction process, sometimes it’s necessary to estimate the
uncertainty in a spectrum using the spectral data itself. Currently
specutils provides the straightforward
noise_region_uncertainty function.
First we build a spectrum like that used in Analysis, but without a known uncertainty:
>>> from astropy.modeling import models
>>> np.random.seed(42)
>>> spectral_axis = np.linspace(10, 1, 200) * u.GHz
>>> spectral_model = models.Gaussian1D(amplitude=3*u.Jy, mean=5*u.GHz, stddev=0.8*u.GHz)
>>> flux = spectral_model(spectral_axis)
>>> flux += np.random.normal(0., 0.2, spectral_axis.shape) * u.Jy
>>> noisy_gaussian = Spectrum1D(spectral_axis=spectral_axis, flux=flux)
Now we estimate the uncertainty from the region that does not contain the line:
>>> from specutils import SpectralRegion
>>> from specutils.manipulation import noise_region_uncertainty
>>> noise_region = SpectralRegion([(10, 7), (3, 0)] * u.GHz)
>>> spec_w_unc = noise_region_uncertainty(noisy_gaussian, noise_region)
>>> spec_w_unc.uncertainty
StdDevUncertainty([0.18823157, ..., 0.18823157])
Or similarly, expressed in pixels:
>>> noise_region = SpectralRegion([(0, 25), (175, 200)]*u.pix)
>>> spec_w_unc = noise_region_uncertainty(noisy_gaussian, noise_region)
>>> spec_w_unc.uncertainty
StdDevUncertainty([0.18739524, ..., 0.18739524])
S/N Threshold Mask¶
It is useful to be able to find all the spaxels in an ND spectrum
in which the signal to noise ratio is greater than some threshold.
This method implements this functionality so that a Spectrum1D
object, SpectrumCollection or an NDData derived
object may be passed in as the first parameter. The second parameter
is a floating point threshold.
For example, first a spectrum with flux and uncertainty is created, and
then call the snr_threshold method:
>>> import numpy as np
>>> from astropy.nddata import StdDevUncertainty
>>> import astropy.units as u
>>> from specutils import Spectrum1D
>>> from specutils.manipulation import snr_threshold
>>> np.random.seed(42)
>>> wavelengths = np.arange(0, 10)*u.um
>>> flux = 100*np.abs(np.random.randn(10))*u.Jy
>>> uncertainty = StdDevUncertainty(np.abs(np.random.randn(10))*u.Jy)
>>> spectrum = Spectrum1D(spectral_axis=wavelengths, flux=flux, uncertainty=uncertainty)
>>> spectrum_masked = snr_threshold(spectrum, 50)
>>> # To create a masked flux array
>>> flux_masked = spectrum_masked.flux
>>> flux_masked[spectrum_masked.mask] = np.nan
The output spectrum_masked is a shallow copy of the input spectrum
with the mask attribute set to False where the S/N is greater than 50
and True elsewhere. It is this way to be consistent with astropy.nddata.
Note
The mask attribute is the only attribute modified by snr_threshold(). To
retrieve the masked flux data use spectrum.masked.flux_masked.
Shifting¶
In addition to resampling, you may sometimes wish to simply shift the
spectral_axis of a spectrum (a la the specshift iraf task).
There is no explicit function for this because it is a basic transform of
the spectral_axis. Therefore one can use a construct like this:
>>> from specutils import Spectrum1D
>>> np.random.seed(42)
>>> wavelengths = np.arange(0, 10) * u.um
>>> flux = 100 * np.abs(np.random.randn(10)) * u.Jy
>>> spectrum = Spectrum1D(spectral_axis=wavelengths, flux=flux)
>>> spectrum
<Spectrum1D(flux=<Quantity [ 49.6714153 , 13.82643012, 64.76885381, 152.30298564,
23.41533747, 23.41369569, 157.92128155, 76.74347292,
46.94743859, 54.25600436] Jy>,
spectral_axis=<SpectralAxis [0., 1., 2., 3., 4., 5., 6., 7., 8., 9.] um,
radial_velocity=0.0 km / s,
redshift=0.0,
doppler_rest=0.0 Angstrom,
doppler_convention=None,
observer=None,
target=None>)>
>>> shift = 12300 * u.AA
>>> new_spec = Spectrum1D(spectral_axis=spectrum.spectral_axis + shift, flux=spectrum.flux)
>>> new_spec
<Spectrum1D(flux=<Quantity [ 49.6714153 , 13.82643012, 64.76885381, 152.30298564,
23.41533747, 23.41369569, 157.92128155, 76.74347292,
46.94743859, 54.25600436] Jy>, spectral_axis=<SpectralAxis [ 1.23, 2.23, 3.23, 4.23, 5.23, 6.23, 7.23, 8.23,
9.23, 10.23] um,
radial_velocity=0.0 km / s,
redshift=0.0,
doppler_rest=0.0 Angstrom,
doppler_convention=None,
observer=None,
target=None>)>