Specutils

A SpectrumCollection is a container that allows a set of spectra to behave like a single Spectrum object. This enables you to perform analysis operations over entire sets of data simultaneously.

Key Requirements & Behaviors:

• Shape Consistency: All spectra within a collection must have the same shape.
• Dispersion Solutions: The collection does not make assumptions about dispersion solutions. Users should ensure spectra are resampled to a common grid or are aware that they may not share the same dispersion solution.
• Data Storage: Unlike a list of objects, a SpectrumCollection stores data as single multi-dimensional arrays (e.g., a 2D flux array). It constructs individual Spectrum objects on-the-fly only when you index the collection.
• Immutability via Indexing: Because indexing creates a new Spectrum object from a slice, performing operations on a spectrum retrieved from the collection (e.g., spec_coll[0].flux *= 2) will not modify the original SpectrumCollection object.

import numpy as np
import astropy.units as u
from astropy.nddata import StdDevUncertainty
from specutils import SpectrumCollection
from specutils.utils.wcs_utils import gwcs_from_array

# Setup data for 5 spectra, each with 10 data points
flux = u.Quantity(np.random.sample((5, 10)), unit='Jy')
spectral_axis = u.Quantity(np.arange(50).reshape((5, 10)), unit='AA')
wcs = np.array([gwcs_from_array(x, x.shape) for x in spectral_axis])
uncertainty = StdDevUncertainty(np.random.sample((5, 10)), unit='Jy')
mask = np.ones((5, 10)).astype(bool)
meta = [{'test': 5, 'info': [1, 2, 3]} for i in range(5)]

# Initialize the collection
spec_coll = SpectrumCollection(
    flux=flux, 
    spectral_axis=spectral_axis, 
    wcs=wcs, 
    uncertainty=uncertainty, 
    mask=mask, 
    meta=meta
)

print(spec_coll.shape)          # Output: (5,)
print(spec_coll.flux.unit)     # Output: Unit("Jy")
print(spec_coll.spectral_axis.shape) # Output: (5, 10)

A Spectrum can represent multi-dimensional data sets (e.g., an (n_spectra, n_pix) array) where multiple flux arrays share a single common spectral_axis.

• Automatic Axis Detection: Specutils attempts to match the shape of the flux and uncertainty arrays to the spectral_axis.
• Manual Axis Specification: If multiple axes have the same length as the spectral axis, you must specify which one is the spectral axis using the spectral_axis_index argument during initialization.
• Reordering Axes: Use the move_spectral_axis argument to force the spectral axis to a specific dimension. Setting move_spectral_axis=-1 or move_spectral_axis='last' reproduces the pre-v2.0 behavior where the spectral axis was always the last axis. This swaps the relevant axes in flux, mask, and uncertainty and updates the WCS accordingly.

Note: If each flux array is related to a different spectral axis, use SpectrumCollection instead.

from specutils import Spectrum
import numpy as np
import astropy.units as u

# Example: (5, 10) flux array with a spectral axis of length 10
spec = Spectrum(spectral_axis=np.arange(5000, 5010) * u.AA,
                flux=np.random.default_rng(12345).random((5, 10)) * u.Jy)

# Slicing a single spectrum from the collection
spec_slice = spec[0]
print(spec_slice.spectral_axis)
print(spec_slice.flux)

A SpectralRegion defines one or more sub-regions in spectral space. They are independent of any specific Spectrum object.

Creation

You can define a region in several ways:

• Single region: Pass two astropy.units.Quantity objects (e.g., SpectralRegion(0.45*u.um, 0.6*u.um)).
• Multiple sub-regions: Pass a list of 2-tuples (e.g., SpectralRegion([(0.45*u.um, 0.6*u.um), (0.8*u.um, 0.9*u.um)])).
• Units: Supports any valid spectral unit or u.pixel for direct indexing.

Operations

• Combining: Use the + operator to combine regions or += for in-place addition.
• Slicing: Use integer indexing or range slicing (e.g., sr[0] or sr[3:5]).
• Deletion: Use del sr[index] to remove a sub-region.
• Iteration: You can iterate over a SpectralRegion to access its individual sub-regions.
• Inversion: Use .invert(lower, upper) to return a region representing the space between existing sub-regions within a specified bound. This is useful for finding baseline/noise regions around lines.

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

