I disagree with that last point.  My main complaint about WCS, and it is a persistent one, is that it does not follow the 'principle of least surprise'.  I expect to be able to specify my units in a FITS header and get those exact units out.

Your point about implicitly calling `with_spectral_unit` is valid, though - I don't think it should ever be possible to instantiate a cube with something other than its default unit type, e.g., `SpectralCube(..., spectral_unit=u.m)` should not be possible if the cube is natively in `u.m/u.s`.  What I would like is for `SpectralCube` to be a little smarter than `WCS` and recognize that a header in `u.km/u.s` is different from one with `u.m/u.s`.

Thinking about it more, I really thing this means that the `io.fits` module could return `SpectralCube(...)_new_cube_with(spectral_unit=header['CUNITn'])` - then we don't need anything in `__init__` and we just need smarter `io` functions.  Is that a reasonable compromise?