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).
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.
>>> 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)
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 :doc:`analysis`, but without a
known uncertainty:
>>> from astropy.modeling import models
>>> np.random.seed(42)
>>> spectral_axis = np.linspace(0., 10., 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
>>> spec_w_unc = noise_region_uncertainty(noisy_gaussian, SpectralRegion(0*u.GHz, 3*u.GHz))
>>> spec_w_unc.uncertainty # doctest: +ELLIPSIS
StdDevUncertainty([0.2, ..., 0.2]
Reference/API¶
Functions¶
box_smooth (spectrum, width) |
Smooth a Spectrum1D instance based on a astropy.convolution.Box1DKernel kernel. |
convolution_smooth (spectrum, kernel) |
Apply a convolution based smoothing to the spectrum. |
extract_region (spectrum, region) |
Extract a region from the input Spectrum1D defined by the lower and upper bounds defined by this SpectralRegion instance. |
gaussian_smooth (spectrum, stddev) |
Smooth a Spectrum1D instance based on a astropy.convolution.Gaussian1DKernel . |
median_smooth (spectrum, width) |
Smoothing based on a median filter. |
noise_region_uncertainty (spectrum, …[, …]) |
Generates a new spectrum with an uncertainty from the noise in a particular region of the spectrum. |
trapezoid_smooth (spectrum, width) |
Smoothing based on a astropy.convolution.Trapezoid1DKernel kernel. |