astropy
Spectrum1D axis order discussion
This is intended to hold space to continue the discussion about Spectrum1D axis ordering from #822.
To summarize the discussion thus far from #822:
I said:
I guess here's what I see as the science-centric perspective thinking about a spectral-spatial-spatial to have a concrete way to think about it. Here's what I see as the result of
spectrum[..., spec_idx]in the two cases hgere:
moveaxis: this is the right thing to reproduce the behavior oforiginal_cube[spec_idx, :, :]and have that be "the spatial plane", which I think is what the science user generally expectsswapaxis: this givesoriginal_cube[spec_idx, :, :].T, which is not as intuitive. And it's much worse with like a spatial-spectral-spatial-time or similar >3D datasetSo based on that it feels to me like
moveaxisis the right answer. But I think you might be right that this is not consistent with WCS? But I could be wrong about that. I am constantly confused about transposes in WCS land. My approahc is usually transpose-until-its-right :shrug: ...
and later
Given how tricky this is I'm going to ask for a few more perspectives: @Cadair @PatrickOgle @keflavich @havok2063 .
The TL;DR version is: given that specutils requires the spectral axis be the last axis, what's the more logical way to re-order N-D cubes: "swapping" the last and the spectral axis, or "moving" the spectral axis to the end (which preserves the order of all the other axes).
The second question is: is whatever the answer is above ok even if the WCS expects it to be transposed or something?
I think everyone then agreed the moveaxis solution was the better one to present to the users. By that I really mean what @PatrickOgle said, that the index order should be preserved from whatever the "in-file" representation is (i.e. astropy.io.fits for fits cubes), but with the last index being the spectral axis. (This doesn't necessarily mean literally using moveaxis - in my post that was just the easiest way at the time to explain the case.)
Also, as a quick note (inspired partly by @keflavich comment in #822), I think we should lean on test-driven design for this one. By that I mean write a few tests that behave the way we want, agree on those, and only after that does anyone work on implementation.
This is somewhat similar to the issue that when reading a 2D image from a FITS file, the Numpy axis order is reversed from the WCS. I would personally be in favour of staying consistent with the 'numpy axes are reversed' convention to be consistent with other tools. Otherwise aren't we ending up with an axis order of y,x,spectral? If so why aren't x and y swapped too? One benefit of keeping the original spectral,y,x ordering is that selecting a flux/image from the spectrum is simply spectrum[ispec]. So personally I think neither moveaxis or swapaxis should be used! 😂
Whatever decision is made should be justified in the specutils documentation - I personally do not understand at all the benefit of having the spectral axis be last. I could understand re-ordering it so that it is always first as done in SpectralCube but I can't see why last would be useful or intuitive.
swapaxis is definitely not the right answer though, as the order will be really whacky for a 4D WCS! (where now it will no longer be simply 'reversed')
I would personally be in favour of staying consistent with the 'numpy axes are reversed' convention
I would go as far as to say that you shouldn't be inheriting from NDCube if you aren't doing this. It's a fundamental design decision of APE 14 that the world coordinates are cartesian ordered and the arrays aren't. We used this as a basic building block for all the NDCube 2 design. This is a hill I will now die on, despite me being against it when the decision was made.
Just to add to my comment about 4D WCS, a dataset in Numpy order (spec, time, dec, ra) will now be (ra, time, dec, spec) which seems very confusing.
