Is axes_order supposed to allow reordering inside a frame or a Composite frame

[Cadair]

I understand it can reorder things in CompositeFrame i.e. to have an array which is (lat, spectral, lon) or similar, but what about in a WCS with just a CelestialFrame which is ordered lat, lon. Currently this dosen't seem to work:

In [14]: import astropy.coordinates as coord                                                                                                                                                                                                                   

In [15]: import gwcs                                                                                                                                                                                                                                           

In [16]: from astropy.modeling import models                                                                                                                                                                                                                   

In [17]: model_2d_shift = models.Shift(1) & models.Shift(2)                                                                                                                                                                                                    

In [18]: out_frame = gwcs.coordinate_frames.CelestialFrame(reference_frame=coord.ICRS(), axes_order=(1, 0))                                                                                                                                                    

In [19]: w = gwcs.WCS(model_2d_shift, output_frame=out_frame)                                                                                                                                                                                                  

In [20]: model_2d_shift(3, 3)                                                                                                                                                                                                                                  
Out[20]: (4.0, 5.0)

In [21]: w(3, 3, with_units=True)                                                                                                                                                                                                                              
Out[21]: 
<SkyCoord (ICRS): (ra, dec) in deg
    (4., 5.)>

Given the flipped output order I would expect this to return

<SkyCoord (ICRS): (ra, dec) in deg
    (5., 4.)>

is this a bug or expected behaviour?

[nden]

The CoordinateFrame is only informational. You can control the ordering through the transforms

model_2d_shift = models.Shift(1) & models.Shift(2)  | models.Mapping((1, 0))

It may become messy to allow this changes on both objects because then we have to think which one takes precedence. I think it's cleaner to have the CoordinateFrame be only informational container and make sure the transform gets the correct result and ordering. Does this sound OK?

[Cadair]

I see the point, but then axes_order basically shouldn't exist? You should just have frame_order on composite frame? (which could just be the ordering of the list?)

this fixture in #217 sets the axes ordering over two frames, where it puts the spectral frame between the two spatial axes (something which is at least possible in real data), and it seems to work fine (GWCS.__call__ returns them in the axes_order order as I can tell). So at least this is inconsistent as it currently stands?

[nden]

The fixture looks fine to me. I think that's what I would expect. The ordering is "made to work" in the transform. output_frame.axes_order is defined to match the transform. And the purpose is to give a user a way to tell which is which. I guess I don't understand what is inconsistent.

 w.output_frame.axes_order                                                                   
Out[18]: (2, 0, 1)

In [19]: w.output_frame.axes_names                                                                   
Out[19]: ('lat', '', 'lon')

[Cadair]

Ah! Ok, I see, so the idea is for the axes_order to match the transform, that makes perfect sense.

[chris-simpson]

I don't think this is resolved. axes_order has some effect in certain circumstances but not in others, and I think that's unwanted behaviour.

Below looks like an example similar to @nden's post above but note what happens when I call the gwcs object with with_units=True.

>>> icrs = cf.CelestialFrame(reference_frame=coord.ICRS(), axes_order=(2, 0))
>>> spec = cf.SpectralFrame(name='wave', unit=[u.m,], axes_order=(1,),
...                         axes_names=('lambda',))
>>> output_frame = cf.CompositeFrame(frames=[icrs, spec])
>>> output_frame.axes_order
(2, 0, 1)
>>> output_frame.axes_names
('lat', 'lambda', 'lon')
>>> w = wcs.WCS(forward_transform=astmodels.Identity(3), output_frame=output_frame)
>>> w(1,2,3)
(1.0, 2.0, 3.0)
>>> w(1,2,3, with_units=True)
[<SkyCoord (ICRS): (ra, dec) in deg
    (3., 1.)>, <Quantity 2. m>]

The CompositeFrame.coordinates() method reorders its values according to the order of its component frames. It's obviously difficult to place a single SkyCoord() object in non-adjacent locations of the output, but the fact that with_units gives the output in a different order seems like too much change. There are several related issues to this. Here I'll create a simple composite frame with a spectral axis and a generic spatial axis (e.g., a spectroscopic longslit).

>>> slit = cf.CoordinateFrame(naxes=1, axes_type='SPATIAL', axes_order=(0,),
...                           axes_names=('slit',), unit=u.pix)
>>> spec = cf.SpectralFrame(name='wave', unit=[u.m,], axes_order=(1,),
...                         axes_names=('lambda',))
>>> output_frame = cf.CompositeFrame(frames=[slit, spec])
>>> w = wcs.WCS(forward_transform=astmodels.Identity(2), output_frame=output_frame)
>>> w(1, 2, with_units=True)
[<Quantity 1. m>, (<Quantity 2. pix>,)]

which is as it should be. But what if I define my CompositeFrame in the other order? After all, why should it matter?

>>> spec = cf.SpectralFrame(name='wave', unit=[u.m,], axes_order=(0,),
...                         axes_names=('lambda',))
>>> slit = cf.CoordinateFrame(naxes=1, axes_type='SPATIAL', axes_order=(1,),
...                           axes_names=('slit',), unit=u.pix)
>>> output_frame = cf.CompositeFrame(frames=[spec, slit])
>>> w = wcs.WCS(forward_transform=astmodels.Identity(2), output_frame=output_frame)
>>> w(2, 1)
(2.0, 1.0)
>>> w(2, 1, with_units=True)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/opt/anaconda3/envs/test/lib/python3.7/site-packages/gwcs/wcs.py", line 276, in __call__
    result = self.output_frame.coordinates(*result)
  File "/opt/anaconda3/envs/test/lib/python3.7/site-packages/gwcs/coordinate_frames.py", line 591, in coordinates
    coo.append(frame.coordinates(arg))
  File "/opt/anaconda3/envs/test/lib/python3.7/site-packages/gwcs/coordinate_frames.py", line 273, in coordinates
    args = [args[i] for i in self.axes_order]
  File "/opt/anaconda3/envs/test/lib/python3.7/site-packages/gwcs/coordinate_frames.py", line 273, in <listcomp>
    args = [args[i] for i in self.axes_order]
