# 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 it's canonical source:

f = fits.open('https://dr14.sdss.org/optical/spectrum/view/data/format=fits/spec=lite?plateid=1323&mjd=52797&fiberid=12')  # doctest: +IGNORE_OUTPUT +REMOTE_DATA
# The spectrum is in the second HDU of this file.
specdata = f[1].data # doctest: +REMOTE_DATA
f.close() # doctest: +REMOTE_DATA

# Then we re-format this dataset into astropy quantities, and create a
# `~specutils.Spectrum1D` object:

from specutils import Spectrum1D
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 = Spectrum1D(spectral_axis=lamb, flux=flux) # doctest: +REMOTE_DATA

# And we plot it:

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