From 846b0cff6493c1c586bd186a984482be94828bfb Mon Sep 17 00:00:00 2001
From: Stuart Mumford <stuart@cadair.com>
Date: Tue, 20 Jun 2023 15:24:42 +0100
Subject: [PATCH] Rewrite coordinate systems tests for APE 14

This highlights some API changes here.
---
 gwcs/api.py                           |   4 +-
 gwcs/coordinate_frames.py             |  64 +++++-----
 gwcs/tests/test_api.py                |   4 +-
 gwcs/tests/test_coordinate_systems.py | 167 +++++++++++++++-----------
 gwcs/tests/test_wcs.py                |  29 +++++
 5 files changed, 160 insertions(+), 108 deletions(-)

diff --git a/gwcs/api.py b/gwcs/api.py
index ca606f4d..ff633143 100644
--- a/gwcs/api.py
+++ b/gwcs/api.py
@@ -256,11 +256,11 @@ def serialized_classes(self):
 
     @property
     def world_axis_object_classes(self):
-        return self.output_frame._world_axis_object_classes
+        return self.output_frame.world_axis_object_classes
 
     @property
     def world_axis_object_components(self):
-        return self.output_frame._world_axis_object_components
+        return self.output_frame.world_axis_object_components
 
     @property
     def pixel_axis_names(self):
diff --git a/gwcs/coordinate_frames.py b/gwcs/coordinate_frames.py
index 9d2073b5..6654e1ff 100644
--- a/gwcs/coordinate_frames.py
+++ b/gwcs/coordinate_frames.py
@@ -249,7 +249,7 @@ def axis_physical_types(self):
 
     @property
     @abc.abstractmethod
-    def _world_axis_object_classes(self):
+    def world_axis_object_classes(self):
         """
         The APE 14 object classes for this frame.
 
@@ -260,7 +260,7 @@ def _world_axis_object_classes(self):
 
     @property
     @abc.abstractmethod
-    def _world_axis_object_components(self):
+    def world_axis_object_components(self):
         """
         The APE 14 object components for this frame.
 
@@ -444,14 +444,14 @@ def axis_physical_types(self):
         return self._axis_physical_types or self._default_axis_physical_types
 
     @property
-    def _world_axis_object_classes(self):
+    def world_axis_object_classes(self):
         return {f"{at}{i}" if i != 0 else at: (u.Quantity,
                      (),
                      {'unit': unit})
                 for i, (at, unit) in enumerate(zip(self._axes_type, self.unit))}
 
     @property
-    def _world_axis_object_components(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)]
 
 
@@ -519,7 +519,7 @@ def _default_axis_physical_types(self):
             return tuple("custom:{}".format(t) for t in self.axes_names)
 
     @property
-    def _world_axis_object_classes(self):
+    def world_axis_object_classes(self):
         return {'celestial': (
             coord.SkyCoord,
             (),
@@ -527,7 +527,7 @@ def _world_axis_object_classes(self):
              'unit': self.unit})}
 
     @property
-    def _world_axis_object_components(self):
+    def world_axis_object_components(self):
         return [('celestial', 0, 'spherical.lon'),
                 ('celestial', 1, 'spherical.lat')]
 
@@ -581,14 +581,14 @@ def _default_axis_physical_types(self):
             return ("custom:{}".format(self.unit[0].physical_type),)
 
     @property
-    def _world_axis_object_classes(self):
+    def world_axis_object_classes(self):
         return {'spectral': (
             coord.SpectralCoord,
             (),
             {'unit': self.unit[0]})}
 
     @property
-    def _world_axis_object_components(self):
+    def world_axis_object_components(self):
         return [('spectral', 0, 'value')]
 
 
