Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix bug in translating 3D Spectrum1D with only spectral axis defined #61

Closed
wants to merge 5 commits into from

Conversation

rosteen
Copy link
Contributor

@rosteen rosteen commented Jan 25, 2022

Description

If a 3D flux is input to a Spectrum1D along with a spectral_axis (rather than a full 3D WCS), the resulting object will have a spectral GWCS of only one dimension. Personally I think this behavior is not ideal and might need rethinking at some point, but in either case this fix will stop the translator from trying to swap axes on the WCS in this case, thus fixing (I hope) the bug @pllim ran into in spacetelescope/jdaviz#1040 (comment).

Fixes #60

@codecov
Copy link

codecov bot commented Jan 25, 2022

Codecov Report

Merging #61 (312a840) into main (7f3229d) will decrease coverage by 1.23%.
The diff coverage is 84.09%.

Impacted file tree graph

@@            Coverage Diff             @@
##             main      #61      +/-   ##
==========================================
- Coverage   95.80%   94.56%   -1.24%     
==========================================
  Files          15       15              
  Lines        1145     1123      -22     
==========================================
- Hits         1097     1062      -35     
- Misses         48       61      +13     
Impacted Files Coverage Δ
glue_astronomy/translators/spectrum1d.py 89.91% <80.00%> (-2.66%) ⬇️
...lue_astronomy/translators/tests/test_spectrum1d.py 92.85% <92.85%> (-7.15%) ⬇️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 7f3229d...312a840. Read the comment docs.

@rosteen rosteen marked this pull request as draft January 25, 2022 16:30
@pllim

This comment was marked as resolved.

@pllim

This comment was marked as resolved.

@rosteen
Copy link
Contributor Author

rosteen commented Jan 27, 2022

@astrofrog I was working on the route you suggested of using NDCube's CompoundLowLevelWCS for a general solution of padding the WCS needed for additional axes, but I'm unsure about how many of the properties you had coded into your 2D solution need to be kept in my version. Seems like we do need a high level WCS object here, but just using HighLevelWCSWrapper around the CompoundLowLevelWCS isn't sufficient. If you have any advice/comments based on where I'm at so far I'd appreciate it!

@pllim

This comment was marked as resolved.

@pllim
Copy link
Contributor

pllim commented Feb 14, 2022

