From f434dad84cef3ce13de69b7f8a376ba3729d01b4 Mon Sep 17 00:00:00 2001
From: Nadia Dencheva <nadia.dencheva@gmail.com>
Date: Tue, 2 Sep 2025 07:52:17 -0400
Subject: [PATCH 01/15] Updates for new legacy API.

fix rebase
---
 gwcs/api.py                       |  14 +-
 gwcs/coordinate_frames/_base.py   |   7 +-
 gwcs/tests/test_api.py            |  22 +--
 gwcs/tests/test_api_consistent.py | 223 +++++++++++++++++++++++++
 gwcs/tests/test_wcs.py            |  11 +-
 gwcs/wcs/_wcs.py                  | 269 +++++++++++++++++++-----------
 6 files changed, 419 insertions(+), 127 deletions(-)
 create mode 100644 gwcs/tests/test_api_consistent.py

diff --git a/gwcs/api.py b/gwcs/api.py
index 9ecd8521..53d38f54 100644
--- a/gwcs/api.py
+++ b/gwcs/api.py
@@ -66,14 +66,13 @@ def world_axis_units(self):
         return tuple(unit.to_string(format="vounit") for unit in self.output_frame.unit)
 
     def _remove_quantity_output(self, result, frame):
-        if self.forward_transform.uses_quantity:
-            if frame.naxes == 1:
-                result = [result]
+        if frame.naxes == 1:
+            result = [result]
 
-            result = tuple(
-                r.to_value(unit) if isinstance(r, u.Quantity) else r
-                for r, unit in zip(result, frame.unit, strict=False)
-            )
+        result = tuple(
+            r.to_value(unit) if isinstance(r, u.Quantity) else r
+            for r, unit in zip(result, frame.unit, strict=False)
+        )
 
         # If we only have one output axes, we shouldn't return a tuple.
         if self.output_frame.naxes == 1 and isinstance(result, tuple):
