diff --git a/gwcs/api.py b/gwcs/api.py
index b698f40d..b23ee868 100644
--- a/gwcs/api.py
+++ b/gwcs/api.py
@@ -51,14 +51,7 @@ def world_axis_physical_types(self):
         arbitrary string.  Alternatively, if the physical type is
         unknown/undefined, an element can be `None`.
         """
-        # A CompositeFrame orders the output correctly based on axes_order.
-        if isinstance(self.output_frame, cf.CompositeFrame):
-            return self.output_frame.axis_physical_types
-
-        # If we don't have a CompositeFrame, where this is taken care of for us,
-        # we need to make sure we re-order the output to match the transform.
-        # The underlying frames don't reorder themselves because axes_order is global.
-        return tuple(self.output_frame.axis_physical_types[i] for i in self.output_frame.axes_order)
+        return self.output_frame.axis_physical_types
 
     @property
     def world_axis_units(self):
diff --git a/gwcs/coordinate_frames.py b/gwcs/coordinate_frames.py
index ab6683fe..f63eb8b1 100644
--- a/gwcs/coordinate_frames.py
+++ b/gwcs/coordinate_frames.py
@@ -481,7 +481,8 @@ def axis_physical_types(self):
 
         These physical types are the types in frame order, not transform order.
         """
-        return self._prop.axis_physical_types or self._default_axis_physical_types
+        apt = self._prop.axis_physical_types or self._default_axis_physical_types
+        return self._sort_property(apt)
 
     @property
     def world_axis_object_classes(self):
@@ -500,7 +501,7 @@ class CelestialFrame(CoordinateFrame):
     Representation of a Celesital coordinate system.
 
     This class has a native order of longitude then latitude, meaning
-    ``axes_names``, ``unit`` should be lon, lat ordered.  If your transform is
+    ``axes_names``, ``unit`` and ``axis_physical_types`` should be lon, lat ordered.  If your transform is
     in a different order this should be specified with ``axes_order``.
 
     Parameters
@@ -515,6 +516,8 @@ class CelestialFrame(CoordinateFrame):
         Names of the axes in this frame.
     name : str
         Name of this frame.
+    axis_physical_types : list
+        The UCD 1+ physical types for the axes, in frame order (lon, lat).
     """
 
     def __init__(self, axes_order=None, reference_frame=None,
diff --git a/gwcs/tests/conftest.py b/gwcs/tests/conftest.py
index 284af501..a7f19182 100644
--- a/gwcs/tests/conftest.py
+++ b/gwcs/tests/conftest.py
@@ -54,7 +54,7 @@ def gwcs_2d_spatial_reordered():
     A simple one step spatial WCS, in ICRS with a 1 and 2 px shift.
     """
     out_frame = cf.CelestialFrame(reference_frame=coord.ICRS(),
-                                   axes_order=(1, 0))
+                                  axes_order=(1, 0))
     return wcs.WCS(model_2d_shift | models.Mapping((1, 0)), input_frame=detector_2d, output_frame=out_frame)
 
 
diff --git a/gwcs/tests/test_coordinate_systems.py b/gwcs/tests/test_coordinate_systems.py
index 2d31603b..a54091e7 100644
--- a/gwcs/tests/test_coordinate_systems.py
+++ b/gwcs/tests/test_coordinate_systems.py
@@ -475,3 +475,29 @@ def test_ucd1_to_ctype(caplog):
         assert ctype_to_ucd[v] == k
 
     assert inv_map['new.repeated.type'] in new_ctype_to_ucd
+
+
+def test_celestial_ordering():
+   c1 = cf.CelestialFrame(
+       reference_frame=coord.ICRS(),
+       axes_order=(0, 1),
+       axes_names=("lon", "lat"),
+       unit=(u.deg, u.arcsec),
+       axis_physical_types=("custom:lon", "custom:lat"),
+   )
+   c2 = cf.CelestialFrame(
+       reference_frame=coord.ICRS(),
+       axes_order=(1, 0),
+       axes_names=("lon", "lat"),
+       unit=(u.deg, u.arcsec),
+       axis_physical_types=("custom:lon", "custom:lat"),
+   )
+
+   assert c1.axes_names == ("lon", "lat")
+   assert c2.axes_names == ("lat", "lon")
+
+   assert c1.unit == (u.deg, u.arcsec)
+   assert c2.unit == (u.arcsec, u.deg)
+
+   assert c1.axis_physical_types == ("custom:lon", "custom:lat")
+   assert c2.axis_physical_types == ("custom:lat", "custom:lon")