Collapsed shape from cube is very confusing

[pllim]

Consider this use case:

>>> import numpy as np
>>> from astropy import units as u
>>> from astropy.wcs import WCS
>>> from specutils import Spectrum1D
>>> wcs = WCS(naxis=3)
>>> wcs.wcs.ctype = 'RA---TAN', 'DEC--TAN', 'FREQ'
>>> wcs.wcs.set()
>>> sp = Spectrum1D(flux=np.ones((3, 4, 5)) * u.nJy, wcs=wcs)  # spectral, y, x
WARNING:specutils.spectra.spectrum1d:Input WCS indicates that the
spectral axis is not last. Reshaping arrays to put spectral axis last.
>>> sp.shape  # x, y, spectral
(5, 4, 3)
>>> collapsed_spec = sp.collapse('sum', axis=-1)
>>> collapsed_spec.shape
(5, 4)  # x, y
>>> collapsed_spec.T.shape  # y, x
(4, 5)

Expectation: Collapsed shape should be (y, x).

Reality: It becomes (x, y). As a result, I need to transpose the collapsed result to get the spatial image in the correct shape.

xref spacetelescope/jdaviz#1006

cc @rosteen @javerbukh