Working with Spectral Cubes¶
Spectral cubes can be read directly with Spectrum1D
.
A specific example of this is demonstrated here. In addition to the functions
demonstrated below, as a subclass of NDCube,
Spectrum1D
also inherits useful methods for e.g. cropping
based on combinations of world and spectral coordinates. Most of the
functionality inherited from NDCube
requires initializing the
Spectrum1D
object with a WCS describing the coordinates for all axes of
the data.
Note that the workflow described here is for spectral cubes that are rectified
such that one of the axes is entirely spectral and all the spaxels have the same
spectral_axis
values (i.e., case 2 in Overview of How Specutils Represents Spectra).
For less-rectified cubes, pre-processing steps (not addressed by specutils at the
time of this writing) will be necessary to rectify the cubes into that form.
Note, though, that such cubes can be stored in specutils data structures (cases
3 and 4 in Overview of How Specutils Represents Spectra), which support some of
these behaviors, even though the fulkl set of tools do not yet apply.
Loading a cube¶
We’ll use a MaNGA cube
for our example, and load the data from the
repository directly into a new Spectrum1D
object:
>>> from astropy.utils.data import download_file
>>> from specutils.spectra import Spectrum1D
>>> filename = "https://stsci.box.com/shared/static/28a88k1qfipo4yxc4p4d40v4axtlal8y.fits"
>>> file = download_file(filename, cache=True)
>>> sc = Spectrum1D.read(file, format='MaNGA cube')
The cube has 74x74 spaxels with 4563 spectral axis points in each one:
>>> sc.shape
(74, 74, 4563)
Print the contents of 3 spectral axis points in a 3x3 spaxel array:
>>> sc[30:33,30:33,2000:2003]
<Spectrum1D(flux=[[[0.4892023205757141 ... 0.5994223356246948]]] 1e-17 erg / (Angstrom s spaxel cm2) (shape=(3, 3, 3), mean=0.54165 1e-17 erg / (Angstrom s spaxel cm2)); spectral_axis=<SpectralAxis
(observer to target:
radial_velocity=0.0 km / s
redshift=0.0)
[5.73984286e-07 5.74116466e-07 5.74248676e-07] m> (length=3); uncertainty=InverseVariance)>
Spectral slab extraction¶
The spectral_slab
function can be used to extract
spectral regions from the cube.
>>> import astropy.units as u
>>> from specutils.manipulation import spectral_slab
>>> ss = spectral_slab(sc, 5000.*u.AA, 5003.*u.AA)
>>> ss.shape
(74, 74, 3)
>>> ss[30:33,30:33,::]
<Spectrum1D(flux=[[[0.6103081107139587 ... 0.936118483543396]]] 1e-17 erg / (Angstrom s spaxel cm2) (shape=(3, 3, 3), mean=0.83004 1e-17 erg / (Angstrom s spaxel cm2)); spectral_axis=<SpectralAxis
(observer to target:
radial_velocity=0.0 km / s
redshift=0.0)
[5.00034537e-07 5.00149688e-07 5.00264865e-07] m> (length=3); uncertainty=InverseVariance)>
Spectral Bounding Region¶
The extract_bounding_spectral_region
function can be used to
extract the bounding region that encompases a set of disjoint SpectralRegion
instances, or a composite instance of SpectralRegion
that contains
disjoint sub-regions.
>>> from specutils import SpectralRegion
>>> from specutils.manipulation import extract_bounding_spectral_region
>>> composite_region = SpectralRegion([(5000*u.AA, 5002*u.AA), (5006*u.AA, 5008.*u.AA)])
>>> sub_spectrum = extract_bounding_spectral_region(sc, composite_region)
>>> sub_spectrum.spectral_axis
<SpectralAxis
(observer to target:
radial_velocity=0.0 km / s
redshift=0.0)
[5.00034537e-07, 5.00149688e-07, 5.00264865e-07, 5.00380068e-07,
5.00495298e-07, 5.00610555e-07, 5.00725838e-07] m>
Moments¶
The moment
function can be used to compute moments of any order
along one of the cube’s axes. By default, axis=-1
, which computes moments
along the spectral axis (remember that the spectral axis is always last in a
Spectrum1D
).
>>> from specutils.analysis import moment
>>> m = moment(sc, order=1)
>>> m.shape
(74, 74)
>>> m[30:33,30:33]
<Quantity [[6.97933331e-07, 6.98926463e-07, 7.00540974e-07],
[6.98959625e-07, 7.00280655e-07, 7.03511823e-07],
[7.00740294e-07, 7.04527986e-07, 7.08245958e-07]] m>
Use Case¶
Example of computing moment maps for specific wavelength ranges in a
cube, using spectral_slab
and
moment
.
import numpy as np
import matplotlib.pyplot as plt
import astropy.units as u
from astropy.utils.data import download_file
from specutils import Spectrum1D, SpectralRegion
from specutils.analysis import moment
from specutils.manipulation import spectral_slab
filename = "https://stsci.box.com/shared/static/28a88k1qfipo4yxc4p4d40v4axtlal8y.fits"
fn = download_file(filename, cache=True)
spec1d = Spectrum1D.read(fn)
# Extract H-alpha sub-cube for moment maps using spectral_slab
subspec = spectral_slab(spec1d, 6745.*u.AA, 6765*u.AA)
ha_wave = subspec.spectral_axis
# Extract wider sub-cube covering H-alpha and [N II] using spectral_slab
subspec_wide = spectral_slab(spec1d, 6705.*u.AA, 6805*u.AA)
ha_wave_wide= subspec_wide.spectral_axis
# Convert flux density to microJy and correct negative flux offset for
# this particular dataset
ha_flux = (np.sum(subspec.flux.value, axis=(0,1)) + 0.0093) * 1.0E-6*u.Jy
ha_flux_wide = (np.sum(subspec_wide.flux.value, axis=(0,1)) + 0.0093) * 1.0E-6*u.Jy
# Compute moment maps for H-alpha line
moment0_halpha = moment(subspec, order=0)
moment1_halpha = moment(subspec, order=1)
# Convert moment1 from AA to velocity
# H-alpha is redshifted to 6755 AA for this galaxy
print(moment1_halpha[40,40])
vel_map = 3.0E5 * (moment1_halpha.value - 6.755E-7) / 6.755E-7
# Plot results in 3 panels (subspec_wide, H-alpha line flux, H-alpha velocity map)
f,(ax1,ax2,ax3) = plt.subplots(1, 3, figsize=(15, 5))
ax1.plot(ha_wave_wide, (ha_flux_wide)*1000.)
ax1.set_xlabel('Angstrom', fontsize=14)
ax1.set_ylabel('uJy', fontsize=14)
ax1.tick_params(axis="both", which='major', labelsize=14, length=8, width=2, direction='in', top=True, right=True)
ax2.imshow(moment0_halpha.value, origin='lower')
ax2.set_title('moment = 0')
ax2.set_xlabel('x pixels', fontsize=14)
ax3.imshow(vel_map, vmin=-100., vmax=100., cmap='rainbow', origin='lower')
ax3.set_title('moment = 1')
ax3.set_xlabel('x pixels', fontsize=14)
(Source code
, png
, hires.png
, pdf
)
