Spectral Regions

A spectral region may be defined and may encompass one, or more, sub-regions. They are defined independently of a Spectrum1D object in the sense that spectral regions like “near the Halpha line rest wavelength” have meaning independent of the details of a particular spectrum.

Spectral regions can be defined either as a single region by passing two Quantity’s or by passing a list of 2-tuples. Note that the units of these quantites can be any valid spectral unit or u.pixel (which indicates to use indexing directly).

>>> from astropy import units as u
>>> from specutils.spectra import SpectralRegion

>>> sr = SpectralRegion(0.45*u.um, 0.6*u.um)
>>> sr_two = SpectralRegion([(0.45*u.um, 0.6*u.um), (0.8*u.um, 0.9*u.um)])

SpectralRegion can be combined by using the ‘+’ operator:

>>> from astropy import units as u
>>> from specutils.spectra import SpectralRegion
>>> sr = SpectralRegion(0.45*u.um, 0.6*u.um) + SpectralRegion(0.8*u.um, 0.9*u.um)

Regions can also be added in place:

>>> from astropy import units as u
>>> from specutils.spectra import SpectralRegion

>>> sr1 = SpectralRegion(0.45*u.um, 0.6*u.um)
>>> sr2 = SpectralRegion(0.8*u.um, 0.9*u.um)
>>> sr1 += sr2

Regions can be sliced by indexing by an integer or by a range:

>>> from astropy import units as u
>>> from specutils.spectra import SpectralRegion

>>> sr = SpectralRegion(0.15*u.um, 0.2*u.um) + SpectralRegion(0.3*u.um, 0.4*u.um) +\
...      SpectralRegion(0.45*u.um, 0.6*u.um) + SpectralRegion(0.8*u.um, 0.9*u.um) +\
...      SpectralRegion(1.0*u.um, 1.2*u.um) + SpectralRegion(1.3*u.um, 1.5*u.um)

>>> # Get one spectral region (returns a SpectralRegion instance)
>>> sone = sr1[0]

>>> # Slice spectral region.
>>> subsr = sr[3:5]
>>> # SpectralRegion: 0.8 um - 0.9 um, 1.0 um - 1.2 um

The lower and upper bounds on a region are accessed by calling lower or upper. The lower bound of a SpectralRegion is the minimum of the lower bounds of each sub-region and the upper bound is the maximum of the upper bounds:

>>> from astropy import units as u
>>> from specutils.spectra import SpectralRegion

>>> sr = (SpectralRegion(0.15*u.um, 0.2*u.um) + SpectralRegion(0.3*u.um, 0.4*u.um) +
...       SpectralRegion(0.45*u.um, 0.6*u.um) + SpectralRegion(0.8*u.um, 0.9*u.um) +
...       SpectralRegion(1.0*u.um, 1.2*u.um) + SpectralRegion(1.3*u.um, 1.5*u.um))

>>> # Bounds on the spectral region (most minimum and maximum bound)
>>> sr.bounds
(<Quantity 0.15 um>, <Quantity 1.5 um>)

>>> # Lower bound on the spectral region (most minimum)
>>> sr.lower
<Quantity 0.15 um>

>>> sr.upper
<Quantity 1.5 um>

>>> # Lower bound on one element of the spectral region.
>>> sr[3].lower
<Quantity 0.8 um>

One can also delete a sub-region:

>>> from astropy import units as u
>>> from specutils.spectra import SpectralRegion

>>> sr = (SpectralRegion(0.15*u.um, 0.2*u.um) + SpectralRegion(0.3*u.um, 0.4*u.um) +
...       SpectralRegion(0.45*u.um, 0.6*u.um) + SpectralRegion(0.8*u.um, 0.9*u.um) +
...       SpectralRegion(1.0*u.um, 1.2*u.um) + SpectralRegion(1.3*u.um, 1.5*u.um))

>>> del sr[1]
>>> sr
Spectral Region, 5 sub-regions:
(0.15 um, 0.2 um)   (0.45 um, 0.6 um)   (0.8 um, 0.9 um)
(1.0 um, 1.2 um)    (1.3 um, 1.5 um)

There is also the ability to iterate:

>>> from astropy import units as u
>>> from specutils.spectra import SpectralRegion

>>> sr = (SpectralRegion(0.15*u.um, 0.2*u.um) + SpectralRegion(0.3*u.um, 0.4*u.um) +
...       SpectralRegion(0.45*u.um, 0.6*u.um) + SpectralRegion(0.8*u.um, 0.9*u.um) +
...       SpectralRegion(1.0*u.um, 1.2*u.um) + SpectralRegion(1.3*u.um, 1.5*u.um))

