Add a test for axes ordering with CelestialFrame

fixes spacetelescope#269

gwcs/coordinate_frames.py

@@ -436,6 +436,11 @@ def world_axis_object_classes(self):
     def world_axis_object_components(self):
         return [(f"{at}{i}" if i != 0 else at, 0, 'value') for i, at in enumerate(self._axes_type)]
 
+    @property
+    def _native_world_axis_object_components(self):
+        """Defines the target component ordering (i.e. not taking into account axes_order)"""
+        return self.world_axis_object_components
+
 
 class CelestialFrame(CoordinateFrame):
     """
@@ -509,10 +514,18 @@ def world_axis_object_classes(self):
              'unit': self.unit})}
 
     @property
-    def world_axis_object_components(self):
+    def _native_world_axis_object_components(self):
         return [('celestial', 0, 'spherical.lon'),
                 ('celestial', 1, 'spherical.lat')]
 
+    @property
+    def world_axis_object_components(self):
+        # Sort the native waoc by the axes order. The axes order may have jumps
+        # in it if there are other frames in between the components.
+        ordered = np.array(self._native_world_axis_object_components,
+                           dtype=object)[np.argsort(self.axes_order)]
+        return list(map(tuple, ordered))
+
 
 class SpectralFrame(CoordinateFrame):
     """
@@ -709,7 +722,7 @@ def _wao_renamed_components_iter(self):
         mapper = self._wao_classes_rename_map
         for frame in self.frames:
             renamed_components = []
-            for comp in frame.world_axis_object_components:
+            for comp in frame._native_world_axis_object_components:
                 comp = list(comp)
                 rename = mapper[frame].get(comp[0])
                 if rename:

gwcs/tests/test_wcs.py

@@ -1272,9 +1272,14 @@ def test_spatial_spectral_stokes():
 
 
 def test_split_frame_wcs():
-    # Setup a model where the pixel & world axes are (lat, wave, lon)
-    spatial = models.Multiply(10*u.arcsec/u.pix) & models.Multiply(15*u.arcsec/u.pix)  # pretend this is a spatial model
+    # Setup a WCS where the pixel & world axes are (lat, wave, lon)
+
+    # We setup a model which is pretending to be a celestial transform. Note
+    # that we are pretending that this model is ordered lon, lat because that's
+    # what the projections require in astropy.
+    spatial = models.Multiply(20*u.arcsec/u.pix) & models.Multiply(15*u.arcsec/u.pix)
     compound = models.Linear1D(intercept=0*u.nm, slope=10*u.nm/u.pix) & spatial
+    # This forward transforms uses mappings to be (lat, wave, lon)
     forward = models.Mapping((1, 2, 0)) | compound | models.Mapping((2, 0, 1))
 
     # Setup the output frame
@@ -1287,14 +1292,52 @@ def test_split_frame_wcs():
                                      axes_order=list(range(3)), unit=[u.pix]*3)
 
     iwcs = wcs.WCS(forward, input_frame, output_frame)
-    input_pixel = [0*u.pix, 1*u.pix, 2*u.pix]
+    input_pixel = [1*u.pix, 2*u.pix, 3*u.pix]
     output_world = iwcs.pixel_to_world_values(*input_pixel)
     output_pixel = iwcs.world_to_pixel_values(*output_world)
     assert_allclose(output_pixel, u.Quantity(input_pixel).to_value(u.pix))
 
+    expected_world = [15*u.arcsec, 20*u.nm, 60*u.arcsec]
+    for expected, output in zip(expected_world, output_world):
+        assert_allclose(output, expected.value)
+
     world_obj = iwcs.pixel_to_world(*input_pixel)
     assert isinstance(world_obj[0], coord.SkyCoord)
     assert isinstance(world_obj[1], coord.SpectralCoord)
 
+    assert u.allclose(world_obj[0].spherical.lat, expected_world[0])
+    assert u.allclose(world_obj[0].spherical.lon, expected_world[2])
+    assert u.allclose(world_obj[1], expected_world[1])
+
     obj_pixel = iwcs.world_to_pixel(*world_obj)
     assert_allclose(obj_pixel, u.Quantity(input_pixel).to_value(u.pix))
+
+
+def test_reordered_celestial():
+    # This is a spatial model which is ordered lat, lon for the purposes of this test.
+    spatial = models.Multiply(20*u.deg/u.pix) & models.Multiply(15*u.deg/u.pix)
+
+    celestial_frame = cf.CelestialFrame(axes_order=(1, 0), unit=(u.deg, u.deg),
+                                        reference_frame=coord.ICRS())
+
+    input_frame = cf.CoordinateFrame(2, ["PIXEL"]*2,
+                                     axes_order=list(range(2)), unit=[u.pix]*2)
+
+    iwcs = wcs.WCS(spatial, input_frame, celestial_frame)
+
+    input_pixel = [1*u.pix, 3*u.pix]
+    output_world = iwcs.pixel_to_world_values(*input_pixel)
+    output_pixel = iwcs.world_to_pixel_values(*output_world)
+    assert_allclose(output_pixel, u.Quantity(input_pixel).to_value(u.pix))
+
+    expected_world = [20*u.deg, 45*u.deg]
+    assert_allclose(output_world, [e.value for e in expected_world])
+
+    world_obj = iwcs.pixel_to_world(*input_pixel)
+    assert isinstance(world_obj, coord.SkyCoord)
+
+    assert u.allclose(world_obj.spherical.lat, expected_world[0])
+    assert u.allclose(world_obj.spherical.lon, expected_world[1])
+
+    obj_pixel = iwcs.world_to_pixel(world_obj)
+    assert_allclose(obj_pixel, u.Quantity(input_pixel).to_value(u.pix))