Okay. I got it. With this PR, it works if I pass in only the flux like this Spectrum1D(flux=np.random.random((10, 10, 5))), where 5 is the length of the spectral axis. If I pass in spectral_axis explicitly, it crashes. Does the shape makes sense to you? Is spectral axis supposed to be in arr.shape[2]? I still get confused by all the swapping that happens (astropy/specutils#923).

Screenshot 2022-02-14 164836

https://github.com/astropy/specutils/blob/7f7b55272b7e36d827c130bb2bfa0ec6078f41a3/specutils/spectra/spectrum1d.py#L265-L269

pllim added a commit to pllim/jdaviz that referenced this pull request Feb 14, 2022
@dhomeier
Copy link
Contributor

https://github.com/glue-viz/glue-astronomy/blob/7f3229d01bdc2bff64fdb1a83738447cf696c59d/glue_astronomy/translators/tests/test_spectrum1d.py#L204-#234
will need to be changed to flux.T/spec.flux.T if that is the intended behaviour now. Still feels a bit weird...

pllim added a commit to pllim/jdaviz that referenced this pull request Feb 16, 2022
pllim added a commit to pllim/jdaviz that referenced this pull request Feb 16, 2022
pllim added a commit to pllim/jdaviz that referenced this pull request Feb 21, 2022
pllim added a commit to pllim/jdaviz that referenced this pull request Feb 21, 2022
@pllim
Copy link
Contributor

pllim commented Feb 21, 2022

Is there a way to push this forward? I need this for spacetelescope/jdaviz#1040 .

@dhomeier
Copy link
Contributor

I could go ahead with adapting the tests, but @astrofrog and I are quite uncertain at this point what the logic of reordering axes here again is. specutils as far as I see it is consistent in using the last dimension as spectral axis, so it's rather confusing to have the glue loader swap this again, and (only?) for [G]WCS spectra.

@pllim
Copy link
Contributor

pllim commented Feb 21, 2022

In our FITS data, the spectral axis seems to be in NAXIS3, which in Python-speak, is axis=0, so maybe it is specutils who is putting it in the wrong end. Though perhaps that is necessary for performance reason? 🤷

@dhomeier
Copy link
Contributor

dhomeier commented Feb 21, 2022

NAXIS3 -> axis=0 is just the conversion from Fortran- to C-ordering, right?
https://github.com/astropy/specutils/blob/af729120269e776e40d8a5e6c2e270617d7a8d6e/specutils/spectra/spectrum1d.py#L199-L205
notes that for the WCS the order is turned around again, but I think that is simply because by convention the WCS is kept in Fortran order. While the two have to match, I so not see a technical reason why the spectral axis has to be last, and performance-wise it would certainly be better to avoid those np.swapaxes.
I suppose the spectral axis was simply found more frequently in NAXIS1, and to avoid even more confusion that was set as a rule. According to the comment in the data_translator the convention for glue is exactly the opposite, so I guess swapping the axes again (and fixing the tests accordingly) is unavoidable at this point.
If the "glue-order" is common for large data cube, I think it is worth considering to allow that in specutils as well (e.g. MaNGA also uses NAXIS3, and the axes swapping seems to introduce a quite noticeable delay).

@pllim
Copy link
Contributor

pllim commented Feb 22, 2022

Yes, I think @rosteen mentioned that it would be nice if specutils is flexible enough to allow specutils in whatever axis but that's unlikely to happen any time soon. :(

Comment on lines -108 to -109
def world_axis_names(self):
return (UCD_TO_SPECTRAL_NAME[self.spectral_wcs.world_axis_physical_types[0]], 'Offset')
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is changing the label of the returned components in the 2d example to 'World 0','World 1', breaking test_spectrum1d_2d_data further down. Maybe the new naming is more logical, but that should be documented somewhere (and at least World 1 does not sound like an improvement over Frequency).

@dhomeier
Copy link
Contributor

@rosteen I have changed the expected labels in the test for now, so it can progress further, but does not necessarily mean they should remain like this.
The current failure however raises already when just dereferencing data['World 1'] as an array-like object ( in the call to _return_single_array), which seems to indicate something is still more fundamentally broken in the implementation (or is exposing a new bug in astropy.wcs?).

@pllim
Copy link
Contributor

pllim commented Feb 22, 2022

What does this error mean?

        assert_equal(data['World 0'], [[10, 10, 10], [10, 10, 10]])
>       assert data['World 1'].shape == (2, 3)

glue_astronomy/translators/tests/test_spectrum1d.py:255: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
.tox/py38-test/lib/python3.8/site-packages/glue/core/data.py:571: in __getitem__
    return self.get_data(key, view=view)
.tox/py38-test/lib/python3.8/site-packages/glue/core/data.py:1366: in get_data
    result = comp.data
.tox/py38-test/lib/python3.8/site-packages/glue/core/component.py:215: in data
    return self._calculate()
.tox/py38-test/lib/python3.8/site-packages/glue/core/component.py:313: in _calculate
    world_coords = pixel2world_single_axis(self._data.coords,
.tox/py38-test/lib/python3.8/site-packages/glue/core/coordinate_helpers.py:56: in pixel2world_single_axis
    result = wcs.pixel_to_world_values(*pixel)
.tox/py38-test/lib/python3.8/site-packages/ndcube/wcs/wrappers/compound_wcs.py:112: in pixel_to_world_values
    world_arrays_sub = w.pixel_to_world_values(*pixel_arrays_sub)
.tox/py38-test/lib/python3.8/site-packages/astropy/wcs/wcsapi/fitswcs.py:322: in pixel_to_world_values
    world = self.all_pix2world(*pixel_arrays, 0)

@dhomeier
Copy link
Contributor

Per my speculation above, it's two levels further down ;-)

../../opt/lib/python3.9/site-packages/astropy/wcs/wcs.py:1351: in all_pix2world
    return self._array_converter(
../../opt/lib/python3.9/site-packages/astropy/wcs/wcs.py:1328: in _array_converter
    return _return_single_array(xy, origin)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

xy = array([[0, 0, 0]]), origin = 0

    def _return_single_array(xy, origin):
        if xy.shape[-1] != self.naxis:
>           raise ValueError(
                "When providing two arguments, the array must be "
                "of shape (N, {})".format(self.naxis))
E           ValueError: When providing two arguments, the array must be of shape (N, 1)

This apparently already when trying to access the shape attribute. Looks wicked. Or is it just that the shape is inconsistent with the no. of dimensions somewhere else?

pllim added a commit to pllim/jdaviz that referenced this pull request Feb 23, 2022
pllim added a commit to pllim/jdaviz that referenced this pull request Feb 23, 2022
@pllim
Copy link
Contributor

pllim commented Feb 23, 2022

Definitely something is still wonky here, unfortunately. Even if you ignore the data['World 1'] parts of the test (it seems to go into ndcube at some point and I got lost there), it crashes later at:

>>> s, o = data.coords.pixel_to_world(1, 2)
AttributeError: 'PaddedSpectrumWCS' object has no attribute 'pixel_to_world'

and

>>> px, py = data.coords.world_to_pixel(s, o)
AttributeError: 'PaddedSpectrumWCS' object has no attribute 'world_to_pixel'

The spec_new stuff works.

@pllim
Copy link
Contributor

pllim commented Feb 23, 2022

Maybe subclassing from CompoundLowLevelWCS is too low of a level?

@dhomeier
Copy link
Contributor

Maybe subclassing from CompoundLowLevelWCS is too low of a level?

Yea, the pixel_to_world etc. are only provided by the HighLevel version, but that in turn is not supported by gwcs afair.
Dunno if it can all be moved to the pixel_to_world_values level...

pllim added a commit to pllim/jdaviz that referenced this pull request Feb 23, 2022
pllim added a commit to pllim/jdaviz that referenced this pull request Feb 23, 2022
Copy link
Contributor

@pllim pllim left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As per PO feedback at sprint review today, they want the spectral axis to be dimensionless rather than pixels. Is it possible to do that upstream here? If not, I'll have to hack something up downstream in Jdaviz.

def __init__(self, spectral_wcs, n_extra_axes):
self.spectral_wcs = spectral_wcs
frame1 = CoordinateFrame(n_extra_axes, ['SPATIAL']*n_extra_axes,
np.arange(n_extra_axes), unit=[u.pix]*n_extra_axes,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
np.arange(n_extra_axes), unit=[u.pix]*n_extra_axes,
np.arange(n_extra_axes), unit=[u.dimensionless_unscaled]*n_extra_axes,

np.arange(n_extra_axes), unit=[u.pix]*n_extra_axes,
name="Dummy1")
frame2 = CoordinateFrame(n_extra_axes, ['SPATIAL']*n_extra_axes,
np.arange(n_extra_axes), unit=[u.pix]*n_extra_axes,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
np.arange(n_extra_axes), unit=[u.pix]*n_extra_axes,
np.arange(n_extra_axes), unit=[u.dimensionless_unscaled]*n_extra_axes,

pllim added a commit to pllim/jdaviz that referenced this pull request Feb 28, 2022
pllim added a commit to pllim/jdaviz that referenced this pull request Feb 28, 2022
@rosteen
Copy link
Contributor Author

rosteen commented Mar 1, 2022

As per PO feedback at sprint review today, they want the spectral axis to be dimensionless rather than pixels. Is it possible to do that upstream here? If not, I'll have to hack something up downstream in Jdaviz.

I seem to remember a decision specifically not to handle dimensionless "spectral" axes in specutils, with the compromise that we'd handle pixel units. I could be remembering wrong though, I'll look back through some notes/git discussion.

Btw, thanks @dhomeier and @pllim for looking at this, I left it in a state where I'd run into the limit of my WCS knowledge. I'm trying to find time this sprint to get back to it.

@pllim
Copy link
Contributor

pllim commented Mar 1, 2022

dimensionless "spectral" axes

Okay, just keep it as pixel here. I dealt with that part downstream. Thanks!

pllim added a commit to pllim/jdaviz that referenced this pull request Mar 1, 2022
pllim added a commit to pllim/jdaviz that referenced this pull request Mar 1, 2022
pllim added a commit to pllim/jdaviz that referenced this pull request Mar 2, 2022
pllim added a commit to pllim/jdaviz that referenced this pull request Mar 2, 2022
pllim added a commit to pllim/jdaviz that referenced this pull request Mar 9, 2022
pllim added a commit to pllim/jdaviz that referenced this pull request Mar 9, 2022
pllim added a commit to pllim/jdaviz that referenced this pull request Mar 10, 2022
pllim added a commit to pllim/jdaviz that referenced this pull request Mar 10, 2022
pllim added a commit to pllim/jdaviz that referenced this pull request Mar 11, 2022
pllim added a commit to pllim/jdaviz that referenced this pull request Mar 11, 2022
pllim added a commit to pllim/jdaviz that referenced this pull request Mar 15, 2022
pllim added a commit to pllim/jdaviz that referenced this pull request Mar 15, 2022
@pllim
Copy link
Contributor

pllim commented Mar 22, 2022

This can be closed without merged. Superseded by #68

@astrofrog astrofrog closed this Mar 22, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Spectrum1D translator cannot handle spectral_axis in 3D cube
4 participants