>>> for s in sr:
...     print(s.lower)
0.15 um
0.3 um
0.45 um
0.8 um
1.0 um
1.3 um

And, lastly, there is the ability to invert a SpectralRegion given a lower and upper bound. For example, if a set of ranges are defined each defining a range around lines, then calling invert will return a SpectralRegion that defines the baseline/noise regions:

>>> from astropy import units as u
>>> from specutils.spectra import SpectralRegion

>>> sr = (SpectralRegion(0.15*u.um, 0.2*u.um) + SpectralRegion(0.3*u.um, 0.4*u.um) +
...       SpectralRegion(0.45*u.um, 0.6*u.um) + SpectralRegion(0.8*u.um, 0.9*u.um) +
...       SpectralRegion(1.0*u.um, 1.2*u.um) + SpectralRegion(1.3*u.um, 1.5*u.um))

>>> sr_inverted = sr.invert(0.05*u.um, 3*u.um)
>>> sr_inverted
Spectral Region, 7 sub-regions:
(0.05 um, 0.15 um)   (0.2 um, 0.3 um)     (0.4 um, 0.45 um)
(0.6 um, 0.8 um)     (0.9 um, 1.0 um)     (1.2 um, 1.3 um)
(1.5 um, 3.0 um)

Region Extraction

Given a SpectralRegion, one can extract a sub-spectrum from a Spectrum1D object. If the SpectralRegion has multiple sub-regions then by default a list of Spectrum1D objects will be returned. If the return_single_spectrum argument is set to True, the resulting spectra will be concatenated together into a single Spectrum1D object instead.

An example of a single sub-region SpectralRegion:

>>> from astropy import units as u
>>> import numpy as np
>>> from specutils import Spectrum1D, SpectralRegion
>>> from specutils.manipulation import extract_region

>>> region = SpectralRegion(8*u.nm, 22*u.nm)
>>> spectrum = Spectrum1D(spectral_axis=np.arange(1, 50) * u.nm, flux=np.random.sample(49)*u.Jy)
>>> sub_spectrum = extract_region(spectrum, region)
>>> sub_spectrum.spectral_axis
<SpectralAxis [ 8.,  9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21.,
   22.] nm>

Extraction also correctly interprets different kinds of spectral region units as would be expected:

>>> from astropy import units as u
>>> import numpy as np
>>> from specutils import Spectrum1D, SpectralRegion
>>> from specutils.manipulation import extract_region

>>> spectrum = Spectrum1D(spectral_axis=np.arange(1, 50) * u.nm, flux=np.random.sample(49)*u.Jy)
>>> region_angstroms = SpectralRegion(80*u.AA, 220*u.AA)
>>> sub_spectrum = extract_region(spectrum, region_angstroms)
>>> sub_spectrum.spectral_axis
<SpectralAxis [ 8.,  9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21.,
   22.] nm>
>>> region_pixels = SpectralRegion(7.5*u.pixel, 21.5*u.pixel)
>>> sub_spectrum = extract_region(spectrum, region_pixels)
>>> sub_spectrum.spectral_axis
<SpectralAxis [ 8.,  9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21.,
   22.] nm>

An example of a multiple sub-region SpectralRegion:

>>> from astropy import units as u
>>> import numpy as np
>>> from specutils import Spectrum1D, SpectralRegion
>>> from specutils.manipulation import extract_region

>>> region = SpectralRegion([(8*u.nm, 22*u.nm), (34*u.nm, 40*u.nm)])
>>> spectrum = Spectrum1D(spectral_axis=np.arange(1, 50) * u.nm, flux=np.random.sample(49)*u.Jy)
>>> sub_spectra = extract_region(spectrum, region)
>>> sub_spectra[0].spectral_axis
<SpectralAxis [ 8.,  9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21.,
   22.] nm>
>>> sub_spectra[1].spectral_axis
<SpectralAxis [34., 35., 36., 37., 38., 39., 40.] nm>

Multiple sub-regions can also be returned as a single concatenated spectrum:

>>> sub_spectrum = extract_region(spectrum, region, return_single_spectrum=True)
>>> sub_spectrum.spectral_axis
<SpectralAxis [ 8.,  9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21.,
    22., 34., 35., 36., 37., 38., 39., 40.] nm>

The bounding region that includes all data, including the ones that lie in between disjointed spectral regions, can be extracted with specutils.manipulation.extract_bounding_spectral_region:

