WIP: Keep full 3D WCS objects for LDOs

[keflavich]

Related to #111, #367, #392, #394:
When slicing or projecting a cube, there is generally still a 3rd dimension,
even if it is not reflected in the data. Currently, if you attempt to do
anything with the WCS of a position-velocity slice, the WCS crashes, and we
want to avoid that.

There are big open questions though:

• for projections, what coordinate applies to operations like .max?
• for projections, can we assume the central pixel (the 'mean' along that
  axis) is the appropriate coordinate?
• What should CDELT be?

in slices, I think all of these questions are easily answered.

We will need to provide convenience functions (WCS wrappers) to get the
coordinates in projections and slices, because most users will never care about
the sliced-over axes, or they will care but only want to see them as metadata.

cc @e-koch, who has contributed significantly to this discussion already.

[keflavich]

Tests are pretty close to passable, but they fail for the aplpy case because now we're writing out HDUs with 3 dimensions for a 2 dimensional object. We just need to tell that to aplpy in the vis calls...

[keflavich]

We really need this, unfortunately.

Here's an example case I had to do to get WCSes attached properly:

for fn in cubefns:
    cube = SpectralCube.read(fn).subcube_from_regions([reg])
    print(cube)
    avg = cube.mean(axis=(1,2))
    hdu = avg.hdu
    y,x = cube.wcs.celestial.world_to_pixel(reg.center)

    slc = getslice()[int(floor(x)):int(ceil(x)), int(floor(y)):int(ceil(y)), :]
    ww = wcs_utils.slice_wcs(cube.wcs, slc, numpy_order=False)
    hdu.header.update(ww.to_header())
    hdu.data = hdu.data[None,None,:]