Improve WCS handling when slicing/extracting regions

As of https://github.com/astropy/specutils/pull/918, slicing the spectral axis of a `Spectrum1D` will drop the existing WCS and recreate a fresh lookup table GWCS from the new spectral axis. It seems like it would be better to propagate any spatial components of the existing WCS for a multi-dimensional `Spectrum1D`. Additionally, for spatial slicing it would be nice to properly handle any spatial components of the WCS so that the resulting `Spectrum1D` has a WCS with as much information as the original. I'm not sure, but NDCube may handle this better than we currently do here for at least spatial dimensions - I'll have to check. It may be that we're overriding something useful.

Start by reviewing specutils PR 918 and the Spectrum1D slicing and spatial-extraction paths described in this issue. Compare their WCS behavior with NDCube, then establish coverage showing that spectral and spatial slicing preserve as much of the original WCS information as appropriate.

Written by the indexing model from the issue text.