@@ -633,8 +633,19 @@ def __init__(self, reference_frame, unit=None, axes_order=(0,),
     def _default_axis_physical_types(self):
         return ("time",)
 
+    def _convert_to_time(self, dt, *, unit, **kwargs):
+        if (not isinstance(dt, time.TimeDelta) and
+                isinstance(dt, time.Time) or
+                isinstance(self.reference_frame.value, np.ndarray)):
+            return time.Time(dt, **kwargs)
+
+        if not hasattr(dt, 'unit'):
+            dt = dt * unit
+
+        return self.reference_frame + dt
+
     @property
-    def _world_axis_object_classes(self):
+    def world_axis_object_classes(self):
         comp = (
             time.Time,
             (),
@@ -644,7 +655,7 @@ def _world_axis_object_classes(self):
         return {'temporal': comp}
 
     @property
-    def _world_axis_object_components(self):
+    def world_axis_object_components(self):
         if isinstance(self.reference_frame.value, np.ndarray):
             return [('temporal', 0, 'value')]
 
@@ -652,17 +663,6 @@ def offset_from_time_and_reference(time):
             return (time - self.reference_frame).sec
         return [('temporal', 0, offset_from_time_and_reference)]
 
-    def _convert_to_time(self, dt, *, unit, **kwargs):
-        if (not isinstance(dt, time.TimeDelta) and
-                isinstance(dt, time.Time) or
-                isinstance(self.reference_frame.value, np.ndarray)):
-            return time.Time(dt, **kwargs)
-
-        if not hasattr(dt, 'unit'):
-            dt = dt * unit
-
-        return self.reference_frame + dt
-
 
 class CompositeFrame(CoordinateFrame):
     """
@@ -699,10 +699,10 @@ def __init__(self, frames, name=None):
                              "axes_order should contain unique numbers, "
                              "got {}.".format(axes_order))
 
-        super(CompositeFrame, self).__init__(naxes, axes_type=axes_type,
-                                             axes_order=axes_order,
-                                             unit=unit, axes_names=axes_names,
-                                             name=name)
+        super().__init__(naxes, axes_type=axes_type,
+                         axes_order=axes_order,
+                         unit=unit, axes_names=axes_names,
+                         name=name)
         self._axis_physical_types = tuple(ph_type)
 
     @property
@@ -719,7 +719,7 @@ def _wao_classes_rename_map(self):
         for frame in self.frames:
             # ensure the frame is in the mapper
             mapper[frame]
-            for key in frame._world_axis_object_classes.keys():
+            for key in frame.world_axis_object_classes.keys():
                 if key in seen_names:
                     new_key = f"{key}{seen_names.count(key)}"
                     mapper[frame][key] = new_key
@@ -731,7 +731,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.world_axis_object_components:
                 comp = list(comp)
                 rename = mapper[frame].get(comp[0])
                 if rename:
@@ -743,14 +743,14 @@ def _wao_renamed_components_iter(self):
     def _wao_renamed_classes_iter(self):
         mapper = self._wao_classes_rename_map
         for frame in self.frames:
-            for key, value in frame._world_axis_object_classes.items():
+            for key, value in frame.world_axis_object_classes.items():
                 rename = mapper[frame].get(key)
                 if rename:
                     key = rename
                 yield key, value
 
     @property
-    def _world_axis_object_components(self):
+    def world_axis_object_components(self):
         """
         We need to generate the components respecting the axes_order.
         """
@@ -764,7 +764,7 @@ def _world_axis_object_components(self):
         return out
 
     @property
-    def _world_axis_object_classes(self):
+    def world_axis_object_classes(self):
         return dict(self._wao_renamed_classes_iter)
 
 
@@ -790,7 +790,7 @@ def _default_axis_physical_types(self):
         return ("phys.polarization.stokes",)
 
     @property
-    def _world_axis_object_classes(self):
+    def world_axis_object_classes(self):
         return {'stokes': (
             StokesCoord,
             (),
@@ -798,7 +798,7 @@ def _world_axis_object_classes(self):
         )}
 
     @property
-    def _world_axis_object_components(self):
+    def world_axis_object_components(self):
         return [('stokes', 0, 'value')]
 
 
diff --git a/gwcs/tests/test_api.py b/gwcs/tests/test_api.py
index 65b886f6..0c5b5c4e 100644
--- a/gwcs/tests/test_api.py
+++ b/gwcs/tests/test_api.py
@@ -491,12 +491,12 @@ def test_composite_many_base_frame():
     q_frame_2 = cf.CoordinateFrame(name='distance', axes_order=(1,), naxes=1, axes_type="SPATIAL", unit=(u.m,))
     frame = cf.CompositeFrame([q_frame_1, q_frame_2])
 
-    wao_classes = frame._world_axis_object_classes
+    wao_classes = frame.world_axis_object_classes
 
     assert len(wao_classes) == 2
     assert not set(wao_classes.keys()).difference({"SPATIAL", "SPATIAL1"})
 
-    wao_components = frame._world_axis_object_components
+    wao_components = frame.world_axis_object_components
 
     assert len(wao_components) == 2
     assert not {c[0] for c in wao_components}.difference({"SPATIAL", "SPATIAL1"})
diff --git a/gwcs/tests/test_coordinate_systems.py b/gwcs/tests/test_coordinate_systems.py
index 967657f8..88035b10 100644
--- a/gwcs/tests/test_coordinate_systems.py
+++ b/gwcs/tests/test_coordinate_systems.py
@@ -10,11 +10,12 @@
 from astropy.tests.helper import assert_quantity_allclose
 from astropy.modeling import models as m
 from astropy.wcs.wcsapi.fitswcs import CTYPE_TO_UCD1
-from astropy.coordinates import StokesCoord
+from astropy.coordinates import StokesCoord, SpectralCoord
 
 from .. import WCS
 from .. import coordinate_frames as cf
 
+from astropy.wcs.wcsapi.high_level_api import values_to_high_level_objects, high_level_objects_to_values
 import astropy
 astropy_version = astropy.__version__
 
@@ -33,7 +34,7 @@
 focal = cf.Frame2D(name='focal', axes_order=(0, 1), unit=(u.m, u.m))
 
 spec1 = cf.SpectralFrame(name='freq', unit=[u.Hz, ], axes_order=(2, ))
-spec2 = cf.SpectralFrame(name='wave', unit=[u.m, ], axes_order=(2, ), axes_names=('lambda', ))
+spec2 = cf.SpectralFrame(name='wave', unit=[u.m, ], axes_order=(2, ), axes_names=('lambda',))
 spec3 = cf.SpectralFrame(name='energy', unit=[u.J, ], axes_order=(2, ))
 spec4 = cf.SpectralFrame(name='pixel', unit=[u.pix, ], axes_order=(2, ))
 spec5 = cf.SpectralFrame(name='speed', unit=[u.m/u.s, ], axes_order=(2, ))
@@ -55,6 +56,19 @@
 inputs3 = [(xscalar, yscalar, xscalar), (xarr, yarr, xarr)]
 
 
+@pytest.fixture(autouse=True, scope="module")
+def serialized_classes():
+    """
+    In the rest of this test file we are passing the CoordinateFrame object to
+    astropy helper functions as if they were a low level WCS object.
+
+    This little patch means that this works.
+    """
+    cf.CoordinateFrame.serialized_classes = False
+    yield
+    del cf.CoordinateFrame.serialized_classes
+
+
 def test_units():
     assert(comp1.unit == (u.deg, u.deg, u.Hz))
     assert(comp2.unit == (u.m, u.m, u.m))
@@ -64,19 +78,34 @@ def test_units():
     assert(comp.unit == (u.deg, u.deg, u.Hz, u.m))
 
 
+# These two functions fake the old methods on CoordinateFrame to reduce the
+# amount of refactoring that needed doing in these tests.
+def coordinates(*inputs, frame):
+    results = values_to_high_level_objects(*inputs, low_level_wcs=frame)
+    if isinstance(results, list) and len(results) == 1:
+        return results[0]
+    return results
+
+
+def coordinate_to_quantity(*inputs, frame):
+    results = high_level_objects_to_values(*inputs, low_level_wcs=frame)
+    results = [r<<unit for r, unit in zip(results, frame.unit)]
+    return results
+
+
 @pytest.mark.parametrize('inputs', inputs2)
 def test_coordinates_spatial(inputs):
-    sky_coo = icrs.coordinates(*inputs)
+    sky_coo = coordinates(*inputs, frame=icrs)
     assert isinstance(sky_coo, coord.SkyCoord)
     assert_allclose((sky_coo.ra.value, sky_coo.dec.value), inputs)
-    focal_coo = focal.coordinates(*inputs)
+    focal_coo = coordinates(*inputs, frame=focal)
     assert_allclose([coo.value for coo in focal_coo], inputs)
     assert [coo.unit for coo in focal_coo] == [u.m, u.m]
 
 
 @pytest.mark.parametrize('inputs', inputs1)
 def test_coordinates_spectral(inputs):
-    wave = spec2.coordinates(inputs)
+    wave = coordinates(inputs, frame=spec2)
     assert_allclose(wave.value, inputs)
     assert wave.unit == 'meter'
     assert isinstance(wave, u.Quantity)
@@ -85,7 +114,7 @@ def test_coordinates_spectral(inputs):
 @pytest.mark.parametrize('inputs', inputs3)
 def test_coordinates_composite(inputs):
     frame = cf.CompositeFrame([icrs, spec2])
-    result = frame.coordinates(*inputs)
+    result = coordinates(*inputs, frame=frame)
     assert isinstance(result[0], coord.SkyCoord)
     assert_allclose((result[0].ra.value, result[0].dec.value), inputs[:2])
     assert_allclose(result[1].value, inputs[2])
@@ -96,7 +125,7 @@ def test_coordinates_composite_order():
     dist = cf.CoordinateFrame(name='distance', naxes=1,
                               axes_type=["SPATIAL"], unit=[u.m, ], axes_order=(1, ))
     frame = cf.CompositeFrame([time, dist])
-    result = frame.coordinates(0, 0)
+    result = coordinates(0, 0, frame=frame)
     assert result[0] == Time("2011-01-01T00:00:00")
     assert u.allclose(result[1], 0*u.m)
 
@@ -104,7 +133,8 @@ def test_coordinates_composite_order():
 def test_bare_baseframe():
     # This is a regression test for the following call:
     frame = cf.CoordinateFrame(1, "SPATIAL", (0,), unit=(u.km,))
-    assert u.allclose(frame.coordinate_to_quantity((1*u.m,)), 1*u.m)
+    quantity = coordinate_to_quantity(1*u.m, frame=frame)
+    assert u.allclose(quantity, 1*u.m)
 
     # Now also setup the same situation through the whole call stack to be safe.
     w = WCS(forward_transform=m.Tabular1D(points=np.arange(10)*u.pix,
@@ -158,55 +188,52 @@ def test_base_coordinate():
     assert frame.name == 'CoordinateFrame'
     frame = cf.CoordinateFrame(name="CustomFrame", naxes=2,
                                axes_type=("SPATIAL", "SPATIAL"),
-                               axes_order=(0, 1))
+                               axes_order=(0, 1),
+                               unit=(u.deg, u.arcsec))
     assert frame.name == 'CustomFrame'
     frame.name = "DeLorean"
     assert frame.name == 'DeLorean'
 
-    q1, q2 = frame.coordinate_to_quantity(12 * u.deg, 3 * u.arcsec)
+    q1, q2 = coordinate_to_quantity(12 * u.deg, 3 * u.arcsec, frame=frame)
     assert_quantity_allclose(q1, 12 * u.deg)
     assert_quantity_allclose(q2, 3 * u.arcsec)
 
-    q1, q2 = frame.coordinate_to_quantity((12 * u.deg, 3 * u.arcsec))
+    q1, q2 = coordinate_to_quantity(*(12 * u.deg, 3 * u.arcsec), frame=frame)
     assert_quantity_allclose(q1, 12 * u.deg)
     assert_quantity_allclose(q2, 3 * u.arcsec)
 
 
 def test_temporal_relative():
     t = cf.TemporalFrame(reference_frame=Time("2018-01-01T00:00:00"), unit=u.s)
-    assert t.coordinates(10) == Time("2018-01-01T00:00:00") + 10 * u.s
-    assert t.coordinates(10 * u.s) == Time("2018-01-01T00:00:00") + 10 * u.s
+    assert coordinates(10, frame=t) == Time("2018-01-01T00:00:00") + 10 * u.s
+    assert coordinates(10 * u.s, frame=t) == Time("2018-01-01T00:00:00") + 10 * u.s
 
-    a = t.coordinates((10, 20))
+    a = coordinates((10, 20), frame=t)
     assert a[0] == Time("2018-01-01T00:00:00") + 10 * u.s
     assert a[1] == Time("2018-01-01T00:00:00") + 20 * u.s
 
     t = cf.TemporalFrame(reference_frame=Time("2018-01-01T00:00:00"))
-    assert t.coordinates(10 * u.s) == Time("2018-01-01T00:00:00") + 10 * u.s
-    assert t.coordinates(TimeDelta(10, format='sec')) == Time("2018-01-01T00:00:00") + 10 * u.s
+    assert coordinates(10 * u.s, frame=t) == Time("2018-01-01T00:00:00") + 10 * u.s
+    assert coordinates(TimeDelta(10, format='sec'), frame=t) == Time("2018-01-01T00:00:00") + 10 * u.s
 
-    a = t.coordinates((10, 20) * u.s)
+    a = coordinates((10, 20) * u.s, frame=t)
     assert a[0] == Time("2018-01-01T00:00:00") + 10 * u.s
     assert a[1] == Time("2018-01-01T00:00:00") + 20 * u.s
 
 
-@pytest.mark.skipif(astropy_version<"4", reason="Requires astropy 4.0 or higher")
 def test_temporal_absolute():
     t = cf.TemporalFrame(reference_frame=Time([], format='isot'))
-    assert t.coordinates("2018-01-01T00:00:00") == Time("2018-01-01T00:00:00")
+    assert coordinates("2018-01-01T00:00:00", frame=t) == Time("2018-01-01T00:00:00")
 
-    a = t.coordinates(("2018-01-01T00:00:00", "2018-01-01T00:10:00"))
+    a = coordinates(("2018-01-01T00:00:00", "2018-01-01T00:10:00"), frame=t)
     assert a[0] == Time("2018-01-01T00:00:00")
     assert a[1] == Time("2018-01-01T00:10:00")
 
     t = cf.TemporalFrame(reference_frame=Time([], scale='tai', format='isot'))
-    assert t.coordinates("2018-01-01T00:00:00") == Time("2018-01-01T00:00:00", scale='tai')
+    assert coordinates("2018-01-01T00:00:00", frame=t) == Time("2018-01-01T00:00:00", scale='tai')
 
 
 @pytest.mark.parametrize('inp', [
-    (10 * u.deg, 20 * u.deg),
-    ((10 * u.deg, 20 * u.deg),),
-    (u.Quantity([10, 20], u.deg),),
     (coord.SkyCoord(10 * u.deg, 20 * u.deg, frame=coord.ICRS),),
     # This is the same as 10,20 in ICRS
     (coord.SkyCoord(119.26936774, -42.79039286, unit=u.deg, frame='galactic'),)
@@ -214,54 +241,40 @@ def test_temporal_absolute():
 def test_coordinate_to_quantity_celestial(inp):
     cel = cf.CelestialFrame(reference_frame=coord.ICRS(), axes_order=(0, 1))
 
-    lon, lat = cel.coordinate_to_quantity(*inp)
+    lon, lat = coordinate_to_quantity(*inp, frame=cel)
     assert_quantity_allclose(lon, 10 * u.deg)
     assert_quantity_allclose(lat, 20 * u.deg)
 
     with pytest.raises(ValueError):
-        cel.coordinate_to_quantity(10 * u.deg, 2 * u.deg, 3 * u.deg)
+        coordinate_to_quantity(10 * u.deg, 2 * u.deg, 3 * u.deg, frame=cel)
 
     with pytest.raises(ValueError):
-        cel.coordinate_to_quantity((1, 2))
+        coordinate_to_quantity((1, 2), frame=cel)
 
 
 @pytest.mark.parametrize('inp', [
-    (100,),
-    (100 * u.nm,),
-    (0.1 * u.um,),
+    (SpectralCoord(100 * u.nm),),
+    (SpectralCoord(0.1 * u.um),),
 ])
 def test_coordinate_to_quantity_spectral(inp):
     spec = cf.SpectralFrame(unit=u.nm, axes_order=(1, ))
-    wav = spec.coordinate_to_quantity(*inp)
+    wav = coordinate_to_quantity(*inp, frame=spec)
     assert_quantity_allclose(wav, 100 * u.nm)
 
 
 @pytest.mark.parametrize('inp', [
     (Time("2011-01-01T00:00:10"),),
-    (10 * u.s,)
 ])
-@pytest.mark.skipif(astropy_version<"4", reason="Requires astropy 4.0 or higher.")
 def test_coordinate_to_quantity_temporal(inp):
     temp = cf.TemporalFrame(reference_frame=Time("2011-01-01T00:00:00"), unit=u.s)
 
-    t = temp.coordinate_to_quantity(*inp)
+    t = coordinate_to_quantity(*inp, frame=temp)
 
     assert_quantity_allclose(t, 10 * u.s)
 
-    temp2 = cf.TemporalFrame(reference_frame=Time([], format='isot'), unit=u.s)
-
-    tt = Time("2011-01-01T00:00:00")
-    t = temp2.coordinate_to_quantity(tt)
-
-    assert t is tt
-
 
 @pytest.mark.parametrize('inp', [
-    (211 * u.AA, 0 * u.s, 0 * u.arcsec, 0 * u.arcsec),
-    (211 * u.AA, 0 * u.s, (0 * u.arcsec, 0 * u.arcsec)),
-    (211 * u.AA, 0 * u.s, (0, 0) * u.arcsec),
-    (211 * u.AA, Time("2011-01-01T00:00:00"), (0, 0) * u.arcsec),
-    (211 * u.AA, Time("2011-01-01T00:00:00"), coord.SkyCoord(0, 0, unit=u.arcsec)),
+    (SpectralCoord(211 * u.AA), Time("2011-01-01T00:00:00"), coord.SkyCoord(0, 0, unit=u.arcsec)),
 ])
 def test_coordinate_to_quantity_composite(inp):
     # Composite
@@ -272,29 +285,46 @@ def test_coordinate_to_quantity_composite(inp):
 
     comp = cf.CompositeFrame([wave_frame, time_frame, sky_frame])
 
-    coords = comp.coordinate_to_quantity(*inp)
+    coords = coordinate_to_quantity(*inp, frame=comp)
 
     expected = (211 * u.AA, 0 * u.s, 0 * u.arcsec, 0 * u.arcsec)
     for output, exp in zip(coords, expected):
         assert_quantity_allclose(output, exp)
 
 
+def test_coordinate_to_quantity_composite_split():
+    inp = (
+        SpectralCoord(211 * u.AA),
+        coord.SkyCoord(0, 0, unit=u.arcsec),
+        Time("2011-01-01T00:00:00"),
+    )
+
+    # Composite
+    wave_frame = cf.SpectralFrame(axes_order=(1, ), unit=u.AA)
+    sky_frame = cf.CelestialFrame(axes_order=(2, 0), reference_frame=coord.ICRS())
+    time_frame = cf.TemporalFrame(
+        axes_order=(3,), unit=u.s, reference_frame=Time("2011-01-01T00:00:00"))
+
+    comp = cf.CompositeFrame([wave_frame, sky_frame, time_frame])
+
+    coords = coordinate_to_quantity(*inp, frame=comp)
+
+    expected = (0 * u.arcsec, 211 * u.AA, 0 * u.arcsec, 0 * u.s)
+    for output, exp in zip(coords, expected):
+        assert_quantity_allclose(output, exp)
+
+
 def test_stokes_frame():
     sf = cf.StokesFrame()
 
-    assert sf.coordinates(1) == 'I'
-    assert sf.coordinates(1 * u.pix) == 'I'
-    assert sf.coordinate_to_quantity(StokesCoord('I')) == 1 * u.one
-    assert sf.coordinate_to_quantity(1) == 1
+    assert coordinates(1, frame=sf) == 'I'
+    assert coordinates(1 * u.one, frame=sf) == 'I'
+    assert coordinate_to_quantity(StokesCoord('I'), frame=sf) == 1 * u.one
+    assert coordinate_to_quantity(StokesCoord(1), frame=sf) == 1 * u.one
 
 
-@pytest.mark.parametrize('inp', [
-    (211 * u.AA, 0 * u.s, 0 * u.one, 0 * u.one),
-    (211 * u.AA, 0 * u.s, (0 * u.one, 0 * u.one)),
-    (211 * u.AA, 0 * u.s, (0, 0) * u.one),
-    (211 * u.AA, Time("2011-01-01T00:00:00"), (0, 0) * u.one)
-])
-def test_coordinate_to_quantity_frame2d_composite(inp):
+def test_coordinate_to_quantity_frame2d_composite():
+    inp = (SpectralCoord(211 * u.AA), Time("2011-01-01T00:00:00"), 0 * u.one, 0 * u.one)
     wave_frame = cf.SpectralFrame(axes_order=(0, ), unit=u.AA)
     time_frame = cf.TemporalFrame(
         axes_order=(1, ), unit=u.s, reference_frame=Time("2011-01-01T00:00:00"))
@@ -303,7 +333,7 @@ def test_coordinate_to_quantity_frame2d_composite(inp):
 
     comp = cf.CompositeFrame([wave_frame, time_frame, frame2d])
 
-    coords = comp.coordinate_to_quantity(*inp)
+    coords = coordinate_to_quantity(*inp, frame=comp)
 
     expected = (211 * u.AA, 0 * u.s, 0 * u.one, 0 * u.one)
     for output, exp in zip(coords, expected):
@@ -312,31 +342,24 @@ def test_coordinate_to_quantity_frame2d_composite(inp):
 
 def test_coordinate_to_quantity_frame_2d():
     frame = cf.Frame2D(unit=(u.one, u.arcsec))
-    inp = (1, 2)
-    expected = (1 * u.one, 2 * u.arcsec)
-    result = frame.coordinate_to_quantity(*inp)
-    for output, exp in zip(result, expected):
-        assert_quantity_allclose(output, exp)
-
-    inp = (1 * u.one, 2)
+    inp = (1 * u.one, 2 * u.arcsec)
     expected = (1 * u.one, 2 * u.arcsec)
-    result = frame.coordinate_to_quantity(*inp)
+    result = coordinate_to_quantity(*inp, frame=frame)
     for output, exp in zip(result, expected):
         assert_quantity_allclose(output, exp)
 
 
-@pytest.mark.skipif(astropy_version<"4", reason="Requires astropy 4.0 or higher.")
 def test_coordinate_to_quantity_error():
     frame = cf.Frame2D(unit=(u.one, u.arcsec))
     with pytest.raises(ValueError):
-        frame.coordinate_to_quantity(1)
+        coordinate_to_quantity(1, frame=frame)
 
     with pytest.raises(ValueError):
-        comp1.coordinate_to_quantity((1, 1), 2)
+        coordinate_to_quantity((1, 1), 2, frame=frame)
 
     frame = cf.TemporalFrame(reference_frame=Time([], format='isot'), unit=u.s)
     with pytest.raises(ValueError):
-        frame.coordinate_to_quantity(1)
+        coordinate_to_quantity(1, frame=frame)
 
 
 def test_axis_physical_type():
@@ -406,7 +429,7 @@ def test_base_frame():
     assert frame.naxes == 1
     assert frame.axes_names == ("x",)
 
-    frame.coordinate_to_quantity(1, 2)
+    coordinate_to_quantity(1*u.one, frame=frame)
 
 
 def test_ucd1_to_ctype_not_out_of_sync(caplog):
diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py
index c2126acc..e7cd7dae 100644
--- a/gwcs/tests/test_wcs.py
+++ b/gwcs/tests/test_wcs.py
@@ -1303,3 +1303,32 @@ def test_spatial_spectral_stokes():
     assert_allclose(gw_sky.data.lat, aw_sky.data.lat)
     assert_allclose(gw_spec.value, aw_spec.value)
     assert_allclose(gw_stokes.value, aw_stokes.value)
+
+
+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
+    compound = models.Linear1D(intercept=0*u.nm, slope=10*u.nm/u.pix) & spatial
+    forward = models.Mapping((1, 2, 0)) | compound | models.Mapping((2, 0, 1))
+
+    # Setup the output frame
+    celestial_frame = cf.CelestialFrame(axes_order=(2, 0), unit=(u.arcsec, u.arcsec),
+                                        reference_frame=coord.ICRS())
+    spectral_frame = cf.SpectralFrame(axes_order=(1,), unit=u.nm)
+    output_frame = cf.CompositeFrame([spectral_frame, celestial_frame])
+
+    input_frame = cf.CoordinateFrame(3, ["PIXEL"]*3,
+                                     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]
+    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))
+
+    world_obj = iwcs.pixel_to_world(*input_pixel)
+    assert isinstance(world_obj[0], coord.SkyCoord)
+    assert isinstance(world_obj[1], coord.SpectralCoord)
+
+    obj_pixel = iwcs.world_to_pixel(*world_obj)
+    assert_allclose(obj_pixel, u.Quantity(input_pixel).to_value(u.pix))