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 :class:`~specutils.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
# :class:`~specutils.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  # doctest: +IGNORE_OUTPUT

f, ax = plt.subplots()  # doctest: +IGNORE_OUTPUT
ax.step(spec.spectral_axis, spec.flux) # doctest: +IGNORE_OUTPUT +REMOTE_DATA