Collapsing Cubes

Currently, the specutils library does not support loading three dimensional spectral cubes and manipulating their non-spectral axes, so parsing a spectral cube to a Spectrum1D object requires some programmatic manipulation using the spectral-cube library. Below is an example of loading a spectral cube object, collapsing the cube, and creating a Spectrum1D object.

Loading a cube

We’ll use the example cube data found in the spectral-cube package, and load the data from the repository directly into a new SpectralCube object:

>>> from astropy.utils.data import download_file
>>> from spectral_cube import SpectralCube
>>> sc = SpectralCube.read(download_file("https://github.com/radio-astro-tools/spectral-cube/raw/master/spectral_cube/tests/data/example_cube.fits"), format='fits')
>>> sc
SpectralCube with shape=(7, 4, 3) and unit=Jy / beam:
n_x:      3  type_x: RA---ARC  unit_x: deg    range:    52.231466 deg:   52.231544 deg
n_y:      4  type_y: DEC--ARC  unit_y: deg    range:    31.243639 deg:   31.243739 deg
n_s:      7  type_s: VRAD      unit_s: m / s  range:    14322.821 m / s:   14944.909 m / s

We see that the example cube is small, with a shape of (7, 3, 4); two spatial axes RA and DEC, and one spectral axis that’s in velocity units.

Extracting the spectral data

Now that we have a data cube, our mission is to collapse this cube along the two spatial axes to get the representative spectral data. We can then parse the result into a Spectrum1D object. SpectralCube can take and apply arbitrary functions over axes of the data, in this case we can apply a simple median function over the spatial axes:

>>> import numpy as np
>>> sc_spec = sc.apply_numpy_function(np.mean, axis=(1,2), projection=True)
>>> sc_spec.shape
(7,)

The projection keyword returns a lower dimesional object that retains the WCS information that we’ll need to use when constructing the Spectrum1D object.

Creating a spectrum object

With our data parsed and extracted, let’s go ahead and create our Spectrum1D object:

>>> from specutils import Spectrum1D
>>> spec = Spectrum1D(flux=sc_spec.data[:] * sc_spec.unit, wcs=sc_spec.wcs)

Let’s plot and see what this resulting spectrum looks like:

>>> from matplotlib import pyplot as plt
>>> from astropy.visualization import quantity_support
>>> quantity_support()  # for getting units on the axes below  

>>> f, ax = plt.subplots()  
>>> ax.step(spec.spectral_axis, spec.flux)