>>> from astropy import units as u
>>> import numpy as np
>>> from specutils import Spectrum1D, SpectralRegion
>>> from specutils.manipulation import extract_bounding_spectral_region

>>> spectrum = Spectrum1D(spectral_axis=np.arange(1, 50) * u.nm,
...                       flux=np.random.default_rng(12345).random(49)*u.Jy)
>>> region = SpectralRegion([(8*u.nm, 12*u.nm), (24*u.nm, 30*u.nm)])
>>> sub_spectrum = extract_bounding_spectral_region(spectrum, region)
>>> sub_spectrum.spectral_axis
<SpectralAxis [ 8.,  9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21.,
    22., 23., 24., 25., 26., 27., 28., 29., 30.] nm>

spectral_slab is basically an alternate entry point for extract_region. Notice the slightly different way to input the spectral axis range to be extracted. This function’s purpose is to facilitate migration of spectral_cube functionality into specutils:

>>> from astropy import units as u
>>> import numpy as np
>>> from specutils import Spectrum1D, SpectralRegion
>>> from specutils.manipulation import spectral_slab

>>> spectrum = Spectrum1D(spectral_axis=np.arange(1, 50) * u.nm,
...                       flux=np.random.default_rng(12345).random(49)*u.Jy)
>>> sub_spectrum = spectral_slab(spectrum, 8*u.nm, 20*u.nm)
>>> sub_spectrum.spectral_axis
<SpectralAxis [ 8.,  9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20.] nm>

Line List Compatibility

SpectralRegion objects can also be created from the QTable object returned from the line finding functions:

>>> from astropy import units as u
>>> import numpy as np
>>> from specutils import Spectrum1D, SpectralRegion
>>> from astropy.modeling.models import Gaussian1D
>>> from specutils.fitting import find_lines_derivative

>>> g1 = Gaussian1D(1, 4.6, 0.2)
>>> g2 = Gaussian1D(2.5, 5.5, 0.1)
>>> g3 = Gaussian1D(-1.7, 8.2, 0.1)

>>> x = np.linspace(0, 10, 200)
>>> y = g1(x) + g2(x) + g3(x)
>>> spectrum = Spectrum1D(flux=y * u.Jy, spectral_axis=x * u.um)

>>> lines = find_lines_derivative(spectrum, flux_threshold=0.01)
>>> spec_reg = SpectralRegion.from_line_list(lines)
>>> spec_reg
Spectral Region, 3 sub-regions:
  (4.072864321608041 um, 5.072864321608041 um)
  (4.977386934673367 um, 5.977386934673367 um)
  (7.690954773869347 um, 8.690954773869347 um)

This can be fed into the exclude_regions argument of the fit_generic_continuum or fit_continuum functions to avoid fitting regions that contain line features. Note that, by default, this uses pythonic slicing, i.e., spectral values greater than or equal to the lower bound and less than the upper bound of the region will be excluded from the fit. For convenience in some cases, the exclude_region_upper_bounds keyword can be set to True to exlude spectral values less than or equal to the upper bound instead.

Conversely, users can also invert the spectral region

>>> inv_spec_reg = spec_reg.invert(spectrum.spectral_axis[0], spectrum.spectral_axis[-1])
>>> inv_spec_reg
Spectral Region, 3 sub-regions:
  (0.0 um, 4.072864321608041 um)
  (5.977386934673367 um, 7.690954773869347 um)
  (8.690954773869347 um, 10.0 um)

and use that result as the exclude_regions argument in the fit_lines function in order to avoid attempting to fit any of the continuum region.

Reading and Writing

SpectralRegion objects can be written to ecsv files, which uses the QTable write machinery:

>>> spec_reg.write("spectral_region.ecsv")

This overwrites existing files by default if a duplicate filename is input. The resulting files can be read back in to create a new SpectralRegion using the read class method:

>>> spec_reg = SpectralRegion.read("spectral_region.ecsv")

The QTable created to write out the SpectralRegion to file can also be accessed directly with the as_table method, and a SpectralRegion can be created directly from a QTable with the appropriate columns (minimally lower_bound and upper_bound) using the from_qtable class method.

>>> spec_reg_table = spec_reg.as_table()
>>> spec_reg_2 = SpectralRegion.from_qtable(spec_reg_table)

Reference/API

Classes

Conf()	Configuration parameters for specutils.
SpectralRegion(*args)	A SpectralRegion is a container class for regions (intervals) along a spectral coordinate.
SpectrumList([iterable])	A list that is used to hold a list of Spectrum1D objects