# We begin with some basic  imports:
#
from astropy.io import fits
from astropy import units as u
import numpy as np
from matplotlib import pyplot as plt
from astropy.visualization import quantity_support
quantity_support()  # for getting units on the axes below  # doctest: +IGNORE_OUTPUT
#
# Now we load the dataset from its canonical source:
#
filename = 'https://data.sdss.org/sas/dr16/sdss/spectro/redux/26/spectra/1323/spec-1323-52797-0012.fits'
with fits.open(filename) as f:  # doctest: +IGNORE_OUTPUT +REMOTE_DATA
    specdata = f[1].data  # doctest: +REMOTE_DATA
#
# Then we re-format this dataset into astropy quantities, and create a
# `~specutils.Spectrum` object:
#
from specutils import Spectrum
lamb = 10**specdata['loglam'] * u.AA # doctest: +REMOTE_DATA
flux = specdata['flux'] * 10**-17 * u.Unit('erg cm-2 s-1 AA-1') # doctest: +REMOTE_DATA
spec = Spectrum(spectral_axis=lamb, flux=flux) # doctest: +REMOTE_DATA
#
# And we plot it:
#
f, ax = plt.subplots()  # doctest: +IGNORE_OUTPUT +REMOTE_DATA
ax.step(spec.spectral_axis, spec.flux) # doctest: +IGNORE_OUTPUT +REMOTE_DATA