# Single region
sr = SpectralRegion(0.45*u.um, 0.6*u.um)

# Multiple sub-regions
sr_two = SpectralRegion([(0.45*u.um, 0.6*u.um), (0.8*u.um, 0.9*u.um)])

# Combining
sr_combined = SpectralRegion(0.45*u.um, 0.6*u.um) + SpectralRegion(0.8*u.um, 0.9*u.um)

# Inverting (finding gaps between regions)
sr_inverted = sr_two.invert(0.05*u.um, 3*u.um)

Multiple SpectralRegion objects can be combined into a CompoundSpectralRegion using bitwise operators:

• | (OR): Combines regions.
• & (AND): Intersection of regions.
• ^ (XOR): Symmetric difference.

Once combined, you can:

1. Check if a specific wavelength is inside the region using .contains(wavelength).
2. Convert the compound region to a boolean mask for a Spectrum using .to_mask(spectrum).
3. Flatten the compound region into a set of simple sub-regions using .to_subregions().

import astropy.units as u
from specutils import Spectrum, SpectralRegion

reg1 = SpectralRegion(5*u.um, 7*u.um)
reg2 = SpectralRegion([(4*u.um, 5.5*u.um), (6.5*u.um, 8*u.um)])

# Combine using OR
compound = reg1 | reg2

# Check containment
print(compound.contains(4.5*u.um))  # True

# Convert to mask
spectrum = Spectrum(spectral_axis=np.arange(1, 10) * u.um, flux=np.ones(9)*u.Jy)
mask = compound.to_mask(spectrum)

# Convert to simple sub-regions
subregions = compound.to_subregions()

Specutils maintains a WCS (World Coordinate System) object for every Spectrum. If not explicitly provided, specutils creates an internal GWCS object. Currently, specutils supports two WCS formats:

1. FITS WCS (via astropy.wcs)
2. GWCS

To provide a FITS-style WCS, pass an astropy.wcs.WCS object to the wcs argument during Spectrum initialization. You can then use spec.wcs.pixel_to_world() to convert pixel coordinates to world coordinates.

from specutils.spectra import Spectrum
import astropy.wcs as fitswcs
import astropy.units as u
import numpy as np

my_wcs = fitswcs.WCS(header={
    'CDELT1': 1, 'CRVAL1': 6562.8, 'CUNIT1': 'Angstrom', 'CTYPE1': 'WAVE',
    'RESTFRQ': 1400000000, 'CRPIX1': 25
})
spec = Spectrum(flux=[5,6,7] * u.Jy, wcs=my_wcs)

# Access spectral axis and convert pixels to world
print(spec.spectral_axis)
print(spec.wcs.pixel_to_world(np.arange(3)))

Specutils uses different data structures depending on the dimensionality and shape of the spectroscopic data. The choice of object depends on whether the flux is 1D, multidimensional, or part of a collection with varying spectral axes.

1. Single or Multidimensional Spectra (Spectrum)

Use the specutils.Spectrum class for:

• 1D Spectra: A single flux of length n with a matched spectral axis (array or WCS).
• Multidimensional Fluxes: A set of fluxes with shape n x m (x ...) where the spectral axis is strictly of length n. Spectrum subclasses NDCube, providing utilities for multidimensional data.

2. Uniform Collections (SpectrumCollection)

Use specutils.SpectrumCollection when you have a set of fluxes of shape n x m (x ...) and a corresponding set of spectral axes of the same shape.

• Performance: This is more performant than a list of individual objects because it stores data in unified arrays.
• Behavior: It does not store individual Spectrum objects. Instead, it constructs a Spectrum object on the fly when indexed. Note: Operations performed on a Spectrum returned from a SpectrumCollection do not modify the original collection.

3. Ragged or Heterogeneous Collections

For collections where spectra have different lengths or different spectral axes (ragged data):

• Standard Approach: Use a standard Python list or a numpy object-array containing individual specutils.Spectrum objects. Iterate over the list to process them.
• SpectrumList: A subclass of list integrated with the Astropy IO registry. It is primarily used by data loaders to return multiple heterogeneous spectra. For most user workflows, a standard list of Spectrum objects is sufficient.