from specutils.fitting import fit_generic_continuum
cont_norm_spec = spec / fit_generic_continuum(spec)(spec.spectral_axis) # doctest: +REMOTE_DATA
lines = plt.step(cont_norm_spec.wavelength, cont_norm_spec.flux) # doctest: +REMOTE_DATA
plt.xlim(654*u.nm, 660*u.nm)  # doctest: +IGNORE_OUTPUT +REMOTE_DATA