IndexError: tuple index out of range

CompositeFrame.coordinates() only sends the value 2 to slit.coordinates() but slit.coordinates() needs to have been sent at least 2 values because its own axes_order=(1,). The original post notes that CelestialFrame ignores its own axes_order so why do CompositeFrame and CoordinateFrame not also ignore it?

I'm not sure what the desired behaviour here is. Personally, I'd like to define both frames with axes_order=(0,) (so each is valid on its own) and have CompositeFrame assign axes numbers according to the order in which the frames are located in its input list. This would prevent a multi-dimensional frame (e.g., CelestialFrame, Frame2D) from having non-adjacent axes but I struggle to see a use case for that.

And, going back to the original post,

>>> icrs = cf.CelestialFrame(reference_frame=coord.ICRS(), axes_order=(1, 0))
>>> icrs.coordinates(90, 20)
<SkyCoord (ICRS): (ra, dec) in deg
    (90., 20.)>

is counterintuitive at best. But if I do this

>>> icrs = cf.CelestialFrame(reference_frame=coord.ICRS(), axes_order=(1, 0))
>>> spec = cf.SpectralFrame(name='wave', unit=[u.m,], axes_order=(2,),
...                         axes_names=('lambda',))
>>> output_frame = cf.CompositeFrame(frames=[icrs, spec])
>>> output_frame.coordinates(90, 20, 1)
[<SkyCoord (ICRS): (ra, dec) in deg
    (20., 90.)>, <Quantity 1. m>]

now it pays attention to axes_order because it's in a CompositeFrame. This behaviour looks inconsistent to me.

[jehturner]

This does seem rather unexpected at best. How is the inverse used in that case with_units? I'm having trouble seeing how to feed those outputs back through the other way.

Personally, I'd like to define both frames with axes_order=(0,) (so each is valid on its own) and have CompositeFrame assign axes numbers according to the order in which the frames are located in its input list. This would prevent a multi-dimensional frame (e.g., CelestialFrame, Frame2D) from having non-adjacent axes but I struggle to see a use case for that.

While agreeing with your sentiment, @chris-simpson, if I understand you correctly there are certainly real use cases for non-adjacent spatial axes, notably visualizing slices through an IFU data cube (I don't think it's the end of the World if that case doesn't work for now, but it should not be precluded by design).

[chris-simpson]

As a pretty fundamental axiom, I don't think the order of the outputs should be affected by the value of with_units. Since with_units will return a SkyCoord instance where appropriate, this requires that the celestial frame axes be adjacent in the world system. It's irrelevant how they're arranged in the data (because you can use Mapping) and the arrangement in the world system appears to be left to the user. So we just constrain that a little in order to be able to predict the order of the outputs, which seems like a requirement. That seems like a reasonable trade-off to me. But I may be missing something, which is why I'm not submitting a PR to "fix" this.

[jehturner]

Yes, constraining the order in the outputs doesn't seem unreasonable (given that you can't enforce a 1-to-1 mapping of pixel to world co-ordinates in general anyway).

[nden]

I looked at this issue a bit in an attempt to define the behavior and tie loose ends. First how we got here.

Gwcs started as a numerical only library but there was always a desire to represent the results as astropy objects and to use units with transforms. This resulted in the additional parameter with_units which calls the corresponding CoordinateFrame.coordinates method to construct the corresponding astropy object.
Later the WCS common API was developed with two methods to return either numerical values or astropy objects.
gwcs.pixel_to_world uses with_units=True to return astropy objects.

In the case of Celestial frame a SkyCoord object is always ordered as (lon, lat) which is why the ordering is different than the numerical result. I think for those who work with astropy objects this is natural and if accessing the results as SkyCoord.lon, SkyCoord.lat they should be able to get the correct results which means they need to be reordered on individual and Composite frames, even if they are not adjacent.

There's another twist to this, as I discovered while looking at this issue. The methods _world_axis_object_components and _world_axis_object_classes are used by the HighLevelWCSWrapper to essentially do the same things as WCS.pixel_to_world and WCS.world_to_pixel taking a different path in the code. IMO it'd be best for the HighLevelWCSWrapper to check if pixel_to_world and world_to_pixel are available and call them instead of using alternate methods. If this is done, I think the rest comes into place easily and works correctly. BTW the SlicedLowLevelWCS tests pass with this change.

I collected use cases from issues and tests in the notebook below and the results show how I think it should work.
I have locally made the change to HighLevelWCSWrapper and some other changes to gwcs.
I would appreciate if we agree on the behavior before I push changes here and to astropy. Also it won't hurt if someone else looks at the results and ordering (there was too much axes swapping in my brain while I did this).
If we agree on this quickly enough, I will attempt to put it in astropy 4.2.

https://github.com/nden/documentation/blob/master/Reorder_axes.ipynb

@Cadair @chris-simpson @jehturner @mcara