@@ -94,7 +93,6 @@ def pixel_to_world_values(self, *pixel_arrays):
         is the vertical coordinate.
         """
         result = self(*pixel_arrays)
-
         return self._remove_quantity_output(result, self.output_frame)
 
     def array_index_to_world_values(self, *index_arrays):
diff --git a/gwcs/coordinate_frames/_base.py b/gwcs/coordinate_frames/_base.py
index eda53c6f..d9873cc6 100644
--- a/gwcs/coordinate_frames/_base.py
+++ b/gwcs/coordinate_frames/_base.py
@@ -128,6 +128,8 @@ def add_units(self, arrays: u.Quantity | np.ndarray | float) -> tuple[u.Quantity
         """
         Add units to the arrays
         """
+        if self.naxes == 1 and np.isscalar(arrays):
+            return u.Quantity(arrays, self.unit[0])
         return tuple(
             u.Quantity(array, unit=unit)
             for array, unit in zip(arrays, self.unit, strict=True)
@@ -139,10 +141,11 @@ def remove_units(
         """
         Remove units from the input arrays
         """
-        if self.naxes == 1:
+        if self.naxes == 1 and (np.isscalar(arrays) or isinstance(arrays, u.Quantity)):
             arrays = (arrays,)
 
-        return tuple(
+        result = tuple(
             array.to_value(unit) if isinstance(array, u.Quantity) else array
             for array, unit in zip(arrays, self.unit, strict=True)
         )
+        return result
diff --git a/gwcs/tests/test_api.py b/gwcs/tests/test_api.py
index 528def79..528462de 100644
--- a/gwcs/tests/test_api.py
+++ b/gwcs/tests/test_api.py
@@ -309,6 +309,7 @@ def test_high_level_wrapper(wcsobj, request):
     wc1 = hlvl.pixel_to_world(*pixel_input)
     wc2 = wcsobj(*pixel_input)
     results = wcsobj._remove_units_input(wc2, wcsobj.output_frame)
+
     wc2 = values_to_high_level_objects(*results, low_level_wcs=wcsobj)
     if len(wc2) == 1:
         wc2 = wc2[0]
@@ -325,19 +326,9 @@ def test_high_level_wrapper(wcsobj, request):
         wc1 = (wc1,)
 
     pix_out1 = hlvl.world_to_pixel(*wc1)
-    pix_out2 = wcsobj.invert(*wc1)
-
-    if not isinstance(pix_out2, list | tuple):
-        pix_out2 = (pix_out2,)
-
-    if wcsobj.forward_transform.uses_quantity:
-        pix_out2 = tuple(
-            p.to_value(unit)
-            for p, unit in zip(pix_out2, wcsobj.input_frame.unit, strict=False)
-        )
-
     np.testing.assert_allclose(pix_out1, pixel_input)
-    np.testing.assert_allclose(pix_out2, pixel_input)
+    with pytest.raises(TypeError) as e:
+        pix_out2 = wcsobj.invert(*wc1)
 
 
 def test_stokes_wrapper(gwcs_stokes_lookup):
@@ -407,7 +398,8 @@ def test_pixel_bounds(wcsobj):
 
     wcsobj.bounding_box = ((-0.5, 2039.5), (-0.5, 1019.5))
     assert_array_equal(wcsobj.pixel_bounds, wcsobj.bounding_box)
-
+    # Reset the bounding box or this will affect other tests
+    wcsobj.bounding_box = None
 
 @wcs_objs
 def test_axis_correlation_matrix(wcsobj):
@@ -598,8 +590,8 @@ def test_coordinate_frame_api():
     pixel = wcs.world_to_pixel(world)
     assert isinstance(pixel, float)
 
-    pixel2 = wcs.invert(world)
-    assert u.allclose(pixel2, 0 * u.pix)
+    with pytest.raises(TypeError):
+        pixel2 = wcs.invert(world)
 
 
 def test_world_axis_object_components_units(gwcs_3d_identity_units):
diff --git a/gwcs/tests/test_api_consistent.py b/gwcs/tests/test_api_consistent.py
new file mode 100644
index 00000000..f91e8a14
--- /dev/null
+++ b/gwcs/tests/test_api_consistent.py
@@ -0,0 +1,223 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+"""
+Test the API is consistent with units and quantities and follows the rules below.
+
+WCS functions which considered for this work are part of the legacy API:
+wcs(x, y)
+wcs.invert(ra, dec)
+wcs.forward_transform(x,y), wcs.backward_transform() and wcs.get_transform(f1, f2)
+wcs.numerical_inverse(ra, dec) - does not support units
+
+Rules:
+
+
+1. Neither transforms nor inputs support units -> the output is clearly numerical for all functions above
+2. Transforms support units but inputs do not -> return quantities assuming the units of the coordinate frame
+  - This should work for the wcs methods (wcs(x,y) and wcs.invert
+  - The methods using transforms should follow modeling rules and will require units
+    on the input and raise an exception if not
+3. Both transforms and inputs support units -> return quantities
+  - Wcs methods return quantities
+  - Transforms work and return quantities
+4. Transforms do not support units but inputs are quantities -> raise an error
+
+"""
+import numbers
+import numpy as np
+from numpy.testing import assert_array_equal, assert_allclose
+
+from astropy import units as u
+from astropy.tests.helper import assert_quantity_allclose
+
+import pytest
+from .conftest import gwcs_with_pipeline_celestial
+
+x = 1
+y = 2
+xq = [1, 1] * u.pix
+yq = 2 * u.pix
+
+
+def is_numerical(args):
+    if isinstance(args, numbers.Number):
+        return True
+    return all([isinstance(arg, numbers.Number) or type(arg) == np.ndarray for arg in args])
+
+
+def is_quantity(args):
+    return all([isinstance(arg, u.Quantity) for arg in args])
+
+
+@pytest.fixture
+def wcsobj(request):
+    return request.getfixturevalue(request.param)
+
+wno_unit_1d = ["gwcs_1d_freq", "gwcs_1d_spectral",]
+
+wno_unit_nd = ["gwcs_2d_shift_scale", "gwcs_3d_spatial_wave", "gwcs_2d_spatial_shift", "gwcs_2d_spatial_reordered",
+    "gwcs_3d_spatial_wave", "gwcs_simple_imaging", "gwcs_3spectral_orders",
+    "gwcs_3d_galactic_spectral", "gwcs_spec_cel_time_4d", "gwcs_romanisim", ]
+
+# "gwcs_7d_complex_mapping" errors in astropy - fix
+# "gwcs_2d_quantity_shift" errors when inputs are quantities. Need to confirm if Qs are HLO
+
+w_unit_1d = ["gwcs_stokes_lookup", "gwcs_1d_freq_quantity"]
+
+w_unit_nd = ["gwcs_2d_shift_scale_quantity", "gwcs_3d_identity_units", "gwcs_3d_identity_units",
+    "gwcs_4d_identity_units", "gwcs_simple_imaging_units", "gwcs_with_pipeline_celestial", ]
+
+w_transform_test = ["gwcs_1d_freq_quantity", "gwcs_2d_quantity_shift"]
+
+wcs_no_unit_1d = pytest.mark.parametrize(("wcsobj"), wno_unit_1d, indirect=True)
+wcs_no_unit_nd = pytest.mark.parametrize(("wcsobj"), wno_unit_nd, indirect=True)
+wcs_with_unit_1d = pytest.mark.parametrize(("wcsobj"), w_unit_1d, indirect=True)
+wcs_with_unit_nd = pytest.mark.parametrize(("wcsobj"), w_unit_nd, indirect=True)
+
+
+@wcs_no_unit_1d
+def test_trnou_inpnou_1d(wcsobj):
+    """ Transforms do not support units. Inputs are numbers."""
+    assert not wcsobj.forward_transform.uses_quantity
+
+    # the case of a scalar input
+    x = 1
+    bbox = wcsobj.bounding_box
+    if bbox is not None:
+        x = np.mean(bbox.bounding_box())
+
+    result_num = wcsobj(x)
+    assert np.isscalar(result_num)
+
+    assert_allclose(wcsobj.invert(result_num), x)
+
+    xq = x * wcsobj.input_frame.unit[0]
+    result = wcsobj(xq)
+    assert_quantity_allclose(result, result_num * wcsobj.output_frame.unit[0])
+
+
+@wcs_no_unit_nd
+def test_no_units_nd(wcsobj):
+    assert not wcsobj.forward_transform.uses_quantity
+
+    n_inputs = wcsobj.input_frame.naxes
+
+    inp = [1] * n_inputs
+    bbox = wcsobj.bounding_box
+    if bbox is not None:
+        bb = bbox.bounding_box()
+        inp = [np.mean(interval) for interval in bb]
+    # Inputs are numbers
+    result = wcsobj(*inp)
+    assert is_numerical(result)
+    if np.isscalar(result):
+        result = [result]
+    inp_new = wcsobj.invert(*result)
+    _ = [assert_allclose(i, j) for i, j in zip(inp_new, inp)]
+
+    # Inputs are quantities; return quantities (except for pixels?)
+    inpq = [coo * un for coo, un in zip(inp, wcsobj.input_frame.unit)]
+    result = wcsobj(*inpq)
+    assert is_quantity(result)
+    inp_new = wcsobj.invert(*result)
+    _ = [assert_allclose(i, j) for i, j in zip(inp_new, inpq)]
+
+    # input is HLO - raise an Error
+    sky = wcsobj.pixel_to_world(*inp)
+    with pytest.raises(TypeError) as e:
+        wcsobj.invert(*sky)
+
+
+@wcs_with_unit_1d
+def test_with_units_1d(wcsobj):
+    """ Transform do not support units."""
+    assert wcsobj.forward_transform.uses_quantity
+
+    # the case of a scalar input
+    x = 1 * wcsobj.input_frame.unit[0]
+
+    result = wcsobj(x)
+    assert isinstance(result, u.Quantity)
+    assert_allclose(wcsobj.invert(result), x)
+
+    x = 1
+    result = wcsobj(x)
+    assert np.isscalar(result)
+    assert_allclose(wcsobj.invert(result), x)
+
+
+@wcs_with_unit_nd
+def test_transform_with_units(wcsobj):
+    """ Transforms support units."""
+    assert wcsobj.forward_transform.uses_quantity
+
+    n_inputs = wcsobj.input_frame.naxes
+    xx = [x] * n_inputs
+
+    # input is numerical; return numbers
+    result_num = wcsobj(*xx)
+    assert is_numerical(result_num)
+
+    inp = wcsobj.invert(*result_num)
+    assert is_numerical(inp)
+
+    # input is quantities; return quantities
+    xxq = [1 * u.pix] * n_inputs
+    result = wcsobj(*xxq)
+    assert all([type(res)==u.Quantity for res in result])
+    assert_allclose([r.value for r in result], result_num)
+
+    # input is HLO; raise an error
+    sky = wcsobj.pixel_to_world(*xxq)
+    with pytest.raises(TypeError) as e:
+        wcsobj.invert(*sky)
+
+
+@wcs_no_unit_1d
+def test_add_units(wcsobj):
+    if wcsobj.input_frame.naxes == 1:
+        assert wcsobj._add_units_input((1,), wcsobj.input_frame) == 1 * wcsobj.input_frame.unit[0]
+        assert_allclose(
+            wcsobj._add_units_input(([1, 1],), wcsobj.input_frame),
+            ([1, 1] * wcsobj.input_frame.unit[0],))
+    elif wcsobj.input_frame.naxes == 2:
+        assert_quantity_allclose(
+            wcsobj._add_units_input((1, 1), wcsobj.input_frame),
+            (1*u.pix, 1*u.pix))
+        assert_quantity_allclose(
+            wcsobj._add_units_input(([1, 1], [1, 1]), wcsobj.input_frame),
+            ([1, 1]*u.pix, [1, 1]*u.pix))
+
+
+@wcs_with_unit_1d
+def test_remove_units(wcsobj):
+    if wcsobj.input_frame.naxes == 1:
+        unit = wcsobj.input_frame.unit[0]
+        assert wcsobj._remove_units_input(1 * unit, wcsobj.input_frame) == (1,)
+        assert_allclose(
+            wcsobj._remove_units_input(([1, 1] * unit,), wcsobj.input_frame),
+            ([1, 1],))
+    elif wcsobj.input_frame.naxes == 2:
+        assert_quantity_allclose(
+            wcsobj._remove_units_input((1*u.pix, 1*u.pix), wcsobj.input_frame),
+            (1, 1)
+            )
+        assert_quantity_allclose(
+            wcsobj._remove_units_input(([1, 1]*u.pix, [1, 1]*u.pix), wcsobj.input_frame),
+            ([1, 1], [1, 1]))
+
+
+def test_transform_multistage_wcs(gwcs_with_pipeline_celestial):
+    wcsobj = gwcs_with_pipeline_celestial
+    frames = wcsobj.available_frames
+    result = wcsobj.transform(frames[0], frames[-1], 1*u.pix, 1*u.pix)
+    assert is_quantity(result)
+    assert_allclose([r.value for r in result], wcsobj(1, 1))
+    final_result = wcsobj.transform(frames[0], frames[-1], 1*u.pix, 1*u.pix)
+    assert is_quantity(final_result)
+    assert_allclose([r.value for r in final_result], wcsobj(1, 1))
+    interm_result = wcsobj.transform(frames[0], frames[1], 1*u.pix, 1*u.pix)
+    assert is_quantity(interm_result)
+    tr = wcsobj.get_transform(frames[0], frames[1])
+    assert_quantity_allclose(interm_result, tr(1*u.pix, 1*u.pix))
+    ninterm_result = wcsobj.transform(frames[0], frames[1], 1, 1)
+    assert_allclose([r.value for r in interm_result], ninterm_result)
diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py
index 7d0b221e..377d031a 100644
--- a/gwcs/tests/test_wcs.py
+++ b/gwcs/tests/test_wcs.py
@@ -516,10 +516,7 @@ def test_bounding_box_eval():
     pipeline = [
         (
             cf.CoordinateFrame(
-                naxes=3,
-                axes_type=("PIXEL", "PIXEL", "PIXEL"),
-                axes_order=(0, 1, 2),
-                name="detector",
+                naxes=3, axes_type=("PIXEL", "PIXEL", "PIXEL",), axes_order=(0, 1, 2), name="detector"
             ),
             trans3,
         ),
@@ -1703,8 +1700,8 @@ def test_quantities_in_pipeline_forward(gwcs_with_pipeline_celestial):
 
     output_world = iwcs(*input_pixel)
 
-    assert output_world[0].unit == u.deg
-    assert output_world[1].unit == u.deg
+    assert output_world[0].unit == u.arcsec
+    assert output_world[1].unit == u.arcsec
     assert u.allclose(output_world[0], 20 * u.arcsec + 1 * u.deg)
     assert u.allclose(output_world[1], 15 * u.deg + 2 * u.deg)
 
@@ -1883,6 +1880,8 @@ def test_parameterless_transform():
     assert gwcs(1 * u.pix, 1 * u.pix) == (1 * u.pix, 1 * u.pix)
 
     assert gwcs.invert(1, 1) == (1, 1)
+    # Strictly speaking it's correct that this fails Because
+    # for this setup the HLO are Quantities
     assert gwcs.invert(1 * u.pix, 1 * u.pix) == (1, 1)
 
 
diff --git a/gwcs/wcs/_wcs.py b/gwcs/wcs/_wcs.py
index 786e9651..1c34bef2 100644
--- a/gwcs/wcs/_wcs.py
+++ b/gwcs/wcs/_wcs.py
@@ -29,6 +29,7 @@
     CelestialFrame,
     CompositeFrame,
     CoordinateFrame,
+    EmptyFrame,
     get_ctype_from_ucd,
 )
 from gwcs.utils import _compute_lon_pole, is_high_level, to_index
@@ -157,41 +158,135 @@ def __call__(
         if transform is None:
             msg = "Transform is not defined."
             raise NotImplementedError(msg)
+        # move this to an evaluate_funciton, called by forward, transform and invert
+        # input_is_quantity = any(isinstance(a, u.Quantity) for a in args)
+        # transform_uses_quantity = not (transform is None or not transform.uses_quantity)
+
+        input_is_quantity, transform_uses_quantity = self._units_are_present(args, transform)
         args = self._make_input_units_consistent(
+            transform,
             *args,
-            transform=transform,
-            from_frame=self.input_frame,
-            to_frame=self.output_frame,
+            frame=self.input_frame,
+            input_is_quantity=input_is_quantity,
+            transform_uses_quantity=transform_uses_quantity,
         )
 
-        return self._call_forward(
-            *args, with_bounding_box=with_bounding_box, fill_value=fill_value, **kwargs
-        )
+        result = transform(*args, with_bounding_box=with_bounding_box, fill_value=fill_value, **kwargs)
+        if self.output_frame is not None:
+            if self.output_frame.naxes == 1:
+                result = (result,)
+            result = self._make_output_units_consistent(
+                transform,
+                *result,
+                frame=self.output_frame,
+                input_is_quantity=input_is_quantity,
+                transform_uses_quantity=transform_uses_quantity,
+            )
+        if self.output_frame is not None and self.output_frame.naxes == 1:
+            return result[0]
+        return result
+
+    def _units_are_present(self, args, transform):
+        """
+        Determing if the inputs to a transform are quantities and the transform
+        supports units.
+
+        Parameters
+        ----------
+        args : a tuple of scalars or ndarray-like objects
+            Inputs to a transform.
+        transform : `~astropy.modeling.Model`
+            Transform to be evaluated.
+
+        Returns
+        -------
+        input_is_quantity, transform_uses_quantity : bool
+
+        """
+        # Validate that the input type matches what the transform expects
+        input_is_quantity = any(isinstance(a, u.Quantity) for a in args)
+        transform_uses_quantity = not (transform is None or not transform.uses_quantity)
+        return input_is_quantity, transform_uses_quantity
+
 
     def _make_input_units_consistent(
         self,
-        *args,
         transform,
-        from_frame: CoordinateFrame | None = None,
-        to_frame: CoordinateFrame | None = None,
+        *args,
+        frame: CoordinateFrame | None = None,
+        input_is_quantity=False,
+        transform_uses_quantity=False,
         **kwargs,
     ):
         """
         Adds or removes units from the arguments as needed so that the transform
         can be successfully evaluated.
         """
-        # Validate that the input type matches what the transform expects
-        input_is_quantity = any(isinstance(a, u.Quantity) for a in args)
-        transform_uses_quantity = not (transform is None or not transform.uses_quantity)
-        if (
+        # # Validate that the input type matches what the transform expects
+        # input_is_quantity = any(isinstance(a, u.Quantity) for a in args)
+        # transform_uses_quantity = not (transform is None or not transform.uses_quantity)
+        if not input_is_quantity and not transform_uses_quantity:
+            return args
+        elif input_is_quantity and transform_uses_quantity:
+            return args
+        elif (
             not input_is_quantity
-            and transform_uses_quantity
-            and transform.parameters.size
+            and (transform_uses_quantity
+            or transform.parameters.size) # possibly remove this, check is in _units_are_present
         ):
-            return self._add_units_input(args, from_frame)
-        if not transform_uses_quantity and input_is_quantity:
-            return self._remove_units_input(args, from_frame)
-        return args
+            return self._add_units_input(args, frame)
+        elif not transform_uses_quantity and input_is_quantity:
+            return self._remove_units_input(args, frame)
+        #return args
+
+    def _make_output_units_consistent(
+        self,
+        transform,
+        *args,
+        frame: CoordinateFrame | None = None,
+        input_is_quantity=False,
+        transform_uses_quantity=False,
+        **kwargs,
+    ):
+        """
+        Adds or removes units from the arguments as needed so that the type of the output
+        matches the input.
+        """
+        if not input_is_quantity and not transform_uses_quantity:
+            return args
+
+        elif input_is_quantity and transform_uses_quantity:
+            # make sure the output is returned in the units of the output frame
+            return self._add_units_input(args, frame)
+        elif (
+            not input_is_quantity
+            and (transform_uses_quantity
+            or transform.parameters.size)
+        ):
+            result = self._remove_units_input(args, frame)
+        elif not transform_uses_quantity and input_is_quantity:
+            result = self._add_units_input(args, frame)
+        return result
+
+    def _get_transform(
+        self,
+        *args,
+        from_frame: CoordinateFrame | None = None,
+        to_frame: CoordinateFrame | None = None,
+        **kwargs,
+    ):
+        """
+        Executes the forward transform, but values only.
+        """
+        if from_frame is None and to_frame is None:
+            transform = self.forward_transform
+        else:
+            transform = self.get_transform(from_frame, to_frame)
+
+        if transform is None:
+            msg = "WCS.forward_transform is not implemented."
+            raise NotImplementedError(msg)
+        return transform
 
     def _evaluate_transform(
         self,
@@ -261,41 +356,6 @@ def _transform(*args):
 
             raise
 
-    def _call_forward(
-        self,
-        *args,
-        from_frame: CoordinateFrame | None = None,
-        to_frame: CoordinateFrame | None = None,
-        with_bounding_box: bool = True,
-        fill_value: float | np.number = np.nan,
-        **kwargs,
-    ):
-        """
-        Executes the forward transform, but values only.
-        """
-        if from_frame is None and to_frame is None:
-            transform = self.forward_transform
-        else:
-            transform = self.get_transform(from_frame, to_frame)
-        if from_frame is None:
-            from_frame = self.input_frame
-        if to_frame is None:
-            to_frame = self.output_frame
-
-        if transform is None:
-            msg = "WCS.forward_transform is not implemented."
-            raise NotImplementedError(msg)
-
-        return self._evaluate_transform(
-            transform,
-            from_frame,
-            to_frame,
-            *args,
-            with_bounding_box=with_bounding_box,
-            fill_value=fill_value,
-            **kwargs,
-        )
-
     def in_image(self, *args, **kwargs):
         """
         This method tests if one or more of the input world coordinates are
@@ -380,46 +440,28 @@ def invert(
             transform = self.backward_transform
         except NotImplementedError:
             transform = None
+            # TODO raise error?
         if is_high_level(*args, low_level_wcs=self):
-            args = high_level_objects_to_values(*args, low_level_wcs=self)
-        args = self._make_input_units_consistent(
-            *args,
-            transform=transform,
-            from_frame=self.output_frame,
-            to_frame=self.input_frame,
-        )
-
-        return self._call_backward(
-            *args, with_bounding_box=with_bounding_box, fill_value=fill_value, **kwargs
-        )
-
-    def _call_backward(
-        self,
-        *args,
-        with_bounding_box: bool = True,
-        fill_value: float | np.number = np.nan,
-        **kwargs,
-    ):
-        try:
-            transform = self.backward_transform
-        except NotImplementedError:
-            transform = None
+            # args = high_level_objects_to_values(*args, low_level_wcs=self)
+            message = "High Level objects are not supported with the native API. \
+                       Please use the `world_to_pixel` method."
+            raise TypeError(message)
 
         if with_bounding_box and self.bounding_box is not None:
             args = self.outside_footprint(args)
 
+        input_is_quantity, transform_uses_quantity = self._units_are_present(args, transform)
+
+        args = self._make_input_units_consistent(
+            transform,
+            *args,
+            frame=self.output_frame,
+            input_is_quantity=input_is_quantity,
+            transform_uses_quantity=transform_uses_quantity,
+        )
         if transform is not None:
-            # remove iterative inverse-specific keyword arguments:
             akwargs = {k: v for k, v in kwargs.items() if k not in _ITER_INV_KWARGS}
-            result = self._evaluate_transform(
-                transform,
-                self.output_frame,
-                self.input_frame,
-                *args,
-                with_bounding_box=with_bounding_box,
-                fill_value=fill_value,
-                **akwargs,
-            )
+            result = transform(*args, with_bounding_box=with_bounding_box, fill_value=fill_value, **akwargs)
         else:
             # Always strip units for numerical inverse
             args = self._remove_units_input(args, self.output_frame)
@@ -430,10 +472,20 @@ def _call_backward(
                 **kwargs,
             )
 
-        # deal with values outside the bounding box
         if with_bounding_box and self.bounding_box is not None:
             result = self.out_of_bounds(result, fill_value=fill_value)
 
+        if self.input_frame.naxes == 1:
+            result = (result,)
+        result = self._make_output_units_consistent(
+            transform,
+            *result,
+            frame=self.input_frame,
+            input_is_quantity=input_is_quantity,
+            transform_uses_quantity=transform_uses_quantity,
+        )
+        if self.input_frame.naxes == 1:
+            return result[0]
         return result
 
     def outside_footprint(self, world_arrays):
@@ -1175,17 +1227,42 @@ def transform(
         to_step = self._get_step(to_frame)
         transform = self.get_transform(from_step.step.frame, to_step.step.frame)
 
-        if not transform.uses_quantity and is_high_level(
-            *args, low_level_wcs=from_step.step.frame
-        ):
-            args = high_level_objects_to_values(
-                *args, low_level_wcs=from_step.step.frame
-            )
+        with_bounding_box = kwargs.get("with_bounding_box", None)
+        fill_value = kwargs.get("fill_value", np.nan)
+
+        # If frames are of type ``str``, set the object to ``None``.
+        from_frame_obj = getattr(self, from_frame) if isinstance(from_frame, str) else from_frame
+        if isinstance(from_frame_obj, EmptyFrame):
+            from_frame_obj = None
+
+        to_frame_obj = getattr(self, to_frame) if isinstance(to_frame, str) else to_frame
+        if isinstance(to_frame_obj, EmptyFrame):
+            to_frame_obj = None
 
-        return self._evaluate_transform(
-            transform, from_step.step.frame, to_step.step.frame, *args, **kwargs
+        input_is_quantity, transform_uses_quantity = self._units_are_present(args, transform)
+        args = self._make_input_units_consistent(
+            transform,
+            *args,
+            frame=from_frame_obj,
+            input_is_quantity=input_is_quantity,
+            transform_uses_quantity=transform_uses_quantity,
         )
 
+        result = transform(*args, with_bounding_box=with_bounding_box, fill_value=fill_value, **kwargs)
+        if to_frame_obj is not None:
+            if to_frame_obj.naxes == 1:
+                result = (result,)
+            result = self._make_output_units_consistent(
+                transform,
+                *result,
+                frame=to_frame_obj,
+                input_is_quantity=input_is_quantity,
+                transform_uses_quantity=transform_uses_quantity,
+            )
+        if self.output_frame is not None and self.output_frame.naxes == 1:
+            return result[0]
+        return result
+
     @property
     def name(self) -> str:
         """Return the name for this WCS."""

From 7b6070bafe5b781a221c3a68952f0997e7c64e98 Mon Sep 17 00:00:00 2001
From: Nadia Dencheva <nadia.dencheva@gmail.com>
Date: Thu, 11 Dec 2025 16:50:23 -0500
Subject: [PATCH 02/15] remove old code; fix style

---
 gwcs/coordinate_frames/_base.py   |  3 +-
 gwcs/tests/test_api.py            |  5 +-
 gwcs/tests/test_api_consistent.py | 34 +++++++-----
 gwcs/tests/test_wcs.py            | 12 ++--
 gwcs/wcs/_wcs.py                  | 92 -------------------------------
 5 files changed, 32 insertions(+), 114 deletions(-)

diff --git a/gwcs/coordinate_frames/_base.py b/gwcs/coordinate_frames/_base.py
index d9873cc6..c31eed6b 100644
--- a/gwcs/coordinate_frames/_base.py
+++ b/gwcs/coordinate_frames/_base.py
@@ -144,8 +144,7 @@ def remove_units(
         if self.naxes == 1 and (np.isscalar(arrays) or isinstance(arrays, u.Quantity)):
             arrays = (arrays,)
 
-        result = tuple(
+        return tuple(
             array.to_value(unit) if isinstance(array, u.Quantity) else array
             for array, unit in zip(arrays, self.unit, strict=True)
         )
-        return result
diff --git a/gwcs/tests/test_api.py b/gwcs/tests/test_api.py
index 528462de..eaff2982 100644
--- a/gwcs/tests/test_api.py
+++ b/gwcs/tests/test_api.py
@@ -328,7 +328,8 @@ def test_high_level_wrapper(wcsobj, request):
     pix_out1 = hlvl.world_to_pixel(*wc1)
     np.testing.assert_allclose(pix_out1, pixel_input)
     with pytest.raises(TypeError) as e:
-        pix_out2 = wcsobj.invert(*wc1)
+        _ = wcsobj.invert(*wc1)
+    assert "High Level objects are not supported with the native" in str(e)
 
 
 def test_stokes_wrapper(gwcs_stokes_lookup):
@@ -591,7 +592,7 @@ def test_coordinate_frame_api():
     assert isinstance(pixel, float)
 
     with pytest.raises(TypeError):
-        pixel2 = wcs.invert(world)
+        _ = wcs.invert(world)
 
 
 def test_world_axis_object_components_units(gwcs_3d_identity_units):
diff --git a/gwcs/tests/test_api_consistent.py b/gwcs/tests/test_api_consistent.py
index f91e8a14..1f8677fb 100644
--- a/gwcs/tests/test_api_consistent.py
+++ b/gwcs/tests/test_api_consistent.py
@@ -11,8 +11,10 @@
 Rules:
 
 
-1. Neither transforms nor inputs support units -> the output is clearly numerical for all functions above
-2. Transforms support units but inputs do not -> return quantities assuming the units of the coordinate frame
+1. Neither transforms nor inputs support units -> the output is clearly numerical
+   for all functions above
+2. Transforms support units but inputs do not -> return quantities assuming the
+   units of the coordinate frame
   - This should work for the wcs methods (wcs(x,y) and wcs.invert
   - The methods using transforms should follow modeling rules and will require units
     on the input and raise an exception if not
@@ -41,7 +43,8 @@
 def is_numerical(args):
     if isinstance(args, numbers.Number):
         return True
-    return all([isinstance(arg, numbers.Number) or type(arg) == np.ndarray for arg in args])
+    #return all([isinstance(arg, numbers.Number) or type(arg) == np.ndarray for arg in args])
+    return all([isinstance(arg, numbers.Number) or arg is np.ndarray for arg in args])
 
 
 def is_quantity(args):
@@ -54,17 +57,20 @@ def wcsobj(request):
 
 wno_unit_1d = ["gwcs_1d_freq", "gwcs_1d_spectral",]
 
-wno_unit_nd = ["gwcs_2d_shift_scale", "gwcs_3d_spatial_wave", "gwcs_2d_spatial_shift", "gwcs_2d_spatial_reordered",
-    "gwcs_3d_spatial_wave", "gwcs_simple_imaging", "gwcs_3spectral_orders",
-    "gwcs_3d_galactic_spectral", "gwcs_spec_cel_time_4d", "gwcs_romanisim", ]
+wno_unit_nd = ["gwcs_2d_shift_scale", "gwcs_3d_spatial_wave", "gwcs_2d_spatial_shift",
+    "gwcs_2d_spatial_reordered", "gwcs_3d_spatial_wave", "gwcs_simple_imaging",
+    "gwcs_3spectral_orders", "gwcs_3d_galactic_spectral", "gwcs_spec_cel_time_4d",
+    "gwcs_romanisim", ]
 
 # "gwcs_7d_complex_mapping" errors in astropy - fix
-# "gwcs_2d_quantity_shift" errors when inputs are quantities. Need to confirm if Qs are HLO
+# "gwcs_2d_quantity_shift" errors when inputs are quantities.
+# Need to confirm if Qs are HLO
 
 w_unit_1d = ["gwcs_stokes_lookup", "gwcs_1d_freq_quantity"]
 
-w_unit_nd = ["gwcs_2d_shift_scale_quantity", "gwcs_3d_identity_units", "gwcs_3d_identity_units",
-    "gwcs_4d_identity_units", "gwcs_simple_imaging_units", "gwcs_with_pipeline_celestial", ]
+w_unit_nd = ["gwcs_2d_shift_scale_quantity", "gwcs_3d_identity_units",
+    "gwcs_3d_identity_units", "gwcs_4d_identity_units", "gwcs_simple_imaging_units",
+    "gwcs_with_pipeline_celestial", ]
 
 w_transform_test = ["gwcs_1d_freq_quantity", "gwcs_2d_quantity_shift"]
 
@@ -75,8 +81,8 @@ def wcsobj(request):
 
 
 @wcs_no_unit_1d
-def test_trnou_inpnou_1d(wcsobj):
-    """ Transforms do not support units. Inputs are numbers."""
+def test_no_units_1d(wcsobj):
+    """ Transforms do not support units."""
     assert not wcsobj.forward_transform.uses_quantity
 
     # the case of a scalar input
@@ -112,14 +118,14 @@ def test_no_units_nd(wcsobj):
     if np.isscalar(result):
         result = [result]
     inp_new = wcsobj.invert(*result)
-    _ = [assert_allclose(i, j) for i, j in zip(inp_new, inp)]
+    _ = [assert_allclose(i, j) for i, j in zip(inp_new, inp, strict=True)]
 
     # Inputs are quantities; return quantities (except for pixels?)
-    inpq = [coo * un for coo, un in zip(inp, wcsobj.input_frame.unit)]
+    inpq = [coo * un for coo, un in zip(inp, wcsobj.input_frame.unit, strict=True)]
     result = wcsobj(*inpq)
     assert is_quantity(result)
     inp_new = wcsobj.invert(*result)
-    _ = [assert_allclose(i, j) for i, j in zip(inp_new, inpq)]
+    _ = [assert_allclose(i, j) for i, j in zip(inp_new, inpq, strict=True)]
 
     # input is HLO - raise an Error
     sky = wcsobj.pixel_to_world(*inp)
diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py
index 377d031a..1792d70e 100644
--- a/gwcs/tests/test_wcs.py
+++ b/gwcs/tests/test_wcs.py
@@ -1728,10 +1728,12 @@ def test_quantities_in_pipeline_backward(gwcs_with_pipeline_celestial):
         20 * u.arcsec + 1 * u.deg,
         15 * u.deg + 2 * u.deg,
     ]
-    pixel = iwcs.invert(*input_world)
+    with pytest.raises(TypeError) as e:
+        pixel = iwcs.invert(*input_world)
+    assert "High Level objects are not supported with the native" in str(e)
 
-    assert all(isinstance(p, u.Quantity) for p in pixel)
-    assert u.allclose(pixel, [1, 1] * u.pix)
+    # assert all(isinstance(p, u.Quantity) for p in pixel)
+    # assert u.allclose(pixel, [1, 1] * u.pix)
 
     intermediate_world = iwcs.transform(
         "output",
@@ -1882,7 +1884,9 @@ def test_parameterless_transform():
     assert gwcs.invert(1, 1) == (1, 1)
     # Strictly speaking it's correct that this fails Because
     # for this setup the HLO are Quantities
-    assert gwcs.invert(1 * u.pix, 1 * u.pix) == (1, 1)
+    with pytest.raises(TypeError) as e:
+        _ = gwcs.invert(1 * u.pix, 1 * u.pix)
+    assert "High Level objects are not supported with the native" in str(e)
 
 
 def test_fitswcs_imaging(fits_wcs_imaging_simple):
diff --git a/gwcs/wcs/_wcs.py b/gwcs/wcs/_wcs.py
index 1c34bef2..9133559b 100644
--- a/gwcs/wcs/_wcs.py
+++ b/gwcs/wcs/_wcs.py
@@ -158,9 +158,6 @@ def __call__(
         if transform is None:
             msg = "Transform is not defined."
             raise NotImplementedError(msg)
-        # move this to an evaluate_funciton, called by forward, transform and invert
-        # input_is_quantity = any(isinstance(a, u.Quantity) for a in args)
-        # transform_uses_quantity = not (transform is None or not transform.uses_quantity)
 
         input_is_quantity, transform_uses_quantity = self._units_are_present(args, transform)
         args = self._make_input_units_consistent(
@@ -237,7 +234,6 @@ def _make_input_units_consistent(
             return self._add_units_input(args, frame)
         elif not transform_uses_quantity and input_is_quantity:
             return self._remove_units_input(args, frame)
-        #return args
 
     def _make_output_units_consistent(
         self,
@@ -268,94 +264,6 @@ def _make_output_units_consistent(
             result = self._add_units_input(args, frame)
         return result
 
-    def _get_transform(
-        self,
-        *args,
-        from_frame: CoordinateFrame | None = None,
-        to_frame: CoordinateFrame | None = None,
-        **kwargs,
-    ):
-        """
-        Executes the forward transform, but values only.
-        """
-        if from_frame is None and to_frame is None:
-            transform = self.forward_transform
-        else:
-            transform = self.get_transform(from_frame, to_frame)
-
-        if transform is None:
-            msg = "WCS.forward_transform is not implemented."
-            raise NotImplementedError(msg)
-        return transform
-
-    def _evaluate_transform(
-        self,
-        transform,
-        from_frame,
-        to_frame,
-        *args,
-        with_bounding_box: bool = True,
-        fill_value: float | np.number = np.nan,
-        **kwargs,
-    ):
-        """
-        Introduces or removes units from the arguments as need so that the transform
-        can be successfully evaluated.
-
-        Notes
-        -----
-        Much of the logic in this method is due to the unfortunate fact that the
-        `uses_quantity` property for models is not reliable for determining if one
-        must pass quantities or not. It instead tells you:
-            1. If it has any parameter that is a quantity
-            2. It defaults to true for parameterless models.
-
-        This is problematic because its entirely possible to construct a model with
-        a parameter that is a quantity but the model itself either doesn't require
-        them or in fact cannot use them. This is a very rare case but it could happen.
-        Currently, this case is not handled, but it is worth noting in case it comes up
-
-        The more problematic case is for parameterless models. `uses_quantity` assumes
-        that if there are no parameters, then the model is agnostic to quantity inputs.
-        This is an incorrect assumption, even with in `astropy.modeling`'s built in
-        models.  The `Tabular1D` model for example has no "parameters" but it can
-        require quantities if its "points" construction input is a quantity. This
-        is the main case for the try/except block in this method.
-
-        Properly dealing with this will require upstream work in `astropy.modeling`
-        which is outside the scope of what GWCS can control.
-
-        to_frame is included as we really ought to be stripping the result of units
-        but we currently are not. API refactor should include addressing this.
-        """
-
-        # Validate that the input type matches what the transform expects
-        input_is_quantity = any(isinstance(a, u.Quantity) for a in args)
-
-        def _transform(*args):
-            """Wrap the transform evaluation"""
-
-            return transform(
-                *args,
-                with_bounding_box=with_bounding_box,
-                fill_value=fill_value,
-                **kwargs,
-            )
-
-        try:
-            return _transform(*args)
-        except u.UnitsError:
-            # In this case we are handling parameterless models that require units
-            # to function correctly.
-            if (
-                not input_is_quantity
-                and transform.uses_quantity
-                and not transform.parameters.size
-            ):
-                return _transform(*self._add_units_input(args, from_frame))
-
-            raise
-
     def in_image(self, *args, **kwargs):
         """
         This method tests if one or more of the input world coordinates are

From 43d80b07c33ddd314297b02d8f0822a8a94b5324 Mon Sep 17 00:00:00 2001
From: Nadia Dencheva <nadia.dencheva@gmail.com>
Date: Fri, 19 Dec 2025 14:28:18 -0500
Subject: [PATCH 03/15] update tests

---
 gwcs/examples.py                  | 10 +++++++++-
 gwcs/tests/conftest.py            |  5 +++++
 gwcs/tests/test_api_consistent.py | 26 ++++++++++++++++++++------
 3 files changed, 34 insertions(+), 7 deletions(-)

diff --git a/gwcs/examples.py b/gwcs/examples.py
index cdcd1b94..ea22120a 100644
--- a/gwcs/examples.py
+++ b/gwcs/examples.py
@@ -410,7 +410,7 @@ def gwcs_spec_cel_time_4d():
     wcslin = models.Mapping((1, 0)) | (offx & offy) | aff
     tan = models.Pix2Sky_TAN(name="tangent_projection")
     n2c = models.RotateNative2Celestial(*crval, 180, name="sky_rotation")
-    cel_model = wcslin | tan | n2c
+    cel_model = wcslin | tan | n2c | models.Mapping((1, 0))
     icrs = cf.CelestialFrame(
         reference_frame=coord.ICRS(), name="sky", axes_order=(2, 1)
     )
@@ -734,3 +734,11 @@ def fits_wcs_imaging_simple(params):
         w.wcs.lonpole = 180
     w.wcs.set()
     return gw, w
+
+
+def gwcs_2d_spatial_shift_reverse():
+    """
+    A simple one step spatial WCS with forward from sky to detector.
+    """
+    pipe = [(ICRC_SKY_FRAME, MODEL_2D_SHIFT), (DETECTOR_2D_FRAME, None)]
+    return wcs.WCS(pipe)
diff --git a/gwcs/tests/conftest.py b/gwcs/tests/conftest.py
index c410bf8a..d2167ccd 100644
--- a/gwcs/tests/conftest.py
+++ b/gwcs/tests/conftest.py
@@ -174,3 +174,8 @@ def gwcs_romanisim():
 def fits_wcs_imaging_simple(request):
     params = request.param
     return examples.fits_wcs_imaging_simple(params)
+
+
+@pytest.fixture
+def gwcs_2d_spatial_shift_reverse():
+    return examples.gwcs_2d_spatial_shift_reverse()
diff --git a/gwcs/tests/test_api_consistent.py b/gwcs/tests/test_api_consistent.py
index 1f8677fb..bd16f5da 100644
--- a/gwcs/tests/test_api_consistent.py
+++ b/gwcs/tests/test_api_consistent.py
@@ -32,7 +32,7 @@
 from astropy.tests.helper import assert_quantity_allclose
 
 import pytest
-from .conftest import gwcs_with_pipeline_celestial
+from .conftest import gwcs_with_pipeline_celestial, gwcs_2d_spatial_shift_reverse
 
 x = 1
 y = 2
@@ -127,10 +127,12 @@ def test_no_units_nd(wcsobj):
     inp_new = wcsobj.invert(*result)
     _ = [assert_allclose(i, j) for i, j in zip(inp_new, inpq, strict=True)]
 
-    # input is HLO - raise an Error
     sky = wcsobj.pixel_to_world(*inp)
+    if not np.iterable(sky):
+        sky=(sky,)
     with pytest.raises(TypeError) as e:
-        wcsobj.invert(*sky)
+        inv_sky = wcsobj.invert(*sky)
+    assert "High Level objects are not supported with the native" in str(e)
 
 
 @wcs_with_unit_1d
@@ -172,10 +174,12 @@ def test_transform_with_units(wcsobj):
     assert all([type(res)==u.Quantity for res in result])
     assert_allclose([r.value for r in result], result_num)
 
-    # input is HLO; raise an error
     sky = wcsobj.pixel_to_world(*xxq)
-    with pytest.raises(TypeError) as e:
-        wcsobj.invert(*sky)
+    if not np.iterable(sky):
+        sky=(sky,)
+    with pytest.raises(TypeError)as e:
+        inv_sky = wcsobj.invert(*sky)
+    assert "High Level objects are not supported with the native" in str(e)
 
 
 @wcs_no_unit_1d
@@ -213,6 +217,7 @@ def test_remove_units(wcsobj):
 
 
 def test_transform_multistage_wcs(gwcs_with_pipeline_celestial):
+    """Tests that the input and output types match for intermediate frames/transforms."""
     wcsobj = gwcs_with_pipeline_celestial
     frames = wcsobj.available_frames
     result = wcsobj.transform(frames[0], frames[-1], 1*u.pix, 1*u.pix)
@@ -227,3 +232,12 @@ def test_transform_multistage_wcs(gwcs_with_pipeline_celestial):
     assert_quantity_allclose(interm_result, tr(1*u.pix, 1*u.pix))
     ninterm_result = wcsobj.transform(frames[0], frames[1], 1, 1)
     assert_allclose([r.value for r in interm_result], ninterm_result)
+
+
+def test_reverse_wcs_direction(gwcs_2d_spatial_shift_reverse):
+    """Test that input quantities are converted to the units of the input frame."""
+    wcsobj = gwcs_2d_spatial_shift_reverse
+    assert_quantity_allclose(
+        wcsobj(1*u.arcsec, 2*u.arcsec),
+        wcsobj(1*u.arcsec.to(u.deg)*u.deg, 2*u.arcsec.to(u.deg)*u.deg)
+        )

From 491bec319b5361f188a7fb12150fba16303e2c90 Mon Sep 17 00:00:00 2001
From: Nadia Dencheva <nadia.dencheva@gmail.com>
Date: Fri, 19 Dec 2025 14:29:38 -0500
Subject: [PATCH 04/15] add a change log

---
 CHANGES.rst | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/CHANGES.rst b/CHANGES.rst
index 4289645f..10a2e925 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -7,6 +7,8 @@
 
 - Adjust and fix tests. Replace ``logging`` with `warnings``. [#659]
 
+- Update the legacy API. [#660]
+
 0.26.1 (2025-11-19)
 -------------------
 - Fix an indexing bug in ``spectroscopy.SellmeierZemax`` where the output ``n`` for array-type wavelength

From a49b40dc93bc971598ac48014f014caf00db33dd Mon Sep 17 00:00:00 2001
From: Nadia Dencheva <nadia.dencheva@gmail.com>
Date: Fri, 2 Jan 2026 16:34:38 -0500
Subject: [PATCH 05/15] address comments

---
 gwcs/wcs/_wcs.py | 12 +++++++-----
 1 file changed, 7 insertions(+), 5 deletions(-)

diff --git a/gwcs/wcs/_wcs.py b/gwcs/wcs/_wcs.py
index 9133559b..73355e1a 100644
--- a/gwcs/wcs/_wcs.py
+++ b/gwcs/wcs/_wcs.py
@@ -348,7 +348,7 @@ def invert(
             transform = self.backward_transform
         except NotImplementedError:
             transform = None
-            # TODO raise error?
+
         if is_high_level(*args, low_level_wcs=self):
             # args = high_level_objects_to_values(*args, low_level_wcs=self)
             message = "High Level objects are not supported with the native API. \
@@ -1112,7 +1112,9 @@ def transform(
         from_frame: str | CoordinateFrame,
         to_frame: str | CoordinateFrame,
         *args,
-        **kwargs,
+        with_bounding_box: bool = True,
+        fill_value: float | np.number = np.nan,
+        **kwargs
     ):
         """
         Transform positions between two frames.
@@ -1128,6 +1130,9 @@ def transform(
         with_bounding_box : bool, optional
              If True(default) values in the result which correspond to any of
              the inputs being outside the bounding_box are set to ``fill_value``.
+        fill_value : float, optional
+            Output value for inputs outside the bounding_box
+            (default is np.nan).
         """
         # Pull the steps and their indices from the pipeline
         # -> this also turns the frame name strings into frame objects
@@ -1135,9 +1140,6 @@ def transform(
         to_step = self._get_step(to_frame)
         transform = self.get_transform(from_step.step.frame, to_step.step.frame)
 
-        with_bounding_box = kwargs.get("with_bounding_box", None)
-        fill_value = kwargs.get("fill_value", np.nan)
-
         # If frames are of type ``str``, set the object to ``None``.
         from_frame_obj = getattr(self, from_frame) if isinstance(from_frame, str) else from_frame
         if isinstance(from_frame_obj, EmptyFrame):

From 7f2e5b3b4d4089311eac080f873703c249669608 Mon Sep 17 00:00:00 2001
From: Nadia Dencheva <nadia.dencheva@gmail.com>
Date: Fri, 2 Jan 2026 16:36:49 -0500
Subject: [PATCH 06/15] run precommit

---
 gwcs/tests/test_api.py            |   1 +
 gwcs/tests/test_api_consistent.py | 106 ++++++++++++++++++------------
 gwcs/tests/test_wcs.py            |   9 ++-
 gwcs/wcs/_wcs.py                  |  61 ++++++++++-------
 4 files changed, 113 insertions(+), 64 deletions(-)

diff --git a/gwcs/tests/test_api.py b/gwcs/tests/test_api.py
index eaff2982..9dc4c4c8 100644
--- a/gwcs/tests/test_api.py
+++ b/gwcs/tests/test_api.py
@@ -402,6 +402,7 @@ def test_pixel_bounds(wcsobj):
     # Reset the bounding box or this will affect other tests
     wcsobj.bounding_box = None
 
+
 @wcs_objs
 def test_axis_correlation_matrix(wcsobj):
     assert_array_equal(wcsobj.axis_correlation_matrix, np.identity(2))
diff --git a/gwcs/tests/test_api_consistent.py b/gwcs/tests/test_api_consistent.py
index bd16f5da..72b510c4 100644
--- a/gwcs/tests/test_api_consistent.py
+++ b/gwcs/tests/test_api_consistent.py
@@ -24,15 +24,14 @@
 4. Transforms do not support units but inputs are quantities -> raise an error
 
 """
+
 import numbers
-import numpy as np
-from numpy.testing import assert_array_equal, assert_allclose
 
+import numpy as np
+import pytest
 from astropy import units as u
 from astropy.tests.helper import assert_quantity_allclose
-
-import pytest
-from .conftest import gwcs_with_pipeline_celestial, gwcs_2d_spatial_shift_reverse
+from numpy.testing import assert_allclose
 
 x = 1
 y = 2
@@ -43,7 +42,7 @@
 def is_numerical(args):
     if isinstance(args, numbers.Number):
         return True
-    #return all([isinstance(arg, numbers.Number) or type(arg) == np.ndarray for arg in args])
+    # return all([isinstance(arg, numbers.Number) or type(arg) == np.ndarray for arg in args])
     return all([isinstance(arg, numbers.Number) or arg is np.ndarray for arg in args])
 
 
@@ -55,12 +54,24 @@ def is_quantity(args):
 def wcsobj(request):
     return request.getfixturevalue(request.param)
 
-wno_unit_1d = ["gwcs_1d_freq", "gwcs_1d_spectral",]
 
-wno_unit_nd = ["gwcs_2d_shift_scale", "gwcs_3d_spatial_wave", "gwcs_2d_spatial_shift",
-    "gwcs_2d_spatial_reordered", "gwcs_3d_spatial_wave", "gwcs_simple_imaging",
-    "gwcs_3spectral_orders", "gwcs_3d_galactic_spectral", "gwcs_spec_cel_time_4d",
-    "gwcs_romanisim", ]
+wno_unit_1d = [
+    "gwcs_1d_freq",
+    "gwcs_1d_spectral",
+]
+
+wno_unit_nd = [
+    "gwcs_2d_shift_scale",
+    "gwcs_3d_spatial_wave",
+    "gwcs_2d_spatial_shift",
+    "gwcs_2d_spatial_reordered",
+    "gwcs_3d_spatial_wave",
+    "gwcs_simple_imaging",
+    "gwcs_3spectral_orders",
+    "gwcs_3d_galactic_spectral",
+    "gwcs_spec_cel_time_4d",
+    "gwcs_romanisim",
+]
 
 # "gwcs_7d_complex_mapping" errors in astropy - fix
 # "gwcs_2d_quantity_shift" errors when inputs are quantities.
@@ -68,9 +79,14 @@ def wcsobj(request):
 
 w_unit_1d = ["gwcs_stokes_lookup", "gwcs_1d_freq_quantity"]
 
-w_unit_nd = ["gwcs_2d_shift_scale_quantity", "gwcs_3d_identity_units",
-    "gwcs_3d_identity_units", "gwcs_4d_identity_units", "gwcs_simple_imaging_units",
-    "gwcs_with_pipeline_celestial", ]
+w_unit_nd = [
+    "gwcs_2d_shift_scale_quantity",
+    "gwcs_3d_identity_units",
+    "gwcs_3d_identity_units",
+    "gwcs_4d_identity_units",
+    "gwcs_simple_imaging_units",
+    "gwcs_with_pipeline_celestial",
+]
 
 w_transform_test = ["gwcs_1d_freq_quantity", "gwcs_2d_quantity_shift"]
 
@@ -82,7 +98,7 @@ def wcsobj(request):
 
 @wcs_no_unit_1d
 def test_no_units_1d(wcsobj):
-    """ Transforms do not support units."""
+    """Transforms do not support units."""
     assert not wcsobj.forward_transform.uses_quantity
 
     # the case of a scalar input
@@ -121,7 +137,7 @@ def test_no_units_nd(wcsobj):
     _ = [assert_allclose(i, j) for i, j in zip(inp_new, inp, strict=True)]
 
     # Inputs are quantities; return quantities (except for pixels?)
-    inpq = [coo * un for coo, un in zip(inp, wcsobj.input_frame.unit, strict=True)]
+    inpq = [coup * un for coup, un in zip(inp, wcsobj.input_frame.unit, strict=True)]
     result = wcsobj(*inpq)
     assert is_quantity(result)
     inp_new = wcsobj.invert(*result)
@@ -129,7 +145,7 @@ def test_no_units_nd(wcsobj):
 
     sky = wcsobj.pixel_to_world(*inp)
     if not np.iterable(sky):
-        sky=(sky,)
+        sky = (sky,)
     with pytest.raises(TypeError) as e:
         inv_sky = wcsobj.invert(*sky)
     assert "High Level objects are not supported with the native" in str(e)
@@ -137,7 +153,7 @@ def test_no_units_nd(wcsobj):
 
 @wcs_with_unit_1d
 def test_with_units_1d(wcsobj):
-    """ Transform do not support units."""
+    """Transform do not support units."""
     assert wcsobj.forward_transform.uses_quantity
 
     # the case of a scalar input
@@ -155,7 +171,7 @@ def test_with_units_1d(wcsobj):
 
 @wcs_with_unit_nd
 def test_transform_with_units(wcsobj):
-    """ Transforms support units."""
+    """Transforms support units."""
     assert wcsobj.forward_transform.uses_quantity
 
     n_inputs = wcsobj.input_frame.naxes
@@ -171,13 +187,13 @@ def test_transform_with_units(wcsobj):
     # input is quantities; return quantities
     xxq = [1 * u.pix] * n_inputs
     result = wcsobj(*xxq)
-    assert all([type(res)==u.Quantity for res in result])
+    assert all([type(res) == u.Quantity for res in result])
     assert_allclose([r.value for r in result], result_num)
 
     sky = wcsobj.pixel_to_world(*xxq)
     if not np.iterable(sky):
-        sky=(sky,)
-    with pytest.raises(TypeError)as e:
+        sky = (sky,)
+    with pytest.raises(TypeError) as e:
         inv_sky = wcsobj.invert(*sky)
     assert "High Level objects are not supported with the native" in str(e)
 
@@ -185,17 +201,22 @@ def test_transform_with_units(wcsobj):
 @wcs_no_unit_1d
 def test_add_units(wcsobj):
     if wcsobj.input_frame.naxes == 1:
-        assert wcsobj._add_units_input((1,), wcsobj.input_frame) == 1 * wcsobj.input_frame.unit[0]
+        assert (
+            wcsobj._add_units_input((1,), wcsobj.input_frame)
+            == 1 * wcsobj.input_frame.unit[0]
+        )
         assert_allclose(
             wcsobj._add_units_input(([1, 1],), wcsobj.input_frame),
-            ([1, 1] * wcsobj.input_frame.unit[0],))
+            ([1, 1] * wcsobj.input_frame.unit[0],),
+        )
     elif wcsobj.input_frame.naxes == 2:
         assert_quantity_allclose(
-            wcsobj._add_units_input((1, 1), wcsobj.input_frame),
-            (1*u.pix, 1*u.pix))
+            wcsobj._add_units_input((1, 1), wcsobj.input_frame), (1 * u.pix, 1 * u.pix)
+        )
         assert_quantity_allclose(
             wcsobj._add_units_input(([1, 1], [1, 1]), wcsobj.input_frame),
-            ([1, 1]*u.pix, [1, 1]*u.pix))
+            ([1, 1] * u.pix, [1, 1] * u.pix),
+        )
 
 
 @wcs_with_unit_1d
@@ -204,32 +225,35 @@ def test_remove_units(wcsobj):
         unit = wcsobj.input_frame.unit[0]
         assert wcsobj._remove_units_input(1 * unit, wcsobj.input_frame) == (1,)
         assert_allclose(
-            wcsobj._remove_units_input(([1, 1] * unit,), wcsobj.input_frame),
-            ([1, 1],))
+            wcsobj._remove_units_input(([1, 1] * unit,), wcsobj.input_frame), ([1, 1],)
+        )
     elif wcsobj.input_frame.naxes == 2:
         assert_quantity_allclose(
-            wcsobj._remove_units_input((1*u.pix, 1*u.pix), wcsobj.input_frame),
-            (1, 1)
-            )
+            wcsobj._remove_units_input((1 * u.pix, 1 * u.pix), wcsobj.input_frame),
+            (1, 1),
+        )
         assert_quantity_allclose(
-            wcsobj._remove_units_input(([1, 1]*u.pix, [1, 1]*u.pix), wcsobj.input_frame),
-            ([1, 1], [1, 1]))
+            wcsobj._remove_units_input(
+                ([1, 1] * u.pix, [1, 1] * u.pix), wcsobj.input_frame
+            ),
+            ([1, 1], [1, 1]),
+        )
 
 
 def test_transform_multistage_wcs(gwcs_with_pipeline_celestial):
     """Tests that the input and output types match for intermediate frames/transforms."""
     wcsobj = gwcs_with_pipeline_celestial
     frames = wcsobj.available_frames
-    result = wcsobj.transform(frames[0], frames[-1], 1*u.pix, 1*u.pix)
+    result = wcsobj.transform(frames[0], frames[-1], 1 * u.pix, 1 * u.pix)
     assert is_quantity(result)
     assert_allclose([r.value for r in result], wcsobj(1, 1))
-    final_result = wcsobj.transform(frames[0], frames[-1], 1*u.pix, 1*u.pix)
+    final_result = wcsobj.transform(frames[0], frames[-1], 1 * u.pix, 1 * u.pix)
     assert is_quantity(final_result)
     assert_allclose([r.value for r in final_result], wcsobj(1, 1))
-    interm_result = wcsobj.transform(frames[0], frames[1], 1*u.pix, 1*u.pix)
+    interm_result = wcsobj.transform(frames[0], frames[1], 1 * u.pix, 1 * u.pix)
     assert is_quantity(interm_result)
     tr = wcsobj.get_transform(frames[0], frames[1])
-    assert_quantity_allclose(interm_result, tr(1*u.pix, 1*u.pix))
+    assert_quantity_allclose(interm_result, tr(1 * u.pix, 1 * u.pix))
     ninterm_result = wcsobj.transform(frames[0], frames[1], 1, 1)
     assert_allclose([r.value for r in interm_result], ninterm_result)
 
@@ -238,6 +262,6 @@ def test_reverse_wcs_direction(gwcs_2d_spatial_shift_reverse):
     """Test that input quantities are converted to the units of the input frame."""
     wcsobj = gwcs_2d_spatial_shift_reverse
     assert_quantity_allclose(
-        wcsobj(1*u.arcsec, 2*u.arcsec),
-        wcsobj(1*u.arcsec.to(u.deg)*u.deg, 2*u.arcsec.to(u.deg)*u.deg)
-        )
+        wcsobj(1 * u.arcsec, 2 * u.arcsec),
+        wcsobj(1 * u.arcsec.to(u.deg) * u.deg, 2 * u.arcsec.to(u.deg) * u.deg),
+    )
diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py
index 1792d70e..e35101d3 100644
--- a/gwcs/tests/test_wcs.py
+++ b/gwcs/tests/test_wcs.py
@@ -516,7 +516,14 @@ def test_bounding_box_eval():
     pipeline = [
         (
             cf.CoordinateFrame(
-                naxes=3, axes_type=("PIXEL", "PIXEL", "PIXEL",), axes_order=(0, 1, 2), name="detector"
+                naxes=3,
+                axes_type=(
+                    "PIXEL",
+                    "PIXEL",
+                    "PIXEL",
+                ),
+                axes_order=(0, 1, 2),
+                name="detector",
             ),
             trans3,
         ),
diff --git a/gwcs/wcs/_wcs.py b/gwcs/wcs/_wcs.py
index 73355e1a..9745b8fa 100644
--- a/gwcs/wcs/_wcs.py
+++ b/gwcs/wcs/_wcs.py
@@ -159,7 +159,9 @@ def __call__(
             msg = "Transform is not defined."
             raise NotImplementedError(msg)
 
-        input_is_quantity, transform_uses_quantity = self._units_are_present(args, transform)
+        input_is_quantity, transform_uses_quantity = self._units_are_present(
+            args, transform
+        )
         args = self._make_input_units_consistent(
             transform,
             *args,
@@ -168,7 +170,9 @@ def __call__(
             transform_uses_quantity=transform_uses_quantity,
         )
 
-        result = transform(*args, with_bounding_box=with_bounding_box, fill_value=fill_value, **kwargs)
+        result = transform(
+            *args, with_bounding_box=with_bounding_box, fill_value=fill_value, **kwargs
+        )
         if self.output_frame is not None:
             if self.output_frame.naxes == 1:
                 result = (result,)
@@ -205,7 +209,6 @@ def _units_are_present(self, args, transform):
         transform_uses_quantity = not (transform is None or not transform.uses_quantity)
         return input_is_quantity, transform_uses_quantity
 
-
     def _make_input_units_consistent(
         self,
         transform,
@@ -222,17 +225,18 @@ def _make_input_units_consistent(
         # # Validate that the input type matches what the transform expects
         # input_is_quantity = any(isinstance(a, u.Quantity) for a in args)
         # transform_uses_quantity = not (transform is None or not transform.uses_quantity)
-        if not input_is_quantity and not transform_uses_quantity:
-            return args
-        elif input_is_quantity and transform_uses_quantity:
+        if (not input_is_quantity and not transform_uses_quantity) or (
+            input_is_quantity and transform_uses_quantity
+        ):
             return args
-        elif (
+        if (
             not input_is_quantity
-            and (transform_uses_quantity
-            or transform.parameters.size) # possibly remove this, check is in _units_are_present
+            and (
+                transform_uses_quantity or transform.parameters.size
+            )  # possibly remove this, check is in _units_are_present
         ):
             return self._add_units_input(args, frame)
-        elif not transform_uses_quantity and input_is_quantity:
+        if not transform_uses_quantity and input_is_quantity:
             return self._remove_units_input(args, frame)
 
     def _make_output_units_consistent(
@@ -251,13 +255,11 @@ def _make_output_units_consistent(
         if not input_is_quantity and not transform_uses_quantity:
             return args
 
-        elif input_is_quantity and transform_uses_quantity:
+        if input_is_quantity and transform_uses_quantity:
             # make sure the output is returned in the units of the output frame
             return self._add_units_input(args, frame)
-        elif (
-            not input_is_quantity
-            and (transform_uses_quantity
-            or transform.parameters.size)
+        if not input_is_quantity and (
+            transform_uses_quantity or transform.parameters.size
         ):
             result = self._remove_units_input(args, frame)
         elif not transform_uses_quantity and input_is_quantity:
@@ -358,7 +360,9 @@ def invert(
         if with_bounding_box and self.bounding_box is not None:
             args = self.outside_footprint(args)
 
-        input_is_quantity, transform_uses_quantity = self._units_are_present(args, transform)
+        input_is_quantity, transform_uses_quantity = self._units_are_present(
+            args, transform
+        )
 
         args = self._make_input_units_consistent(
             transform,
@@ -369,7 +373,12 @@ def invert(
         )
         if transform is not None:
             akwargs = {k: v for k, v in kwargs.items() if k not in _ITER_INV_KWARGS}
-            result = transform(*args, with_bounding_box=with_bounding_box, fill_value=fill_value, **akwargs)
+            result = transform(
+                *args,
+                with_bounding_box=with_bounding_box,
+                fill_value=fill_value,
+                **akwargs,
+            )
         else:
             # Always strip units for numerical inverse
             args = self._remove_units_input(args, self.output_frame)
@@ -1114,7 +1123,7 @@ def transform(
         *args,
         with_bounding_box: bool = True,
         fill_value: float | np.number = np.nan,
-        **kwargs
+        **kwargs,
     ):
         """
         Transform positions between two frames.
@@ -1141,15 +1150,21 @@ def transform(
         transform = self.get_transform(from_step.step.frame, to_step.step.frame)
 
         # If frames are of type ``str``, set the object to ``None``.
-        from_frame_obj = getattr(self, from_frame) if isinstance(from_frame, str) else from_frame
+        from_frame_obj = (
+            getattr(self, from_frame) if isinstance(from_frame, str) else from_frame
+        )
         if isinstance(from_frame_obj, EmptyFrame):
             from_frame_obj = None
 
-        to_frame_obj = getattr(self, to_frame) if isinstance(to_frame, str) else to_frame
+        to_frame_obj = (
+            getattr(self, to_frame) if isinstance(to_frame, str) else to_frame
+        )
         if isinstance(to_frame_obj, EmptyFrame):
             to_frame_obj = None
 
-        input_is_quantity, transform_uses_quantity = self._units_are_present(args, transform)
+        input_is_quantity, transform_uses_quantity = self._units_are_present(
+            args, transform
+        )
         args = self._make_input_units_consistent(
             transform,
             *args,
@@ -1158,7 +1173,9 @@ def transform(
             transform_uses_quantity=transform_uses_quantity,
         )
 
-        result = transform(*args, with_bounding_box=with_bounding_box, fill_value=fill_value, **kwargs)
+        result = transform(
+            *args, with_bounding_box=with_bounding_box, fill_value=fill_value, **kwargs
+        )
         if to_frame_obj is not None:
             if to_frame_obj.naxes == 1:
                 result = (result,)

From a0ba1361335ab34cd80cbd287b1fa3a05a25e331 Mon Sep 17 00:00:00 2001
From: Nadia Dencheva <nadia.dencheva@gmail.com>
Date: Fri, 2 Jan 2026 16:54:05 -0500
Subject: [PATCH 07/15] precommit run

---
 gwcs/tests/test_api_consistent.py | 12 +++++++-----
 gwcs/tests/test_wcs.py            |  3 ---
 gwcs/wcs/_wcs.py                  | 25 ++++++++++---------------
 3 files changed, 17 insertions(+), 23 deletions(-)

diff --git a/gwcs/tests/test_api_consistent.py b/gwcs/tests/test_api_consistent.py
index 72b510c4..437e67ee 100644
--- a/gwcs/tests/test_api_consistent.py
+++ b/gwcs/tests/test_api_consistent.py
@@ -42,7 +42,6 @@
 def is_numerical(args):
     if isinstance(args, numbers.Number):
         return True
-    # return all([isinstance(arg, numbers.Number) or type(arg) == np.ndarray for arg in args])
     return all([isinstance(arg, numbers.Number) or arg is np.ndarray for arg in args])
 
 
@@ -177,17 +176,17 @@ def test_transform_with_units(wcsobj):
     n_inputs = wcsobj.input_frame.naxes
     xx = [x] * n_inputs
 
-    # input is numerical; return numbers
+    # input is numerical, return numbers
     result_num = wcsobj(*xx)
     assert is_numerical(result_num)
 
     inp = wcsobj.invert(*result_num)
     assert is_numerical(inp)
 
-    # input is quantities; return quantities
+    # input is quantities, return quantities
     xxq = [1 * u.pix] * n_inputs
     result = wcsobj(*xxq)
-    assert all([type(res) == u.Quantity for res in result])
+    assert all([type(res) is u.Quantity for res in result])
     assert_allclose([r.value for r in result], result_num)
 
     sky = wcsobj.pixel_to_world(*xxq)
@@ -241,7 +240,10 @@ def test_remove_units(wcsobj):
 
 
 def test_transform_multistage_wcs(gwcs_with_pipeline_celestial):
-    """Tests that the input and output types match for intermediate frames/transforms."""
+    """
+    Tests that the input and output types match for
+    intermediate frames/transforms.
+    """
     wcsobj = gwcs_with_pipeline_celestial
     frames = wcsobj.available_frames
     result = wcsobj.transform(frames[0], frames[-1], 1 * u.pix, 1 * u.pix)
diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py
index e35101d3..4fc64a92 100644
--- a/gwcs/tests/test_wcs.py
+++ b/gwcs/tests/test_wcs.py
@@ -1739,9 +1739,6 @@ def test_quantities_in_pipeline_backward(gwcs_with_pipeline_celestial):
         pixel = iwcs.invert(*input_world)
     assert "High Level objects are not supported with the native" in str(e)
 
-    # assert all(isinstance(p, u.Quantity) for p in pixel)
-    # assert u.allclose(pixel, [1, 1] * u.pix)
-
     intermediate_world = iwcs.transform(
         "output",
         "celestial",
diff --git a/gwcs/wcs/_wcs.py b/gwcs/wcs/_wcs.py
index 9745b8fa..44549cdf 100644
--- a/gwcs/wcs/_wcs.py
+++ b/gwcs/wcs/_wcs.py
@@ -222,22 +222,18 @@ def _make_input_units_consistent(
         Adds or removes units from the arguments as needed so that the transform
         can be successfully evaluated.
         """
-        # # Validate that the input type matches what the transform expects
-        # input_is_quantity = any(isinstance(a, u.Quantity) for a in args)
-        # transform_uses_quantity = not (transform is None or not transform.uses_quantity)
+        # Validate that the input type matches what the transform expects
         if (not input_is_quantity and not transform_uses_quantity) or (
             input_is_quantity and transform_uses_quantity
         ):
             return args
-        if (
-            not input_is_quantity
-            and (
-                transform_uses_quantity or transform.parameters.size
-            )  # possibly remove this, check is in _units_are_present
+        if not input_is_quantity and (
+            transform_uses_quantity or transform.parameters.size
         ):
             return self._add_units_input(args, frame)
         if not transform_uses_quantity and input_is_quantity:
             return self._remove_units_input(args, frame)
+        return args
 
     def _make_output_units_consistent(
         self,
@@ -249,8 +245,8 @@ def _make_output_units_consistent(
         **kwargs,
     ):
         """
-        Adds or removes units from the arguments as needed so that the type of the output
-        matches the input.
+        Adds or removes units from the arguments as needed so that
+        the type of the output matches the input.
         """
         if not input_is_quantity and not transform_uses_quantity:
             return args
@@ -261,10 +257,10 @@ def _make_output_units_consistent(
         if not input_is_quantity and (
             transform_uses_quantity or transform.parameters.size
         ):
-            result = self._remove_units_input(args, frame)
-        elif not transform_uses_quantity and input_is_quantity:
-            result = self._add_units_input(args, frame)
-        return result
+            return self._remove_units_input(args, frame)
+        if not transform_uses_quantity and input_is_quantity:
+            return self._add_units_input(args, frame)
+        return args
 
     def in_image(self, *args, **kwargs):
         """
@@ -352,7 +348,6 @@ def invert(
             transform = None
 
         if is_high_level(*args, low_level_wcs=self):
-            # args = high_level_objects_to_values(*args, low_level_wcs=self)
             message = "High Level objects are not supported with the native API. \
                        Please use the `world_to_pixel` method."
             raise TypeError(message)

From eef7645b00cdfadf97de5d511428eeccf048eefb Mon Sep 17 00:00:00 2001
From: Nadia Dencheva <nadia.dencheva@gmail.com>
Date: Fri, 2 Jan 2026 17:01:09 -0500
Subject: [PATCH 08/15] clarify rules

---
 gwcs/tests/test_api_consistent.py | 24 +++++++++++-------------
 1 file changed, 11 insertions(+), 13 deletions(-)

diff --git a/gwcs/tests/test_api_consistent.py b/gwcs/tests/test_api_consistent.py
index 437e67ee..1d029ee9 100644
--- a/gwcs/tests/test_api_consistent.py
+++ b/gwcs/tests/test_api_consistent.py
@@ -2,26 +2,24 @@
 """
 Test the API is consistent with units and quantities and follows the rules below.
 
-WCS functions which considered for this work are part of the legacy API:
-wcs(x, y)
-wcs.invert(ra, dec)
-wcs.forward_transform(x,y), wcs.backward_transform() and wcs.get_transform(f1, f2)
-wcs.numerical_inverse(ra, dec) - does not support units
+WCS functions considered part of the legacy API:
+wcs()
+wcs.invert()
+wcs.transform()
 
 Rules:
 
 
 1. Neither transforms nor inputs support units -> the output is clearly numerical
    for all functions above
-2. Transforms support units but inputs do not -> return quantities assuming the
-   units of the coordinate frame
-  - This should work for the wcs methods (wcs(x,y) and wcs.invert
-  - The methods using transforms should follow modeling rules and will require units
-    on the input and raise an exception if not
+2. Transforms support units but inputs do not -> return numbers/arrays
+   Attach the units of the input coordinate frame to the inputs.
+   Evaluate the transforms and strip the output of units
 3. Both transforms and inputs support units -> return quantities
-  - Wcs methods return quantities
-  - Transforms work and return quantities
-4. Transforms do not support units but inputs are quantities -> raise an error
+4. Transforms do not support units but inputs are quantities -> return quantities
+   Strip the units from the inputs after converting to the units of the input frame.
+   Evaluate the transform and attach the units of the output frame.
+5. Inputs are High Level Objects - raise an error
 
 """
 

From b2792a6ed4d982df57dcfebf71d4123130ed058f Mon Sep 17 00:00:00 2001
From: William Jamieson <william@jamieson.phd>
Date: Mon, 5 Jan 2026 11:21:10 -0500
Subject: [PATCH 09/15] Fix the style checks

---
 gwcs/tests/test_api_consistent.py | 24 +++++++++++++-----------
 gwcs/tests/test_wcs.py            |  7 ++++---
 gwcs/wcs/_wcs.py                  |  2 +-
 3 files changed, 18 insertions(+), 15 deletions(-)

diff --git a/gwcs/tests/test_api_consistent.py b/gwcs/tests/test_api_consistent.py
index 1d029ee9..8d4279f3 100644
--- a/gwcs/tests/test_api_consistent.py
+++ b/gwcs/tests/test_api_consistent.py
@@ -38,13 +38,13 @@
 
 
 def is_numerical(args):
-    if isinstance(args, numbers.Number):
-        return True
-    return all([isinstance(arg, numbers.Number) or arg is np.ndarray for arg in args])
+    return isinstance(args, numbers.Number) or all(
+        isinstance(arg, numbers.Number) or arg is np.ndarray for arg in args
+    )
 
 
 def is_quantity(args):
-    return all([isinstance(arg, u.Quantity) for arg in args])
+    return all(isinstance(arg, u.Quantity) for arg in args)
 
 
 @pytest.fixture
@@ -143,9 +143,10 @@ def test_no_units_nd(wcsobj):
     sky = wcsobj.pixel_to_world(*inp)
     if not np.iterable(sky):
         sky = (sky,)
-    with pytest.raises(TypeError) as e:
-        inv_sky = wcsobj.invert(*sky)
-    assert "High Level objects are not supported with the native" in str(e)
+    with pytest.raises(
+        TypeError, match=r"High Level objects are not supported with the native"
+    ):
+        wcsobj.invert(*sky)
 
 
 @wcs_with_unit_1d
@@ -184,15 +185,16 @@ def test_transform_with_units(wcsobj):
     # input is quantities, return quantities
     xxq = [1 * u.pix] * n_inputs
     result = wcsobj(*xxq)
-    assert all([type(res) is u.Quantity for res in result])
+    assert all(type(res) is u.Quantity for res in result)
     assert_allclose([r.value for r in result], result_num)
 
     sky = wcsobj.pixel_to_world(*xxq)
     if not np.iterable(sky):
         sky = (sky,)
-    with pytest.raises(TypeError) as e:
-        inv_sky = wcsobj.invert(*sky)
-    assert "High Level objects are not supported with the native" in str(e)
+    with pytest.raises(
+        TypeError, match=r"High Level objects are not supported with the native"
+    ):
+        wcsobj.invert(*sky)
 
 
 @wcs_no_unit_1d
diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py
index 4fc64a92..fe2f3d21 100644
--- a/gwcs/tests/test_wcs.py
+++ b/gwcs/tests/test_wcs.py
@@ -1735,9 +1735,10 @@ def test_quantities_in_pipeline_backward(gwcs_with_pipeline_celestial):
         20 * u.arcsec + 1 * u.deg,
         15 * u.deg + 2 * u.deg,
     ]
-    with pytest.raises(TypeError) as e:
-        pixel = iwcs.invert(*input_world)
-    assert "High Level objects are not supported with the native" in str(e)
+    with pytest.raises(
+        TypeError, match=r"High Level objects are not supported with the native"
+    ):
+        iwcs.invert(*input_world)
 
     intermediate_world = iwcs.transform(
         "output",
diff --git a/gwcs/wcs/_wcs.py b/gwcs/wcs/_wcs.py
index 44549cdf..54808ee0 100644
--- a/gwcs/wcs/_wcs.py
+++ b/gwcs/wcs/_wcs.py
@@ -189,7 +189,7 @@ def __call__(
 
     def _units_are_present(self, args, transform):
         """
-        Determing if the inputs to a transform are quantities and the transform
+        Determining if the inputs to a transform are quantities and the transform
         supports units.
 
         Parameters

From bc76b633581dd0216194da2ff9d9566ded1bc470 Mon Sep 17 00:00:00 2001
From: Nadia Dencheva <nadia.dencheva@gmail.com>
Date: Wed, 7 Jan 2026 12:32:07 -0500
Subject: [PATCH 10/15] fix return values in transform method when frame is 1D

---
 gwcs/examples.py                  | 15 +++++++++++++++
 gwcs/tests/test_api_consistent.py |  5 +++++
 gwcs/wcs/_wcs.py                  |  4 ++--
 3 files changed, 22 insertions(+), 2 deletions(-)

diff --git a/gwcs/examples.py b/gwcs/examples.py
index ea22120a..420cf641 100644
--- a/gwcs/examples.py
+++ b/gwcs/examples.py
@@ -742,3 +742,18 @@ def gwcs_2d_spatial_shift_reverse():
     """
     pipe = [(ICRC_SKY_FRAME, MODEL_2D_SHIFT), (DETECTOR_2D_FRAME, None)]
     return wcs.WCS(pipe)
+
+
+def gwcs_multi_stage():
+    """
+    A 3-step pipeline where the intermediate step is 1D and the final is 2D.
+    """
+    tr1 = models.Shift(10)
+    tr2 = models.Mapping((0, 0)) | models.Scale(-2) & models.Scale(-1)
+    det=cf.CoordinateFrame(name='detector', naxes=1, unit=('pix',),
+        axes_type='SPATIAL', axes_order=(0,))
+    interm = cf.CoordinateFrame(name='interm', naxes=1, unit=('m',),
+        axes_type='SPATIAL', axes_order=(0,))
+    cel = cf.CelestialFrame(name='sky', axes_names=('ra', 'dec'))
+    w = wcs.WCS([(det, tr1), (interm, tr2),(cel, None)])
+    return w
diff --git a/gwcs/tests/test_api_consistent.py b/gwcs/tests/test_api_consistent.py
index 8d4279f3..0d27b798 100644
--- a/gwcs/tests/test_api_consistent.py
+++ b/gwcs/tests/test_api_consistent.py
@@ -267,3 +267,8 @@ def test_reverse_wcs_direction(gwcs_2d_spatial_shift_reverse):
         wcsobj(1 * u.arcsec, 2 * u.arcsec),
         wcsobj(1 * u.arcsec.to(u.deg) * u.deg, 2 * u.arcsec.to(u.deg) * u.deg),
     )
+
+
+def test_transfrom_intermediate_1d(gwcs_multi_stage):
+    wcsobj = gwcs_multi_stage
+    assert wcsobj.transform("detector", "interm", 1) == 11.0
diff --git a/gwcs/wcs/_wcs.py b/gwcs/wcs/_wcs.py
index 54808ee0..87fb8a31 100644
--- a/gwcs/wcs/_wcs.py
+++ b/gwcs/wcs/_wcs.py
@@ -189,7 +189,7 @@ def __call__(
 
     def _units_are_present(self, args, transform):
         """
-        Determining if the inputs to a transform are quantities and the transform
+        Determine if the inputs to a transform are quantities and the transform
         supports units.
 
         Parameters
@@ -1181,7 +1181,7 @@ def transform(
                 input_is_quantity=input_is_quantity,
                 transform_uses_quantity=transform_uses_quantity,
             )
-        if self.output_frame is not None and self.output_frame.naxes == 1:
+        if to_frame_obj is not None and to_frame_obj.naxes == 1:
             return result[0]
         return result
 

From 9c72435ec703e5686f76c69a9d18fb828895408a Mon Sep 17 00:00:00 2001
From: Nadia Dencheva <nadia.dencheva@gmail.com>
Date: Wed, 7 Jan 2026 12:38:06 -0500
Subject: [PATCH 11/15] add fixture

---
 gwcs/tests/conftest.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/gwcs/tests/conftest.py b/gwcs/tests/conftest.py
index d2167ccd..3908a502 100644
--- a/gwcs/tests/conftest.py
+++ b/gwcs/tests/conftest.py
@@ -179,3 +179,8 @@ def fits_wcs_imaging_simple(request):
 @pytest.fixture
 def gwcs_2d_spatial_shift_reverse():
     return examples.gwcs_2d_spatial_shift_reverse()
+
+
+@pytest.fixture
+def gwcs_multi_stage():
+    return examples.gwcs_multi_stage()

From 9cd89d8e9b16b097e094b508710e036abd03c51d Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Wed, 7 Jan 2026 17:44:19 +0000
Subject: [PATCH 12/15] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 gwcs/examples.py | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/gwcs/examples.py b/gwcs/examples.py
index 420cf641..a2269548 100644
--- a/gwcs/examples.py
+++ b/gwcs/examples.py
@@ -750,10 +750,12 @@ def gwcs_multi_stage():
     """
     tr1 = models.Shift(10)
     tr2 = models.Mapping((0, 0)) | models.Scale(-2) & models.Scale(-1)
-    det=cf.CoordinateFrame(name='detector', naxes=1, unit=('pix',),
-        axes_type='SPATIAL', axes_order=(0,))
-    interm = cf.CoordinateFrame(name='interm', naxes=1, unit=('m',),
-        axes_type='SPATIAL', axes_order=(0,))
-    cel = cf.CelestialFrame(name='sky', axes_names=('ra', 'dec'))
-    w = wcs.WCS([(det, tr1), (interm, tr2),(cel, None)])
+    det = cf.CoordinateFrame(
+        name="detector", naxes=1, unit=("pix",), axes_type="SPATIAL", axes_order=(0,)
+    )
+    interm = cf.CoordinateFrame(
+        name="interm", naxes=1, unit=("m",), axes_type="SPATIAL", axes_order=(0,)
+    )
+    cel = cf.CelestialFrame(name="sky", axes_names=("ra", "dec"))
+    w = wcs.WCS([(det, tr1), (interm, tr2), (cel, None)])
     return w

From 730132ab39541b3bea8bf99afed1f1975fbe3dcd Mon Sep 17 00:00:00 2001
From: Nadia Dencheva <nadia.dencheva@gmail.com>
Date: Wed, 7 Jan 2026 12:49:38 -0500
Subject: [PATCH 13/15] make pre-commit happy

---
 gwcs/examples.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/gwcs/examples.py b/gwcs/examples.py
index a2269548..00078e88 100644
--- a/gwcs/examples.py
+++ b/gwcs/examples.py
@@ -757,5 +757,4 @@ def gwcs_multi_stage():
         name="interm", naxes=1, unit=("m",), axes_type="SPATIAL", axes_order=(0,)
     )
     cel = cf.CelestialFrame(name="sky", axes_names=("ra", "dec"))
-    w = wcs.WCS([(det, tr1), (interm, tr2), (cel, None)])
-    return w
+    return wcs.WCS([(det, tr1), (interm, tr2), (cel, None)])

From 1f23184d0b30962c911891f9c0de6553a9dfffd4 Mon Sep 17 00:00:00 2001
From: Nadia Dencheva <nadia.dencheva@gmail.com>
Date: Wed, 7 Jan 2026 11:50:38 -0500
Subject: [PATCH 14/15] add an option to generate an hlo from an intermediate
 frame

---
 gwcs/wcs/_wcs.py | 13 ++++++++++++-
 1 file changed, 12 insertions(+), 1 deletion(-)

diff --git a/gwcs/wcs/_wcs.py b/gwcs/wcs/_wcs.py
index 87fb8a31..ce69e4e5 100644
--- a/gwcs/wcs/_wcs.py
+++ b/gwcs/wcs/_wcs.py
@@ -1118,6 +1118,7 @@ def transform(
         *args,
         with_bounding_box: bool = True,
         fill_value: float | np.number = np.nan,
+        hlo: bool = False,
         **kwargs,
     ):
         """
@@ -1132,11 +1133,13 @@ def transform(
         args : float or array-like
             Inputs in ``from_frame``, separate inputs for each dimension.
         with_bounding_box : bool, optional
-             If True(default) values in the result which correspond to any of
+             If True (default) values in the result which correspond to any of
              the inputs being outside the bounding_box are set to ``fill_value``.
         fill_value : float, optional
             Output value for inputs outside the bounding_box
             (default is np.nan).
+        hlo : bool, optional
+            If True, return a high level object.
         """
         # Pull the steps and their indices from the pipeline
         # -> this also turns the frame name strings into frame objects
@@ -1181,6 +1184,14 @@ def transform(
                 input_is_quantity=input_is_quantity,
                 transform_uses_quantity=transform_uses_quantity,
             )
+
+        if hlo and to_frame_obj is not None:
+            result = values_to_high_level_objects(
+                *result,
+                low_level_wcs=self,
+                object_classes=to_frame_obj.world_axis_object_classes,
+                object_components=to_frame_obj.world_axis_object_components)
+
         if to_frame_obj is not None and to_frame_obj.naxes == 1:
             return result[0]
         return result

From 2264cf32f307692b23a5f27e2368ef2f22d7ff4a Mon Sep 17 00:00:00 2001
From: Nadia Dencheva <nadia.dencheva@gmail.com>
Date: Wed, 7 Jan 2026 12:57:52 -0500
Subject: [PATCH 15/15] add a test

---
 gwcs/tests/test_api_consistent.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/gwcs/tests/test_api_consistent.py b/gwcs/tests/test_api_consistent.py
index 0d27b798..37154b61 100644
--- a/gwcs/tests/test_api_consistent.py
+++ b/gwcs/tests/test_api_consistent.py
@@ -272,3 +272,6 @@ def test_reverse_wcs_direction(gwcs_2d_spatial_shift_reverse):
 def test_transfrom_intermediate_1d(gwcs_multi_stage):
     wcsobj = gwcs_multi_stage
     assert wcsobj.transform("detector", "interm", 1) == 11.0
+
+    assert wcsobj.transform("detector", "interm", 1, hlo=False) == 11.0
+    wcsobj.transform("detector", "interm", 1, hlo=True) == 11.0 * u.m