Allow spectral axis to be anywhere, instead of forcing it to be last

[rosteen]

Currently under heavy development (in other words, super broken 😆 ).

Edit to provide context in the initial comment:

The motivation for this is that the insistence on the spectral axis being last hinders integration with other packages that make different (or no) assumptions on the location of the spectral axis. The two biggest examples, currently, are that spectral-cube assumes that the spectral axis is first, and they would be interested in using specutils objects as the basis of their code but are hindered by our assumption. We don't need to switch to that assumption, I think being flexible here is the right call. jdaviz also had to do a lot of coding around Glue assumptions about where the spectral axis is as well, since they differ from ours.

I have also been concerned for a while that specutils putting the flux axis into a different order than a user would get out from astropy.io.fits is potentially confusing and/or a barrier to uptake of specutils. In my mind at least, this combined with the above outweighs the convenience of numpy broadcasting properly automatically if you index solely on the spectral axis when it is last for a multi-dimensional Spectrum1D, which if I'm recalling correctly was the motivation for forcing the spectral axis to be last.

I do think that it would be useful to provide a spectral_axis_location keyword when initializing a Spectrum1D so that users could either get the legacy behavior (force it to be last) or do the opposite (force first) if desired, which I haven't yet implemented here. I could get that in before merging this, but I also think that I'm going to end up closing this PR and reopening it to merge into a new v2.0-dev branch to keep these breaking changes out of main for now, so I could also open a followup PR to that branch after this is merged.

[codecov]

Codecov Report

Merging #999 (e38c008) into main (6eb7f96) will increase coverage by 0.37%.
The diff coverage is 83.89%.

@@            Coverage Diff             @@
##             main     #999      +/-   ##
==========================================
+ Coverage   70.09%   70.47%   +0.37%     
==========================================
  Files          64       64              
  Lines        4344     4396      +52     
==========================================
+ Hits         3045     3098      +53     
+ Misses       1299     1298       -1     

Impacted Files	Coverage Δ
specutils/fitting/fitmodels.py	68.30% <ø> (ø)
specutils/io/default_loaders/sdss.py	47.19% <ø> (ø)
specutils/manipulation/resample.py	94.94% <ø> (ø)
specutils/spectra/spectrum_mixin.py	67.71% <75.00%> (+0.51%)	⬆️
specutils/spectra/spectrum1d.py	82.53% <80.68%> (+1.60%)	⬆️
specutils/analysis/moment.py	94.44% <92.85%> (-1.71%)	⬇️
specutils/io/default_loaders/wcs_fits.py	35.22% <100.00%> (+0.26%)	⬆️
specutils/manipulation/extract_spectral_region.py	94.17% <100.00%> (+3.90%)	⬆️

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

[rosteen]

I finally have all tests passing for this, so I'm moving it to ready for review. I still need to update some of the narrative docs, and I'm sure there are a couple things that I haven't considered/weren't caught by the tests, so extra eyes are welcome.

[keflavich]

I think this looks good, but there's also a lot of change to uncertainty things - is that unrelated cleanup being smuggled in? Fine by me if it is, but I didn't review that carefully.

The addition to the narrative docs will be needed, so approval waits on that. Spectrum1D.spectral_axis_index is the key variable, right?

[rosteen]

I think this looks good, but there's also a lot of change to uncertainty things - is that unrelated cleanup being smuggled in? Fine by me if it is, but I didn't review that carefully.

The uncertainty stuff was mostly deletions of the handling for making sure the uncertainty also got reordered when moving the spectral axis to last. The rest of the diffs with uncertainty in them are just adding another argument after uncertainty in function calls.

[rosteen]

I pushed a small update to the docs to change from specifying that the spectral axis is always last to saying that it can be anywhere.

[pllim]

If we do this, how do we know how to access the wcs.pixel_to_world method? That is, how do we know the axes order of the spectrum object if we have no prior knowledge of how it was loaded? A lot of code in Cubeviz assumes spectral axis got moved to last.

[keflavich]

It should be by explicitly identifying the axis type, e.g., spectrum.wcs.spectral.world_to_pixel(value).

[pllim]

Re: #999 (comment) -- How would the spatial component work in terms of pixel_to_world for spectral cube?

[keflavich]

re: #999 (comment), should be cube.wcs.celestial.pixel_to_world(...)

[pllim]

.celestial isn't really part of APE 14 and I think it only works for FITS WCS.

[keflavich]

Ah, right. Then I'm not sure. That's the API I'd like to see & use, and it's what I've had in mind during our conversations on the topic, but maybe it's a little more complicated because we need to use the more general APE14 version that can handle spatial axes unassociated with celestial coordinates.

[pllim]

@rosteen , any input on #999 (comment) ?

[rosteen]

If we do this, how do we know how to access the wcs.pixel_to_world method? That is, how do we know the axes order of the spectrum object if we have no prior knowledge of how it was loaded? A lot of code in Cubeviz assumes spectral axis got moved to last.

@pllim I've done some investigating on this to refresh my memory, since I still don't mess with GWCS that often. For the spectral axis at least, with this PR you can use spec.spectral_axis_index to find out where the spectral axis is. You can also always use spec.wcs.world_axis_physical_types and look for the em.[X] axis (for a cube the output will be something like ('pos.eq.ra', 'pos.eq.dec', 'em.wl')). But in general, you can always just provide all three axis locations to pixel_to_world and then take either the SkyCoord or SpectralCoord from the output as desired - I imagine you're considered about mouseover in Cubeviz, in which case we would just grab the SkyCoord from (for example) [<SkyCoord (ICRS): (ra, dec) in deg (339.01582773, 33.97526219)>, <SpectralCoord 4.9124001 um>]. It will definitely need some updates to the Jdaviz code, but I think it won't be too hard to work around.

Is that a sufficient answer to your question? (and did I actually answer the question you were asking? 😆 )

[rosteen]

FYI I edited the initial description to include the context of why I think we should do this - I'll also add some more information to the narrative docs about the change.

[pllim]

did I actually answer the question you were asking?

Not sure. Yes, it is about the mouseover. I just fixed a subtle bug over there (spacetelescope/jdaviz#2009). It is very subtle because even if you mess up the axes order in your pixel_to_world, it still give you something that looks valid but is actually wrong. So my question is, how do I map event (x, y) from the frontend to the actual pixel_to_world call if the order can be arbitrary?

Indeed, spec.spectral_axis_index might solve the problem of grabbing the spectral axis, but I am still not sure about spatial. Right now, we assume it is always (x, y, z) (since specutils forces z to be last). So, if spec.spectral_axis_index comes back to be 1 for a data cube, what am I supposed to think? Is that possible? What if it is 0, is it always (z, y, x) for that case, or can it also be (z, x, y)?

[rosteen]

Closing this in favor of #1033 (merging into a new specutils 2.0 development branch for now rather than main).