From 8407bd07424fe1a7e535db383136da1df869f01c Mon Sep 17 00:00:00 2001 From: William Jamieson Date: Mon, 2 Mar 2026 15:30:03 -0500 Subject: [PATCH 1/7] Make evaluate and invert use the same code path for unit handling --- docs/gwcs/points_to_wcs.rst | 4 +- gwcs/coordinate_frames/_base.py | 58 ++++++- gwcs/tests/test_api.py | 15 +- gwcs/tests/test_api_consistent.py | 10 +- gwcs/tests/test_coordinate_systems.py | 3 + gwcs/tests/test_wcs.py | 15 +- gwcs/wcs/_wcs.py | 223 +++++++++++++++++++------- 7 files changed, 240 insertions(+), 88 deletions(-) diff --git a/docs/gwcs/points_to_wcs.rst b/docs/gwcs/points_to_wcs.rst index 5c3060ab..3004489c 100644 --- a/docs/gwcs/points_to_wcs.rst +++ b/docs/gwcs/points_to_wcs.rst @@ -73,8 +73,10 @@ easily work in both pixel and sky space, and transform between frames. The GWCS object, which by default when called executes for forward transformation, can be used to convert coordinates from pixel to world. +.. doctest-requires:: numpy>=2.0 + >>> gwcs_obj(36.235,642.215) # doctest: +FLOAT_CMP - (246.72158004206716, 43.46075091731673) + (np.float64(246.72158004206716), np.float64(43.46075091731673)) Or using the common WCS API diff --git a/gwcs/coordinate_frames/_base.py b/gwcs/coordinate_frames/_base.py index 284ff68f..5193883c 100644 --- a/gwcs/coordinate_frames/_base.py +++ b/gwcs/coordinate_frames/_base.py @@ -282,7 +282,57 @@ def remove_units( for array in self.add_units(arrays) ) - def to_high_level_coordinates(self, *values): + def is_high_level(self, *args) -> bool: + """ + Return `True` if the input coordinates are already high level objects + described by this frame. + + This is used by the low level WCS API in Astropy to determine whether + to call ``to_high_level_coordinates`` or not. + """ + + if (world_axis_object_classes := self.world_axis_object_classes) is None or len( + args + ) != len(world_axis_object_classes): + return False + + type_match = [] + for arg, world_axis_object_class in zip( + args, world_axis_object_classes.values(), strict=True + ): + if isinstance(class_object := world_axis_object_class.class_object, str): + type_match.append( + type(arg).__name__ == class_object + and class_object != u.Quantity.__name__ + ) + else: + type_match.append( + isinstance(arg, class_object) and class_object is not u.Quantity + ) + + if all(type_match): + return True + + if any(type_match): + types = [ + ( + type(arg).__name__, + c.class_object + if isinstance(c.class_object, str) + else c.class_object.__name__, + ) + for arg, c in zip(args, world_axis_object_classes.values(), strict=True) + ] + msg = ( + "Invalid types were passed, got " + f"({', '.join(t[0] for t in types)}), but expected " + f"({', '.join(t[1] for t in types)})." + ) + raise TypeError(msg) + + return False + + def to_high_level_coordinates(self, *values, correct_1d=True): """ Convert "values" to high level coordinate objects described by this frame. @@ -309,11 +359,11 @@ def to_high_level_coordinates(self, *values): raise TypeError(msg) high_level = values_to_high_level_objects(*values, low_level_wcs=self) - if len(high_level) == 1: + if correct_1d and len(high_level) == 1: high_level = high_level[0] return high_level - def from_high_level_coordinates(self, *high_level_coords): + def from_high_level_coordinates(self, *high_level_coords, correct_1d=True): """ Convert high level coordinate objects to "values" as described by this frame. @@ -332,7 +382,7 @@ def from_high_level_coordinates(self, *high_level_coords): ``naxis`` number of coordinates as scalars or arrays. """ values = high_level_objects_to_values(*high_level_coords, low_level_wcs=self) - if len(values) == 1: + if correct_1d and self.naxes == 1: values = values[0] return values diff --git a/gwcs/tests/test_api.py b/gwcs/tests/test_api.py index e46f47a7..2121e6a1 100644 --- a/gwcs/tests/test_api.py +++ b/gwcs/tests/test_api.py @@ -326,10 +326,15 @@ def test_high_level_wrapper(wcsobj, request): wc1 = (wc1,) pix_out1 = hlvl.world_to_pixel(*wc1) + pix_out2 = wcsobj.invert(*wc1) + + pix_out2 = wcsobj.input_frame.remove_units(pix_out2) + + if not isinstance(pix_out2, list | tuple): + pix_out2 = (pix_out2,) + np.testing.assert_allclose(pix_out1, pixel_input) - with pytest.raises(TypeError) as e: - _ = wcsobj.invert(*wc1) - assert "High Level objects are not supported with the native" in str(e) + np.testing.assert_allclose(pix_out2, pixel_input) def test_stokes_wrapper(gwcs_stokes_lookup): @@ -592,8 +597,8 @@ def test_coordinate_frame_api(): pixel = wcs.world_to_pixel(world) assert isinstance(pixel, float) - with pytest.raises(TypeError): - _ = wcs.invert(world) + pixel2 = wcs.invert(world) + assert u.allclose(pixel2, 0 * u.pix) 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 35cbc288..f6c8cbfe 100644 --- a/gwcs/tests/test_api_consistent.py +++ b/gwcs/tests/test_api_consistent.py @@ -143,10 +143,7 @@ def test_no_units_nd(wcsobj): sky = wcsobj.pixel_to_world(*inp) if not np.iterable(sky): sky = (sky,) - with pytest.raises( - TypeError, match=r"High Level objects are not supported with the native" - ): - wcsobj.invert(*sky) + assert u.allclose(inpq, wcsobj.invert(*sky)) @wcs_with_unit_1d @@ -191,10 +188,7 @@ def test_transform_with_units(wcsobj): sky = wcsobj.pixel_to_world(*xxq) if not np.iterable(sky): sky = (sky,) - with pytest.raises( - TypeError, match=r"High Level objects are not supported with the native" - ): - wcsobj.invert(*sky) + assert u.allclose(xxq, wcsobj.invert(*sky)) @wcs_no_unit_1d diff --git a/gwcs/tests/test_coordinate_systems.py b/gwcs/tests/test_coordinate_systems.py index 3e603fab..65143515 100644 --- a/gwcs/tests/test_coordinate_systems.py +++ b/gwcs/tests/test_coordinate_systems.py @@ -646,6 +646,9 @@ def add_units(self, arrays): def remove_units(self, arrays): return arrays + def is_high_level(self, *inputs): + return False + def to_high_level_coordinates(self, *inputs): return inputs diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py index 21e37152..5ad7ae63 100644 --- a/gwcs/tests/test_wcs.py +++ b/gwcs/tests/test_wcs.py @@ -1788,10 +1788,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, match=r"High Level objects are not supported with the native" - ): - iwcs.invert(*input_world) + pixel = iwcs.invert(*input_world) + + assert all(isinstance(p, u.Quantity) for p in pixel) + assert u.allclose(pixel, [1, 1] * u.pix) intermediate_world = iwcs.transform( "output", @@ -1939,12 +1939,9 @@ def test_parameterless_transform(): assert gwcs(1, 1) == (1, 1) assert gwcs(1 * u.pix, 1 * u.pix) == (1 * u.pix, 1 * u.pix) + # No units introduced by the inverse transform assert gwcs.invert(1, 1) == (1, 1) - # Strictly speaking it's correct that this fails Because - # for this setup the HLO are Quantities - 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) + assert gwcs.invert(1 * u.pix, 1 * u.pix) == (1 * u.pix, 1 * u.pix) def test_fitswcs_imaging(fits_wcs_imaging_simple): diff --git a/gwcs/wcs/_wcs.py b/gwcs/wcs/_wcs.py index d6006edf..aa5b60fb 100644 --- a/gwcs/wcs/_wcs.py +++ b/gwcs/wcs/_wcs.py @@ -92,6 +92,155 @@ def __init__(self, axis, frame, world_axis_order, cunit, ctype, input_axes): self.input_axes = input_axes +class _UnitHandler: + """ + Class to ensure that the input-output type consistency for the GWCS native API + + Parameters + ---------- + args : + The input arguments to be passed to the transform. + _is_high_level : bool + Whether or not the inputs to the WCS API are high level objects. + _add_units : bool + If true, we add (or convert) units of the outputs of the transform to + match the input type. + If false, we remove (and convert if necessary) units of the outputs of + the transform to match the input type. + """ + + args: tuple[LowLevelInput, ...] + _is_high_level: bool + _add_units: bool + + def __init__( + self, + inputs: tuple[LowLevelInput, ...], + transform: Model, + frame: CoordinateFrameProtocol, + ): + # Handle the case where a HLO is passed into the Native API + self._is_high_level = frame.is_high_level(*inputs) + if self._is_high_level: + inputs = frame.from_high_level_coordinates(*inputs, correct_1d=False) + + # Determine if the inputs are quantities + input_is_quantity = any(isinstance(a, u.Quantity) for a in inputs) + + # Determine if the transform uses quantities + # This complexity is due to the fact that Tabular models have a bug + # in astropy + if isinstance(transform, (Tabular1D, Tabular2D)): + transform_uses_quantity = ( + isinstance(transform.lookup_table, u.Quantity) + or transform.input_units_equivalencies is not None + ) + + elif transform is not None and len(transform.parameters) == 0: + transform_uses_quantity = False + + else: + transform_uses_quantity = not ( + transform is None or not transform.uses_quantity + ) + + # Using input_is_quantity and transform_uses_quantity, we determine if we + # need to add or remove units from the inputs to make them compatible + # with the transform and we determine if we need to add or remove units + # from the outputs to make them consistent with the input's type. + + # Case 1: No units are involved + # input -> no change + # (we assume the user has passed numerical values + # in the correct units) + # output -> remove units + # (we convert the output to the correct units if necessary and + # then remove them giving numerical values as output) + if not input_is_quantity and not transform_uses_quantity: + self.args = inputs + self._add_units = False + + # Case 2: Units are involved and the transform supports them + # input -> add units + # (we add -- meaning convert -- units to the input values) + # output -> add units + # (we convert the output to the correct units if necessary and + elif input_is_quantity and transform_uses_quantity: + self.args = frame.add_units(inputs) + self._add_units = True + + # Case 3: no input units, but transform needs units + # input -> add units + # (we assume that the user has passed correct numerical values + # and formally add units to them for the transform to work + # correctly) + # output -> remove units + # (we convert the output to the correct units if necessary and + # then remove them giving numerical values as output) + elif not input_is_quantity and transform_uses_quantity: + self.args = frame.add_units(inputs) + self._add_units = False + + # Case 4: input has units, but transform does not use them + # input -> remove units + # (we convert the input to the correct units if necessary and + # then remove them giving numerical values as input) + # output -> add units + # (we convert the output to the correct units if necessary and + # then add them giving numerical values as output) + # This is the only other case: + # input_is_quantity and not transform_uses_quantity + else: + self.args = frame.remove_units(inputs) + self._add_units = True + + def _handle_output_type( + self, *outputs: LowLevelInput, frame: CoordinateFrameProtocol + ) -> tuple[LowLevelInput, ...]: + # Return a high level object if the input was a high level object + if self._is_high_level: + return frame.to_high_level_coordinates(*outputs, correct_1d=False) # type: ignore[no-any-return] + + # If the input had units return units + if self._add_units: + return frame.add_units(outputs) + + # If the input had no units, return unitless values + return frame.remove_units(outputs) + + def handle_output( + self, + outputs: tuple[LowLevelInput, ...] | LowLevelInput, + frame: CoordinateFrameProtocol, + ) -> tuple[LowLevelInput, ...] | LowLevelInput: + """ + Handle the output of a transform by adding or removing units as needed + and converting to high level objects if the input was high level. + + Parameters + ---------- + outputs : tuple of low level inputs + The outputs of a transform to be handled. + frame : CoordinateFrameProtocol + The frame to use for adding/removing units and converting to high level + objects. + + Returns + ------- + tuple of low level inputs or high level objects + The handled outputs. + """ + if frame.naxes == 1: + outputs = (outputs,) + + outputs = self._handle_output_type(*outputs, frame=frame) + + if frame.naxes == 1: + return outputs[0] + + return outputs + + class WCS(Pipeline, WCSAPIMixin): """ Basic WCS class. @@ -160,34 +309,16 @@ def evaluate( # after each call to it. transform = self.forward_transform - input_is_quantity, transform_uses_quantity = self._units_are_present( - args, transform - ) - args = self._make_input_units_consistent( - transform, - *args, - frame=self.input_frame, - input_is_quantity=input_is_quantity, - transform_uses_quantity=transform_uses_quantity, - ) + unit_handler = _UnitHandler(args, transform, self.input_frame) result = transform( - *args, with_bounding_box=with_bounding_box, fill_value=fill_value, **kwargs - ) - 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, + *unit_handler.args, + with_bounding_box=with_bounding_box, + fill_value=fill_value, + **kwargs, ) - if self.output_frame.naxes == 1: - return result[0] - return result + return unit_handler.handle_output(result, frame=self.output_frame) def __call__( self, @@ -266,9 +397,7 @@ def _make_input_units_consistent( input_is_quantity and transform_uses_quantity ): return args - if not input_is_quantity and ( - transform_uses_quantity or transform.parameters.size - ): + if not input_is_quantity and transform_uses_quantity: return frame.add_units(args) if not transform_uses_quantity and input_is_quantity: return frame.remove_units(args) @@ -394,40 +523,23 @@ def invert( except NotImplementedError: transform = None - if is_high_level(*args, low_level_wcs=self): - msg = ( - "High Level objects are not supported with the native API. " - "Please use the `world_to_pixel` method." - ) - raise TypeError(msg) + unit_handler = _UnitHandler(args, transform, self.output_frame) + inputs = unit_handler.args if with_bounding_box and self.bounding_box is not None: - args = self.outside_footprint(args) + inputs = self.outside_footprint(inputs) - 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: akwargs = {k: v for k, v in kwargs.items() if k not in _ITER_INV_KWARGS} result = transform( - *args, + *inputs, with_bounding_box=with_bounding_box, fill_value=fill_value, **akwargs, ) else: - # Always strip units for numerical inverse - args = self.output_frame.remove_units(args) result = self._numerical_inverse( - *args, + *inputs, with_bounding_box=with_bounding_box, fill_value=fill_value, **kwargs, @@ -436,18 +548,7 @@ def invert( 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 + return unit_handler.handle_output(result, frame=self.input_frame) def outside_footprint(self, world_arrays): world_arrays = [copy(array) for array in world_arrays] From 83de925b64a1554f219f8c9efbc1757355df7531 Mon Sep 17 00:00:00 2001 From: William Jamieson Date: Wed, 17 Jun 2026 16:32:42 -0400 Subject: [PATCH 2/7] Refactor pipeline and wcs.transform to use same internal code path as evaluate and invert --- gwcs/wcs/__init__.py | 3 + gwcs/wcs/_pipeline.py | 186 ++++++++++++++++++++++++++---------------- gwcs/wcs/_wcs.py | 44 +++------- 3 files changed, 130 insertions(+), 103 deletions(-) diff --git a/gwcs/wcs/__init__.py b/gwcs/wcs/__init__.py index 60045198..9bf22ced 100644 --- a/gwcs/wcs/__init__.py +++ b/gwcs/wcs/__init__.py @@ -1,10 +1,13 @@ from ._exception import GwcsBoundingBoxWarning, NoConvergence +from ._pipeline import DirectionalPipeline, Pipeline from ._step import Step from ._wcs import WCS __all__ = [ "WCS", + "DirectionalPipeline", "GwcsBoundingBoxWarning", "NoConvergence", + "Pipeline", "Step", ] diff --git a/gwcs/wcs/_pipeline.py b/gwcs/wcs/_pipeline.py index 3020dfb5..86ec8991 100644 --- a/gwcs/wcs/_pipeline.py +++ b/gwcs/wcs/_pipeline.py @@ -1,7 +1,8 @@ from __future__ import annotations import warnings -from typing import TypeAlias, Union, overload +from collections.abc import Sequence +from typing import Generic, NamedTuple, Self, TypeAlias, TypeVar, Union, overload from astropy.modeling import Model from astropy.modeling.bounding_box import CompoundBoundingBox, ModelBoundingBox @@ -13,10 +14,21 @@ from ._exception import GwcsBoundingBoxWarning, GwcsFrameExistsError from ._step import IndexedStep, Mdl, Step, StepTuple -__all__ = ["ForwardTransform", "Pipeline"] +__all__ = ["DirectionalPipeline", "ForwardTransform", "Pipeline"] # Type aliases due to the use of the `|` for type hints not working with Model -ForwardTransform: TypeAlias = Union[Model, list[Step | StepTuple]] # noqa: UP007 +ForwardTransform: TypeAlias = Union[Model, Sequence[Step | StepTuple]] # noqa: UP007 + +_T = TypeVar("_T") + + +class DirectionalPipeline(NamedTuple, Generic[_T]): + """ + A named tuple to hold a pipeline and the direction it is going in. + """ + + wcs: _T + forward: bool class Pipeline: @@ -40,7 +52,7 @@ def __init__( @overload def __init__( self, - forward_transform: list[Step | StepTuple], + forward_transform: Sequence[Step | StepTuple], *, input_frame: None = None, output_frame: None = None, @@ -103,7 +115,9 @@ def _initialize_pipeline( Step(input_frame, forward_transform.copy()), Step(output_frame, None), ] - case list(): + case Sequence() if not isinstance( + forward_transform, str | bytes | bytearray + ): if input_frame is not None and ( ( isinstance(forward_transform[0], Step) @@ -135,6 +149,9 @@ def _initialize_pipeline( "forward_transform pipeline." ) raise CoordinateFrameError(msg) + + # Normalize supported sequence inputs for downstream list-based APIs. + forward_transform = list(forward_transform) case _: msg = ( "Expected forward_transform to be a None, model, or a " @@ -319,6 +336,90 @@ def _get_step(self, frame: str | CoordinateFrameProtocol) -> IndexedStep: return IndexedStep(index, self._pipeline[index]) + @property + def forward_transform(self) -> Model: + """ + Return the forward transform of the pipeline. + """ + transform = combine_transforms([step.transform for step in self._pipeline[:-1]]) + + if self.bounding_box is not None: + # Currently compound models do not attempt to combine individual model + # bounding boxes. Get the forward transform and assign the bounding_box + # to it before evaluating it. The order Model.bounding_box is reversed. + transform.bounding_box = self.bounding_box + + return transform + + @property + def backward_transform(self) -> Model: + """ + Return the total backward transform if available - from output to input + coordinate system. + + Raises + ------ + NotImplementedError : + An analytical inverse does not exist. + + """ + try: + backward = self.forward_transform.inverse + except NotImplementedError: + # Try to get the backward transform by combining the inverses of the + # individual steps. This will work even if the forward transform does + # not have an analytical inverse, as long as each step does. + backward = combine_transforms( + [step.inverse for step in self._pipeline[:-1][::-1]] + ) + + try: + _ = backward.inverse + except NotImplementedError: # means "hasattr" won't work + backward.inverse = self.forward_transform + return backward + + def pipeline_between( + self, + from_frame: str | CoordinateFrameProtocol, + to_frame: str | CoordinateFrameProtocol, + ) -> DirectionalPipeline[Self] | None: + """ + Return a pipeline between the two given frames. + + Parameters + ---------- + from_frame : str or `~gwcs.coordinate_frames.CoordinateFrame` + Initial coordinate frame name of object. + to_frame : str or `~gwcs.coordinate_frames.CoordinateFrame` + End coordinate frame name or object. + + Returns + ------- + pipeline : `~gwcs.wcs.Pipeline` + A NamedTuple (pipeline, forward) where pipeline is a Pipeline object + the pipeline between the two frames and forward is a boolean indicating + if the pipeline is forward (from_frame to to_frame) or + backward (to_frame to from_frame). + """ + from_index = self._frame_index(from_frame) + to_index = self._frame_index(to_frame) + + # Moving backwards over the pipeline + if to_index < from_index: + pipeline = self._pipeline[to_index:from_index] + pipeline.append(Step(self._pipeline[from_index].frame, None)) + + return DirectionalPipeline(wcs=type(self)(pipeline), forward=False) + + if from_index < to_index: + pipeline = self._pipeline[from_index:to_index] + pipeline.append(Step(self._pipeline[to_index].frame, None)) + + return DirectionalPipeline(wcs=type(self)(pipeline), forward=True) + + return None # from and to are the same frame, so no pipeline needed + def get_transform( self, from_frame: str | CoordinateFrameProtocol, @@ -339,26 +440,18 @@ def get_transform( transform : `~astropy.modeling.Model` Transform between two frames. """ - from_index = self._frame_index(from_frame) - to_index = self._frame_index(to_frame) + direction = self.pipeline_between(from_frame, to_frame) - # Moving backwards over the pipeline - if to_index < from_index: - transforms = [ - step.inverse for step in self._pipeline[to_index:from_index][::-1] - ] - - # from and to are the same - elif to_index == from_index: + # from and to are the same frame, so no transform needed + if direction is None: return None - # Moving forwards over the pipeline - else: - transforms = [ - step.transform for step in self._pipeline[from_index:to_index] - ] + # If the pipeline is forward, return the forward transform + if direction.forward: + return direction.wcs.forward_transform - return combine_transforms(transforms) + # Otherwise it is backward, so return the backward transform + return direction.wcs.backward_transform def set_transform( self, @@ -505,8 +598,7 @@ def bounding_box(self) -> ModelBoundingBox | CompoundBoundingBox | None: """ # Pull the first transform of the pipeline which is what controls the # bounding_box - frames = self.available_frames - transform = self.get_transform(frames[0], frames[1]) + transform = self._pipeline[0].transform if transform is None: return None @@ -556,8 +648,7 @@ def bounding_box( value : tuple or None Tuple of tuples with ("low", high") values for the range. """ - frames = self.available_frames - transform = self.get_transform(frames[0], frames[1]) + transform = self._pipeline[0].transform if transform is None: msg = ( @@ -577,8 +668,6 @@ def bounding_box( transform.bounding_box = bbox - self.set_transform(frames[0], frames[1], transform) - def attach_compound_bounding_box( self, cbbox: dict[tuple[str, ...], tuple], selector_args: tuple[str, ...] ) -> None: @@ -594,47 +683,6 @@ def attach_compound_bounding_box( selector_args: Argument names to the model that are used to select the bounding box """ - frames = self.available_frames - transform_0 = self.get_transform(frames[0], frames[1]) - self.bounding_box = CompoundBoundingBox.validate( - transform_0, cbbox, selector_args=selector_args, order="F" + self.pipeline[0].transform, cbbox, selector_args=selector_args, order="F" ) - - @property - def forward_transform(self) -> Model: - """ - Return the forward transform of the pipeline. - """ - transform = combine_transforms([step.transform for step in self._pipeline[:-1]]) - - if self.bounding_box is not None: - # Currently compound models do not attempt to combine individual model - # bounding boxes. Get the forward transform and assign the bounding_box - # to it before evaluating it. The order Model.bounding_box is reversed. - transform.bounding_box = self.bounding_box - - return transform - - @property - def backward_transform(self): - """ - Return the total backward transform if available - from output to input - coordinate system. - - Raises - ------ - NotImplementedError : - An analytical inverse does not exist. - - """ - try: - backward = self.forward_transform.inverse - except NotImplementedError as err: - msg = f"Could not construct backward transform. \n{err}" - raise NotImplementedError(msg) from err - try: - _ = backward.inverse - except NotImplementedError: # means "hasattr" won't work - backward.inverse = self.forward_transform - return backward diff --git a/gwcs/wcs/_wcs.py b/gwcs/wcs/_wcs.py index aa5b60fb..bd64fbc3 100644 --- a/gwcs/wcs/_wcs.py +++ b/gwcs/wcs/_wcs.py @@ -1294,47 +1294,23 @@ def transform( 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 - from_step = self._get_step(from_frame) - to_step = self._get_step(to_frame) - transform = self.get_transform(from_step.step.frame, to_step.step.frame) + direction = self.pipeline_between(from_frame, to_frame) - if transform is None: + if direction is None: msg = f"No transformation found from {from_frame} to {to_frame}." raise ValueError(msg) - # Get the frame objects from the wcs pipeline - from_frame_obj = self.get_frame(from_frame) - to_frame_obj = self.get_frame(to_frame) - - 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, - ) + if direction.forward: + return direction.wcs.evaluate( + *args, + with_bounding_box=with_bounding_box, + fill_value=fill_value, + **kwargs, + ) - result = transform( + return direction.wcs.invert( *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 to_frame_obj is not None and to_frame_obj.naxes == 1: - return result[0] - return result @property def name(self) -> str: From 47d9d78086018aa01f119aacce01d700005f04cd Mon Sep 17 00:00:00 2001 From: William Jamieson Date: Mon, 2 Mar 2026 20:16:22 -0500 Subject: [PATCH 3/7] Remove no longer needed unit handling methods from WCS class --- gwcs/wcs/_wcs.py | 85 ------------------------------------------------ 1 file changed, 85 deletions(-) diff --git a/gwcs/wcs/_wcs.py b/gwcs/wcs/_wcs.py index bd64fbc3..29e4d4fe 100644 --- a/gwcs/wcs/_wcs.py +++ b/gwcs/wcs/_wcs.py @@ -345,91 +345,6 @@ def __call__( *args, with_bounding_box=with_bounding_box, fill_value=fill_value, **kwargs ) - def _units_are_present( - self, args: tuple[LowLevelInput, ...], transform: Model - ) -> tuple[bool, bool]: - """ - Determine 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) - if isinstance(transform, (Tabular1D, Tabular2D)): - transform_uses_quantity = ( - isinstance(transform.lookup_table, u.Quantity) - or transform.input_units_equivalencies is not None - ) - elif transform is not None and len(transform.parameters) == 0: - transform_uses_quantity = False - else: - 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: Model, - *args: LowLevelInput, - frame: CoordinateFrameProtocol, - input_is_quantity: bool = False, - transform_uses_quantity: bool = False, - **kwargs, - ) -> tuple[LowLevelInput, ...]: - """ - 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 - 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: - return frame.add_units(args) - if not transform_uses_quantity and input_is_quantity: - return frame.remove_units(args) - return args - - def _make_output_units_consistent( - self, - transform: Model, - *args: LowLevelInput, - frame: CoordinateFrameProtocol, - input_is_quantity=False, - transform_uses_quantity=False, - **kwargs, - ) -> tuple[LowLevelInput, ...]: - """ - 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 - - if input_is_quantity and transform_uses_quantity: - # make sure the output is returned in the units of the output frame - return frame.add_units(args) - if not input_is_quantity and ( - transform_uses_quantity or transform.parameters.size - ): - return frame.remove_units(args) - if not transform_uses_quantity and input_is_quantity: - return frame.add_units(args) - return args - def in_image( self, *args: LowLevelInput, From 1944b02e2691b061a10633e39c4f49e11edb9c6c Mon Sep 17 00:00:00 2001 From: William Jamieson Date: Mon, 2 Mar 2026 20:26:05 -0500 Subject: [PATCH 4/7] Centralize the high_level conversion API to use the same internal logic --- gwcs/utils.py | 28 ++-------------------------- gwcs/wcs/_wcs.py | 17 +++++++---------- 2 files changed, 9 insertions(+), 36 deletions(-) diff --git a/gwcs/utils.py b/gwcs/utils.py index 40a55f2d..65b8343e 100644 --- a/gwcs/utils.py +++ b/gwcs/utils.py @@ -495,34 +495,10 @@ def create_projection_transform(projcode): return projklass(**projparams) +# ToDo: Should this be deprecated? def is_high_level(*args, low_level_wcs): """ Determine if args matches the high level classes as defined by ``low_level_wcs``. """ - if low_level_wcs.world_axis_object_classes is None or len(args) != len( - low_level_wcs.world_axis_object_classes - ): - return False - - type_match = [ - (type(arg), waoc[0]) - for arg, waoc in zip( - args, low_level_wcs.world_axis_object_classes.values(), strict=False - ) - ] - - types_are_high_level = [argt is t for argt, t in type_match] - - if all(types_are_high_level): - return True - - if any(types_are_high_level): - msg = ( - "Invalid types were passed, got " - f"({', '.join(tm[0].__name__ for tm in type_match)}) expected " - f"({', '.join(tm[1].__name__ for tm in type_match)})." - ) - raise TypeError(msg) - - return False + return low_level_wcs.output_frame.is_high_level(*args) diff --git a/gwcs/wcs/_wcs.py b/gwcs/wcs/_wcs.py index 29e4d4fe..ab776739 100644 --- a/gwcs/wcs/_wcs.py +++ b/gwcs/wcs/_wcs.py @@ -23,10 +23,6 @@ ) from astropy.modeling.parameters import _tofloat from astropy.wcs.utils import celestial_frame_to_wcs, proj_plane_pixel_scales -from astropy.wcs.wcsapi.high_level_api import ( - high_level_objects_to_values, - values_to_high_level_objects, -) from numpy import typing as npt from scipy import optimize @@ -40,7 +36,7 @@ LowLevelInput, get_ctype_from_ucd, ) -from gwcs.utils import _compute_lon_pole, is_high_level, to_index +from gwcs.utils import _compute_lon_pole, to_index from ._exception import NoConvergence from ._pipeline import ForwardTransform, Pipeline @@ -472,11 +468,12 @@ def outside_footprint(self, world_arrays): axes_phys_types = self.world_axis_physical_types footprint = self.footprint() not_numerical = False - if is_high_level(world_arrays[0], low_level_wcs=self): + if self.output_frame.is_high_level(*world_arrays): not_numerical = True - world_arrays = high_level_objects_to_values( - *world_arrays, low_level_wcs=self + world_arrays = self.output_frame.from_high_level_coordinates( + *world_arrays, correct_1d=False ) + for axtyp in axes_types: ind = np.asarray(np.asarray(self.output_frame.axes_type) == axtyp) @@ -520,8 +517,8 @@ def outside_footprint(self, world_arrays): coord[outside] = np.nan world_arrays[idim] = coord if not_numerical: - world_arrays = values_to_high_level_objects( - *world_arrays, low_level_wcs=self + world_arrays = self.output_frame.to_high_level_coordinates( + *world_arrays, correct_1d=False ) return world_arrays From 98196cdf7750e215e26233b66e76ccfb81c23aea Mon Sep 17 00:00:00 2001 From: William Jamieson Date: Wed, 18 Mar 2026 14:20:37 -0400 Subject: [PATCH 5/7] Bugfix for ndcube Need to handle Identity directly rather than trying to infer it. --- gwcs/wcs/_wcs.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/gwcs/wcs/_wcs.py b/gwcs/wcs/_wcs.py index ab776739..6aa5186e 100644 --- a/gwcs/wcs/_wcs.py +++ b/gwcs/wcs/_wcs.py @@ -14,6 +14,7 @@ from astropy.modeling import Model, fix_inputs, projections from astropy.modeling.bounding_box import ModelBoundingBox as Bbox from astropy.modeling.models import ( + Identity, Mapping, RotateCelestial2Native, Shift, @@ -132,7 +133,7 @@ def __init__( or transform.input_units_equivalencies is not None ) - elif transform is not None and len(transform.parameters) == 0: + elif isinstance(transform, Identity): transform_uses_quantity = False else: From aea34fbcdcbb39839018594e9dc351974e153e97 Mon Sep 17 00:00:00 2001 From: William Jamieson Date: Mon, 9 Feb 2026 16:31:26 -0500 Subject: [PATCH 6/7] Add change log --- changes/696.feature.rst | 8 ++++++++ docs/conf.py | 2 -- 2 files changed, 8 insertions(+), 2 deletions(-) create mode 100644 changes/696.feature.rst diff --git a/changes/696.feature.rst b/changes/696.feature.rst new file mode 100644 index 00000000..852c9994 --- /dev/null +++ b/changes/696.feature.rst @@ -0,0 +1,8 @@ +Clean up internal code flow so that the entire GWCS Native API follows common code +paths. This enables the Native API to properly handle various forms of inputs in +a consistent manner so that evaluation and inverse have the same input handling +behavior. This means HLO and quantities should function properly as they are passed +into these portions of the Native API. + +Also allows for quantities to be passed into the "low-level APE 14 API" without +raising errors. diff --git a/docs/conf.py b/docs/conf.py index 53ce0d7e..35aa488c 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -105,8 +105,6 @@ # Enable nitpicky mode - which ensures that all references in the docs resolve. nitpicky = True nitpick_ignore = [ - ("py:class", "gwcs.api.GWCSAPIMixin"), - ("py:class", "gwcs.wcs._pipeline.Pipeline"), ("py:obj", "astropy.modeling.projections.projcodes"), ("py:attr", "gwcs.WCS.bounding_box"), ("py:meth", "gwcs.WCS.footprint"), From 443ddddbdfe7af3be3b4aabfc9cebbb8e228f045 Mon Sep 17 00:00:00 2001 From: William Jamieson Date: Thu, 16 Jul 2026 15:28:31 -0400 Subject: [PATCH 7/7] Address automated review tool findings Fixes for 1. Inconsistencies in the API. 2. Backwards compatibility 3. Broken edge cases 4. Deprecation warnings for legacy code --- gwcs/coordinate_frames/__init__.py | 3 + gwcs/coordinate_frames/_base.py | 143 ++++---- gwcs/coordinate_frames/_empty.py | 4 +- gwcs/tests/test_api.py | 8 +- gwcs/tests/test_coordinate_systems.py | 41 ++- gwcs/tests/test_utils.py | 8 + gwcs/tests/test_wcs.py | 353 ++++++++++++++++++ gwcs/utils.py | 27 +- gwcs/wcs/__init__.py | 5 +- gwcs/wcs/_pipeline.py | 496 ++++++++++++++++---------- gwcs/wcs/_step.py | 62 +++- gwcs/wcs/_wcs.py | 126 +++++-- gwcs/wcstools.py | 3 +- 13 files changed, 1000 insertions(+), 279 deletions(-) diff --git a/gwcs/coordinate_frames/__init__.py b/gwcs/coordinate_frames/__init__.py index 2bc59dab..5033a70f 100644 --- a/gwcs/coordinate_frames/__init__.py +++ b/gwcs/coordinate_frames/__init__.py @@ -218,6 +218,7 @@ LowLevelInput, WorldAxisObjectClasses, WorldAxisObjectComponent, + _LegacyCoordinateFrameProtocol, ) from ._celestial import CelestialFrame from ._composite import CompositeFrame @@ -247,5 +248,7 @@ "TemporalFrame", "WorldAxisObjectClasses", "WorldAxisObjectComponent", + # See comment in _base.py about why this is included in __all__. + "_LegacyCoordinateFrameProtocol", "get_ctype_from_ucd", ] diff --git a/gwcs/coordinate_frames/_base.py b/gwcs/coordinate_frames/_base.py index 5193883c..c4b70390 100644 --- a/gwcs/coordinate_frames/_base.py +++ b/gwcs/coordinate_frames/_base.py @@ -24,6 +24,8 @@ values_to_high_level_objects, ) +from gwcs.utils import correct_1d_output + from ._axis import AxesType __all__ = [ @@ -35,6 +37,11 @@ "WorldAxisObjectClass", "WorldAxisObjectClassConverter", "WorldAxisObjectComponent", + # Adding _LegacyCoordinateFrameProtocol to __all__ in order for it to be + # picked up by Sphinx for the API documentation. This should be removed + # when the deprecation period of coordinate frames missing is_high_level + # is over. + "_LegacyCoordinateFrameProtocol", ] _DtypeGeneric = TypeVar("_DtypeGeneric", bound=np.generic) @@ -142,9 +149,10 @@ def from_tuple(cls, tup: tuple[str, str | int, str]) -> Self: @runtime_checkable -class CoordinateFrameProtocol(Protocol): +class _LegacyCoordinateFrameProtocol(Protocol): """ - API Definition for a Coordinate frame + Original API definition for a coordinate frame. This is used for deprecation + warnings and to identify when to patch the ``is_high_level`` method to the frame. """ @property @@ -282,56 +290,7 @@ def remove_units( for array in self.add_units(arrays) ) - def is_high_level(self, *args) -> bool: - """ - Return `True` if the input coordinates are already high level objects - described by this frame. - - This is used by the low level WCS API in Astropy to determine whether - to call ``to_high_level_coordinates`` or not. - """ - - if (world_axis_object_classes := self.world_axis_object_classes) is None or len( - args - ) != len(world_axis_object_classes): - return False - - type_match = [] - for arg, world_axis_object_class in zip( - args, world_axis_object_classes.values(), strict=True - ): - if isinstance(class_object := world_axis_object_class.class_object, str): - type_match.append( - type(arg).__name__ == class_object - and class_object != u.Quantity.__name__ - ) - else: - type_match.append( - isinstance(arg, class_object) and class_object is not u.Quantity - ) - - if all(type_match): - return True - - if any(type_match): - types = [ - ( - type(arg).__name__, - c.class_object - if isinstance(c.class_object, str) - else c.class_object.__name__, - ) - for arg, c in zip(args, world_axis_object_classes.values(), strict=True) - ] - msg = ( - "Invalid types were passed, got " - f"({', '.join(t[0] for t in types)}), but expected " - f"({', '.join(t[1] for t in types)})." - ) - raise TypeError(msg) - - return False - + @correct_1d_output def to_high_level_coordinates(self, *values, correct_1d=True): """ Convert "values" to high level coordinate objects described by this frame. @@ -358,11 +317,9 @@ def to_high_level_coordinates(self, *values, correct_1d=True): msg = "All values should be a scalar number or a numpy array." raise TypeError(msg) - high_level = values_to_high_level_objects(*values, low_level_wcs=self) - if correct_1d and len(high_level) == 1: - high_level = high_level[0] - return high_level + return values_to_high_level_objects(*values, low_level_wcs=self) + @correct_1d_output def from_high_level_coordinates(self, *high_level_coords, correct_1d=True): """ Convert high level coordinate objects to "values" as described by this frame. @@ -381,10 +338,76 @@ def from_high_level_coordinates(self, *high_level_coords, correct_1d=True): values : `numbers.Number` or `numpy.ndarray` ``naxis`` number of coordinates as scalars or arrays. """ - values = high_level_objects_to_values(*high_level_coords, low_level_wcs=self) - if correct_1d and self.naxes == 1: - values = values[0] - return values + return high_level_objects_to_values(*high_level_coords, low_level_wcs=self) + + +def _is_high_level(self: CoordinateFrameProtocol, *args) -> bool: + """ + Return `True` if the input coordinates are already high level objects + described by this frame. + + This is used by the low level WCS API in Astropy to determine whether + to call ``to_high_level_coordinates`` or not. + """ + + if (world_axis_object_classes := self.world_axis_object_classes) is None or len( + args + ) != len(world_axis_object_classes): + return False + + type_match = [] + for arg, world_axis_object_class in zip( + args, world_axis_object_classes.values(), strict=True + ): + if isinstance(class_object := world_axis_object_class.class_object, str): + type_match.append( + type(arg).__name__ == class_object + and class_object != u.Quantity.__name__ + ) + else: + type_match.append( + isinstance(arg, class_object) and class_object is not u.Quantity + ) + + if all(type_match): + return True + + if any(type_match): + types = [ + ( + type(arg).__name__, + c.class_object + if isinstance(c.class_object, str) + else c.class_object.__name__, + ) + for arg, c in zip(args, world_axis_object_classes.values(), strict=True) + ] + msg = ( + "Invalid types were passed, got " + f"({', '.join(t[0] for t in types)}), but expected " + f"({', '.join(t[1] for t in types)})." + ) + raise TypeError(msg) + + return False + + +@runtime_checkable +class CoordinateFrameProtocol(_LegacyCoordinateFrameProtocol, Protocol): + """ + API Definition for a Coordinate frame + """ + + def is_high_level(self, *args) -> bool: + """ + Return `True` if the input coordinates are already high level objects + described by this frame. + + This is used by the low level WCS API in Astropy to determine whether + to call ``to_high_level_coordinates`` or not. + """ + + return _is_high_level(self, *args) class BaseCoordinateFrame(CoordinateFrameProtocol): diff --git a/gwcs/coordinate_frames/_empty.py b/gwcs/coordinate_frames/_empty.py index a3c1ade1..3d42a165 100644 --- a/gwcs/coordinate_frames/_empty.py +++ b/gwcs/coordinate_frames/_empty.py @@ -144,10 +144,10 @@ def world_axis_object_components(self) -> list[WorldAxisObjectComponent]: for i, at in enumerate(self.axes_type) ] - def to_high_level_coordinates(self, *values): + def to_high_level_coordinates(self, *values, correct_1d: bool = True): self._raise_error() - def from_high_level_coordinates(self, *high_level_coords): + def from_high_level_coordinates(self, *high_level_coords, correct_1d: bool = True): self._raise_error() def add_units( diff --git a/gwcs/tests/test_api.py b/gwcs/tests/test_api.py index 2121e6a1..614102d7 100644 --- a/gwcs/tests/test_api.py +++ b/gwcs/tests/test_api.py @@ -597,8 +597,12 @@ 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) + # invert (native API) must accept the high-level/quantity world value and + # round-trip back to the original pixel, not merely match a hardcoded zero. + input_pixel = 3.0 + world_nonzero = wcs.pixel_to_world(input_pixel) + pixel2 = wcs.invert(world_nonzero) + assert u.allclose(pixel2, input_pixel * u.pix) def test_world_axis_object_components_units(gwcs_3d_identity_units): diff --git a/gwcs/tests/test_coordinate_systems.py b/gwcs/tests/test_coordinate_systems.py index 65143515..60edb719 100644 --- a/gwcs/tests/test_coordinate_systems.py +++ b/gwcs/tests/test_coordinate_systems.py @@ -16,6 +16,7 @@ from gwcs import WCS from gwcs import coordinate_frames as cf +from gwcs.coordinate_frames._base import _LegacyCoordinateFrameProtocol from gwcs.coordinate_frames._utils import ( _ALLOWED_UCD_DUPLICATES, _ucd1_to_ctype_name_mapping, @@ -649,14 +650,50 @@ def remove_units(self, arrays): def is_high_level(self, *inputs): return False - def to_high_level_coordinates(self, *inputs): + def to_high_level_coordinates(self, *inputs, correct_1d=True): + if correct_1d and len(inputs) == 1: + return inputs[0] return inputs - def from_high_level_coordinates(self, *inputs): + def from_high_level_coordinates(self, *inputs, correct_1d=True): + if correct_1d and len(inputs) == 1: + return inputs[0] return inputs frame = MyFrame() assert isinstance(frame, cf.CoordinateFrameProtocol) + assert frame.to_high_level_coordinates(1, correct_1d=False) == (1,) + assert frame.from_high_level_coordinates(1, correct_1d=False) == (1,) + + +def test_LegacyCoordinateFrameProtocol(): + class LegacyFrame: + naxes = 1 + name = "legacy" + unit = (u.m,) + axes_names = ("x",) + axes_order = (0,) + reference_frame = None + axes_type = ("SPATIAL",) + axis_physical_types = ("custom:x",) + world_axis_object_classes = (u.Quantity,) + world_axis_object_components = ("custom:x",) + + def add_units(self, arrays): + return arrays + + def remove_units(self, arrays): + return arrays + + def to_high_level_coordinates(self, *inputs): + return inputs + + def from_high_level_coordinates(self, *inputs): + return inputs + + frame = LegacyFrame() + assert isinstance(frame, _LegacyCoordinateFrameProtocol) + assert not isinstance(frame, cf.CoordinateFrameProtocol) def test_BaseCoordinateFrame(): diff --git a/gwcs/tests/test_utils.py b/gwcs/tests/test_utils.py index fdc62729..7b7ddec9 100644 --- a/gwcs/tests/test_utils.py +++ b/gwcs/tests/test_utils.py @@ -136,3 +136,11 @@ def test_get_values(): res = gwutils.get_values(None, args) assert res == [2] + + +def test_is_high_level_deprecated(gwcs_simple_2d): + """Test that the deprecated `is_high_level` function works and issues a warning.""" + with pytest.warns( + DeprecationWarning, match=r"The use of `is_high_level` is deprecated" + ): + assert gwutils.is_high_level(1, 2, low_level_wcs=gwcs_simple_2d) is False diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py index 5ad7ae63..48a8baf1 100644 --- a/gwcs/tests/test_wcs.py +++ b/gwcs/tests/test_wcs.py @@ -57,6 +57,70 @@ pipe = [wcs.Step(detector, m1), wcs.Step(focal, m2), wcs.Step(icrs, None)] pipe_copy = pipe.copy() + +class LegacyFrameAdapter: + """ + Duck-typed adapter that exposes the historical coordinate-frame API + (without ``is_high_level``) by delegating to a wrapped frame. Used to + exercise the legacy-frame compatibility paths. + """ + + def __init__(self, frame): + self._frame = frame + + @property + def naxes(self): + return self._frame.naxes + + @property + def name(self): + return self._frame.name + + @property + def unit(self): + return self._frame.unit + + @property + def axes_names(self): + return self._frame.axes_names + + @property + def axes_order(self): + return self._frame.axes_order + + @property + def reference_frame(self): + return self._frame.reference_frame + + @property + def axes_type(self): + return self._frame.axes_type + + @property + def axis_physical_types(self): + return self._frame.axis_physical_types + + @property + def world_axis_object_classes(self): + return self._frame.world_axis_object_classes + + @property + def world_axis_object_components(self): + return self._frame.world_axis_object_components + + def add_units(self, arrays): + return self._frame.add_units(arrays) + + def remove_units(self, arrays): + return self._frame.remove_units(arrays) + + def to_high_level_coordinates(self, *values): + return self._frame.to_high_level_coordinates(*values) + + def from_high_level_coordinates(self, *coords): + return self._frame.from_high_level_coordinates(*coords) + + # Create some data. nx, ny = (5, 2) x = np.linspace(0, 1, nx) @@ -1944,6 +2008,295 @@ def test_parameterless_transform(): assert gwcs.invert(1 * u.pix, 1 * u.pix) == (1 * u.pix, 1 * u.pix) +def test_parameterless_transform_mapping(): + """ + Test that a transform with no parameters other than ``Identity`` correctly + handles units. + + ``Mapping`` also has zero parameters (like ``Identity``), but is not given + the same special-cased treatment as ``Identity`` when determining whether + the transform "uses quantity". ``Model.uses_quantity`` returns `True` for + *any* model with zero parameters (regardless of whether it actually cares + about units), so a ``Mapping``-only step must not be allowed to have units + silently added to its inputs. + + Regression test for #558 (same root cause as ``test_parameterless_transform``, + but for a different zero-parameter model). + """ + + in_frame = cf.Frame2D(name="in_frame", unit=(u.pix, u.pix)) + out_frame = cf.Frame2D(name="out_frame", unit=(u.pix, u.pix)) + + gwcs = wcs.WCS( + [ + (in_frame, models.Mapping((1, 0))), + (out_frame, None), + ] + ) + + # No units introduced by the forward transform + assert gwcs(1, 1) == (1, 1) + assert gwcs(1 * u.pix, 1 * u.pix) == (1 * u.pix, 1 * u.pix) + + # No units introduced by the inverse transform + assert gwcs.invert(1, 1) == (1, 1) + + # Unitless evaluation must still work when a bounding_box is present: + # incorrectly wrapping unitless inputs in Quantity before calling the + # transform causes bounding-box comparisons to raise UnitConversionError. + gwcs.bounding_box = ((0, 10), (0, 10)) + assert np.isnan(gwcs(20, 5)).all() + assert gwcs(5, 5) == (5, 5) + + +def test_invert_with_legacy_frame_conversion_signatures(monkeypatch): + """ + Regression test for legacy custom frames where conversion methods do not + accept the ``correct_1d`` keyword argument. + """ + + gwcs = wcs.WCS([(detector, m), (icrs, None)]) + world = gwcs.pixel_to_world(1, 2) + expected_pixel = gwcs.world_to_pixel(world) + + original_from = gwcs.output_frame.from_high_level_coordinates + original_to = gwcs.input_frame.to_high_level_coordinates + + # Legacy signatures accepted only positional coordinates. + def legacy_from_high_level_coordinates(*coords): + return original_from(*coords) + + def legacy_to_high_level_coordinates(*coords): + return original_to(*coords) + + monkeypatch.setattr( + gwcs.output_frame, + "from_high_level_coordinates", + legacy_from_high_level_coordinates, + ) + monkeypatch.setattr( + gwcs.input_frame, + "to_high_level_coordinates", + legacy_to_high_level_coordinates, + ) + + pixel = gwcs.invert(world) + pixel = gwcs.input_frame.remove_units(pixel) + if not isinstance(pixel, list | tuple): + pixel = (pixel,) + + assert_allclose(pixel, expected_pixel) + + +def test_invert_with_legacy_frame_1d_signatures(monkeypatch): + """ + 1D legacy frame whose conversion methods lack `correct_1d`: invert must + not double-unwrap the single output. Regression for correct_1d being + dropped in the legacy shim. + """ + in_frame = cf.CoordinateFrame( + naxes=1, + axes_type="SPATIAL", + axes_order=(0,), + unit=(u.pix,), + axes_names=("x",), + name="in1d", + ) + out_frame = cf.SpectralFrame(unit=u.Hz, axes_order=(0,), name="out1d") + gw = wcs.WCS([(in_frame, models.Scale(2.0)), (out_frame, None)]) + + world = gw.pixel_to_world(3) + expected_pixel = gw.world_to_pixel(world) + + original_from = gw.output_frame.from_high_level_coordinates + original_to = gw.input_frame.to_high_level_coordinates + + def func_from(*args): + return original_from(*args) + + def func_to(*args): + return original_to(*args) + + # Legacy signatures: positional only, no correct_1d. + monkeypatch.setattr(gw.output_frame, "from_high_level_coordinates", func_from) + monkeypatch.setattr(gw.input_frame, "to_high_level_coordinates", func_to) + + pixel = gw.invert(world) + pixel = gw.input_frame.remove_units(pixel) + assert_allclose(pixel, expected_pixel) + + +def test_transform_same_frame_raises_value_error(): + """ + ``transform`` between a frame and itself has no pipeline (``pipeline_between`` + returns ``None``) and must raise a ``ValueError``. + """ + gw = wcs.WCS([wcs.Step(detector, m1), wcs.Step(focal, m2), wcs.Step(icrs, None)]) + + with pytest.raises(ValueError, match=r"No transformation found from"): + gw.transform("focal", "focal", 1, 2) + + +def test_construction_from_base_pipeline_skips_validation(monkeypatch): + """ + Constructing a ``WCS``/``Pipeline`` from a `~gwcs.wcs.BasePipeline` must take + the no-copy fast path: it must NOT call ``_initialize_pipeline`` and must + share the underlying ``Step`` objects by identity with the source pipeline. + + This locks in the performance contract of the ``pipeline_between`` rework so + a future change cannot silently reintroduce per-call validation/copying on + the hot evaluation path. + """ + gw = wcs.WCS([wcs.Step(detector, m1), wcs.Step(focal, m2), wcs.Step(icrs, None)]) + + base = wcs.BasePipeline(gw.pipeline) + + called = False + original_initialize = wcs.Pipeline._initialize_pipeline + + def _tracking_initialize(self, *args, **kwargs): + nonlocal called + called = True + return original_initialize(self, *args, **kwargs) + + monkeypatch.setattr(wcs.Pipeline, "_initialize_pipeline", _tracking_initialize) + + fast = wcs.WCS(base) + + # The fast path must not re-run pipeline initialization/validation. + assert called is False + + # The steps must be shared by identity, not copied. + assert len(fast.pipeline) == len(gw.pipeline) + for fast_step, src_step in zip(fast.pipeline, gw.pipeline, strict=True): + assert fast_step is src_step + + # Sanity check: the fast-path WCS still evaluates correctly. + assert_allclose(fast(1, 2), gw(1, 2)) + + +def test_transform_does_not_corrupt_parent_pipeline(): + """ + ``transform``/``get_transform`` build ephemeral sub-pipelines that share the + parent's ``Step`` objects (documented no-copy behavior). Because these + operations are read-only, repeated calls must leave the parent's own + evaluation results unchanged. + """ + gw = wcs.WCS([wcs.Step(detector, m1), wcs.Step(focal, m2), wcs.Step(icrs, None)]) + + expected_forward = gw(1, 2) + expected_inverse = gw.invert(*expected_forward) + + for _ in range(5): + gw.transform("detector", "icrs", 1, 2) + gw.transform("icrs", "detector", *expected_forward) + gw.get_transform("detector", "focal") + + assert_allclose(gw(1, 2), expected_forward) + assert_allclose(gw.invert(*expected_forward), expected_inverse) + + +def test_correct_1d_output_decorator(): + """ + Direct coverage of the ``correct_1d_output`` decorator: single-element + sequences collapse to a scalar only when ``correct_1d`` is True, longer + sequences pass through, and non-sequence returns are never indexed. + """ + from gwcs.utils import correct_1d_output + + def returns_one(*args): + return (args[0],) + + def returns_two(*args): + return (args[0], args[1]) + + def returns_scalar(*args): + return args[0] + + # pass_correct_1d=False path (legacy-style functions without the kwarg) + collapse = correct_1d_output(returns_one, pass_correct_1d=False) + assert collapse(42, correct_1d=True) == 42 + assert collapse(42, correct_1d=False) == (42,) + + passthrough = correct_1d_output(returns_two, pass_correct_1d=False) + assert passthrough(1, 2, correct_1d=True) == (1, 2) + + # A non-sequence return must never be indexed, regardless of correct_1d. + scalar = correct_1d_output(returns_scalar, pass_correct_1d=False) + assert scalar(7, correct_1d=True) == 7 + + # pass_correct_1d=True path forwards correct_1d into the wrapped function. + def uses_correct_1d(*args, correct_1d=True): + return (args[0], correct_1d) + + forwarded = correct_1d_output(uses_correct_1d, pass_correct_1d=True) + # Two-element result so it is not collapsed; confirms correct_1d was passed. + assert forwarded(9, correct_1d=False) == (9, False) + + +def test_legacy_coordinate_frame_protocol_warns_and_functions(): + """ + Legacy frames that do not define ``is_high_level`` should still work, + while issuing a deprecation warning. + """ + + legacy_input = LegacyFrameAdapter(cf.Frame2D(name="legacy_detector")) + assert not hasattr(legacy_input, "is_high_level") + + with pytest.warns(DeprecationWarning, match=r"is_high_level"): + gw = wcs.WCS([(legacy_input, m), (icrs, None)]) + + # Check that legacy_input has not been mutated to have is_high_level + assert not hasattr(legacy_input, "is_high_level") + + # Frame should be patched for compatibility and remain callable. + assert hasattr(gw.input_frame, "is_high_level") + assert gw.input_frame.is_high_level(1, 2) is False + assert_allclose(gw(1, 2), m(1, 2)) + + legacy_output = LegacyFrameAdapter(icrs) + with pytest.warns(DeprecationWarning, match=r"is_high_level"): + gw_legacy_output = wcs.WCS([(detector, m), (legacy_output, None)]) + + world = gw_legacy_output.pixel_to_world(1, 2) + assert gw_legacy_output.output_frame.is_high_level(world) is True + + +def test_standard_coordinate_frames_do_not_warn_as_legacy(): + with warnings.catch_warnings(record=True) as captured: + warnings.simplefilter("always", DeprecationWarning) + _ = wcs.WCS([(detector, m), (icrs, None)]) + + assert not any("is_high_level" in str(w.message) for w in captured) + + +def test_legacy_frame_validation_does_not_recheck_after_patch(monkeypatch): + """ + Regression test for Python 3.11 legacy-frame flow where legacy detection + should not be recomputed after patching ``is_high_level`` onto the frame. + """ + + from gwcs.wcs import _step as step_mod + + legacy_input = LegacyFrameAdapter(cf.Frame2D(name="legacy_detector_recheck")) + + # Emulate the 3.11-style behavior where legacy status depends on the + # presence of is_high_level and nominal coordinate-frame checks fail. + monkeypatch.setattr(step_mod, "_is_coordinate_frame", lambda frame: False) + + def _legacy_check(frame): + return hasattr(frame, "to_high_level_coordinates") and not hasattr( + frame, "is_high_level" + ) + + monkeypatch.setattr(step_mod, "_is_legacy_coordinate_frame", _legacy_check) + + with pytest.warns(DeprecationWarning, match=r"is_high_level"): + step = wcs.Step(legacy_input, m) + + assert hasattr(step.frame, "is_high_level") + + def test_fitswcs_imaging(fits_wcs_imaging_simple): """Test simple FITS type imaging WCS.""" gwcs, astwcs = fits_wcs_imaging_simple diff --git a/gwcs/utils.py b/gwcs/utils.py index 65b8343e..345aa4ba 100644 --- a/gwcs/utils.py +++ b/gwcs/utils.py @@ -495,10 +495,35 @@ def create_projection_transform(projcode): return projklass(**projparams) -# ToDo: Should this be deprecated? +def correct_1d_output(func, pass_correct_1d: bool = True): + """ + Decorator to correct the output of a function to be 1D if the input is 1D. + """ + + @functools.wraps(func) + def wrapper(*args, correct_1d: bool = True, **kwargs): + if pass_correct_1d: + value = func(*args, correct_1d=correct_1d, **kwargs) + else: + value = func(*args, **kwargs) + + if correct_1d and isinstance(value, (tuple, list)) and len(value) == 1: + return value[0] + + return value + + return wrapper + + def is_high_level(*args, low_level_wcs): """ Determine if args matches the high level classes as defined by ``low_level_wcs``. """ + warnings.warn( + "The use of `is_high_level` is deprecated. Use the `is_high_level` method " + "of the low level WCS object's output_frame instead.", + DeprecationWarning, + stacklevel=2, + ) return low_level_wcs.output_frame.is_high_level(*args) diff --git a/gwcs/wcs/__init__.py b/gwcs/wcs/__init__.py index 9bf22ced..782e1a5c 100644 --- a/gwcs/wcs/__init__.py +++ b/gwcs/wcs/__init__.py @@ -1,11 +1,12 @@ from ._exception import GwcsBoundingBoxWarning, NoConvergence -from ._pipeline import DirectionalPipeline, Pipeline +from ._pipeline import BasePipeline, DirectionalWCS, Pipeline from ._step import Step from ._wcs import WCS __all__ = [ "WCS", - "DirectionalPipeline", + "BasePipeline", + "DirectionalWCS", "GwcsBoundingBoxWarning", "NoConvergence", "Pipeline", diff --git a/gwcs/wcs/_pipeline.py b/gwcs/wcs/_pipeline.py index 86ec8991..17eaf83e 100644 --- a/gwcs/wcs/_pipeline.py +++ b/gwcs/wcs/_pipeline.py @@ -2,7 +2,8 @@ import warnings from collections.abc import Sequence -from typing import Generic, NamedTuple, Self, TypeAlias, TypeVar, Union, overload +from dataclasses import dataclass +from typing import Generic, Self, TypeAlias, TypeVar, Union, overload from astropy.modeling import Model from astropy.modeling.bounding_box import CompoundBoundingBox, ModelBoundingBox @@ -14,30 +15,321 @@ from ._exception import GwcsBoundingBoxWarning, GwcsFrameExistsError from ._step import IndexedStep, Mdl, Step, StepTuple -__all__ = ["DirectionalPipeline", "ForwardTransform", "Pipeline"] +__all__ = ["BasePipeline", "DirectionalWCS", "ForwardTransform", "Pipeline"] + + +class BasePipeline: + """ + Base class for the Pipeline functionality. + + This class is intended to be an internal implementation detail for the + Pipeline class that is used to handle pipeline functionality without + the overhead of the checks and protections that are in place in the + pipeline class. + + This class only has the minimum functionality to handle the building of the + transforms. + """ + + _pipeline: list[Step] + + def __init__(self, pipeline: list[Step] | BasePipeline | None = None) -> None: + if isinstance(pipeline, BasePipeline): + self._pipeline = pipeline.pipeline + else: + self._pipeline = pipeline or [] + + @property + def pipeline(self) -> list[Step]: + """ + Allow direct access to the raw pipeline steps. + """ + + # TODO: This can still allow direct modification of the pipeline list + # without any of the checks and handling that have been put in + # place in order to ensure the pipeline is functional. + # -> Maybe we should return a copy? + return self._pipeline + + @property + def available_frames(self) -> list[str]: + """ + List of all the frame names in this WCS in their order in the pipeline + """ + return [step.frame.name for step in self._pipeline] + + @staticmethod + def _frame_name(frame: str | CoordinateFrameProtocol) -> str: + """ + Return the name of the frame. + + Parameters + ---------- + frame : str or `~gwcs.coordinate_frames.CoordinateFrameProtocol` + Name of the frame or the frame object. + + Returns + ------- + Name of the frame. + """ + return frame.name if isinstance(frame, CoordinateFrameProtocol) else frame + + def _frame_index(self, frame: str | CoordinateFrameProtocol) -> int: + """ + Return the index of the given frame in the pipeline. + + Parameters + ---------- + frame : str or `~gwcs.coordinate_frames.CoordinateFrameProtocol` + Name of the frame or the frame object. + + Returns + ------- + Index of the frame in the pipeline. + """ + try: + return self.available_frames.index(self._frame_name(frame)) + except ValueError as err: + msg = f"Frame {self._frame_name(frame)} is not in the available frames" + raise CoordinateFrameError(msg) from err + + @property + def bounding_box(self) -> ModelBoundingBox | CompoundBoundingBox | None: + """ + Return the bounding box of the pipeline. + """ + # Pull the first transform of the pipeline which is what controls the + # bounding_box + transform = self._pipeline[0].transform + + if transform is None: + return None + + try: + return transform.bounding_box + except NotImplementedError: + return None + + @property + def forward_transform(self) -> Model: + """ + Return the forward transform of the pipeline. + """ + transform = combine_transforms([step.transform for step in self._pipeline[:-1]]) + + if self.bounding_box is not None: + # Currently compound models do not attempt to combine individual model + # bounding boxes. Get the forward transform and assign the bounding_box + # to it before evaluating it. The order Model.bounding_box is reversed. + transform.bounding_box = self.bounding_box + + return transform + + @property + def backward_transform(self) -> Model: + """ + Return the total backward transform if available - from output to input + coordinate system. + + Raises + ------ + NotImplementedError : + An analytical inverse does not exist. + + """ + try: + backward = self.forward_transform.inverse + except NotImplementedError: + # Try to get the backward transform by combining the inverses of the + # individual steps. This will work even if the forward transform does + # not have an analytical inverse, as long as each step does. + backward = combine_transforms( + [step.inverse for step in self._pipeline[:-1][::-1]] + ) + + try: + _ = backward.inverse + except NotImplementedError: # means "hasattr" won't work + backward.inverse = self.forward_transform + return backward + + def _pipeline_between( + self, + from_frame: str | CoordinateFrameProtocol, + to_frame: str | CoordinateFrameProtocol, + ) -> DirectionalWCS[BasePipeline] | None: + """ + Return a pipeline between the two given frames. + + This is the internal fast path: it builds an unvalidated + `~gwcs.wcs.BasePipeline` (no frame validation or step copies). Callers + that only need the combined transform (e.g. ``get_transform``) use this + to avoid the overhead of constructing a full `~gwcs.wcs.Pipeline` / + `~gwcs.wcs.WCS`. The public `pipeline_between` wraps the result in an + instance of the calling type. + + Validation is unnecessary here because the steps are sliced from an + already-valid pipeline (only a terminal ``Step(frame, None)`` is added). + The shared ``Step`` objects are safe because the returned pipeline is + ephemeral: callers use its built-in functionality (``forward_transform`` + etc.) rather than mutating its raw list of steps. + + Parameters + ---------- + from_frame : str or `~gwcs.coordinate_frames.CoordinateFrameProtocol` + Initial coordinate frame name of object. + to_frame : str or `~gwcs.coordinate_frames.CoordinateFrameProtocol` + End coordinate frame name or object. + + Returns + ------- + pipeline : `~gwcs.wcs.DirectionalWCS` or None + A ``DirectionalWCS`` (wcs, forward) where ``wcs`` is a + `~gwcs.wcs.BasePipeline` for the steps between the two frames and + ``forward`` indicates the direction (True: from_frame -> to_frame, + False: to_frame -> from_frame). Note that if from_frame and to_frame + are the same, then None is returned. + """ + from_index = self._frame_index(from_frame) + to_index = self._frame_index(to_frame) + + # Moving backwards over the pipeline + if to_index < from_index: + pipeline = self._pipeline[to_index:from_index] + pipeline.append(Step(self._pipeline[from_index].frame, None)) + + return DirectionalWCS(wcs=BasePipeline(pipeline), forward=False) + + if from_index < to_index: + pipeline = self._pipeline[from_index:to_index] + pipeline.append(Step(self._pipeline[to_index].frame, None)) + + return DirectionalWCS(wcs=BasePipeline(pipeline), forward=True) + + return None # from and to are the same frame, so no pipeline needed + + def pipeline_between( + self, + from_frame: str | CoordinateFrameProtocol, + to_frame: str | CoordinateFrameProtocol, + ) -> DirectionalWCS[Self] | None: + """ + Return a pipeline between the two given frames. + + Parameters + ---------- + from_frame : str or `~gwcs.coordinate_frames.CoordinateFrame` + Initial coordinate frame name of object. + to_frame : str or `~gwcs.coordinate_frames.CoordinateFrame` + End coordinate frame name or object. + + Returns + ------- + pipeline : + A `~gwcs.wcs.DirectionalWCS` (wcs, forward) where wcs is of the type + of the calling object. The forward attribute indicates if the + wcs/pipeline between the two frames is forward (True, from_frame -> + to_frame) or backward (False, to_frame -> from_frame). + Note that if from_frame and to_frame are the same, then None is + returned. + """ + direction = self._pipeline_between(from_frame, to_frame) + + if direction is None: + return None + + return DirectionalWCS(wcs=type(self)(direction.wcs), forward=direction.forward) + + def get_transform( + self, + from_frame: str | CoordinateFrameProtocol, + to_frame: str | CoordinateFrameProtocol, + ) -> Mdl: + """ + Return a transform between two coordinate frames. + + Parameters + ---------- + from_frame : str or `~gwcs.coordinate_frames.CoordinateFrameProtocol` + Initial coordinate frame name of object. + to_frame : str or `~gwcs.coordinate_frames.CoordinateFrameProtocol` + End coordinate frame name or object. + + Returns + ------- + transform : `~astropy.modeling.Model` + Transform between two frames. + """ + direction = self._pipeline_between(from_frame, to_frame) + + # from and to are the same frame, so no transform needed + if direction is None: + return None + + # If the pipeline is forward, return the forward transform + if direction.forward: + return direction.wcs.forward_transform + + # Otherwise it is backward, so return the backward transform + return direction.wcs.backward_transform + # Type aliases due to the use of the `|` for type hints not working with Model -ForwardTransform: TypeAlias = Union[Model, Sequence[Step | StepTuple]] # noqa: UP007 +ForwardTransform: TypeAlias = Union[Model, Sequence[Step | StepTuple] | BasePipeline] # noqa: UP007 + -_T = TypeVar("_T") +_T = TypeVar("_T", bound=BasePipeline) -class DirectionalPipeline(NamedTuple, Generic[_T]): +@dataclass(frozen=True, slots=True) +class DirectionalWCS(Generic[_T]): """ - A named tuple to hold a pipeline and the direction it is going in. + Dataclass to hold the WCS and the direction of the WCS's pipeline between + two frames. + + Note that the public use of this class is intended only for use on the full + `~gwcs.wcs.WCS` class, which is a subclass of `~gwcs.wcs.Pipeline`. + + Attributes + ---------- + wcs : _T + The WCS object that represents the pipeline between two frames. + + forward : bool + A boolean indicating if the pipeline is forward, True, (from_frame to to_frame) + or backward (to_frame to from_frame), False. """ wcs: _T forward: bool -class Pipeline: +class Pipeline(BasePipeline): """ Class to handle a sequence of WCS transformations. - This is intended to act line a list of steps, but with built in protections + This is intended to act like a list of steps, but with built in protections for things like duplicate frames. In addition, this handles all the logic for handling steps and their frames/transforms. + + This class is intended to be an organizational class for the WCS which handles + the sequence of steps and their transforms. It is not intended to be used + directly by outside users, but rather through the `~gwcs.wcs.WCS` class, + which subclasses this class and adds the additional evaluation functionality + and user interface for the WCS. + + .. warning:: + + This class supports passing any `~gwcs.wcs.BasePipeline` subclass as the + ``forward_transform`` argument. In this case, no validation checks or + copies of the steps are made. This means that when the ``forward_transform`` + is a type of `~gwcs.wcs.BasePipeline`, the user is responsible for ensuring + that the pipeline is valid. Moreover, since no copies of the steps are made, + any changes made to the steps contained within the resulting instance will + back-propagate to the original `~gwcs.wcs.BasePipeline` subclass instance. + This behavior is intended to allow for the efficient creation of pipelines + (or WCSs) from existing pipelines (or WCSs) without the expensive overhead + of validation and copies. """ @overload @@ -52,7 +344,7 @@ def __init__( @overload def __init__( self, - forward_transform: Sequence[Step | StepTuple], + forward_transform: Sequence[Step | StepTuple] | BasePipeline, *, input_frame: None = None, output_frame: None = None, @@ -65,8 +357,11 @@ def __init__( input_frame: str | CoordinateFrameProtocol | None = None, output_frame: str | CoordinateFrameProtocol | None = None, ) -> None: - self._pipeline: list[Step] = [] - self._initialize_pipeline(forward_transform, input_frame, output_frame) + if isinstance(forward_transform, BasePipeline): + super().__init__(forward_transform) + else: + super().__init__() + self._initialize_pipeline(forward_transform, input_frame, output_frame) def _initialize_pipeline( self, @@ -175,25 +470,6 @@ def _initialize_pipeline( if isinstance(self._pipeline[-1].frame, EmptyFrame): self._pipeline[-1].frame.naxes = self.forward_transform.n_outputs - @property - def pipeline(self) -> list[Step]: - """ - Allow direct access to the raw pipeline steps. - """ - - # TODO: This can still allow direct modification of the pipeline list - # without any of the checks and handling that have been put in - # place in order to ensure the pipeline is functional. - # -> Maybe we should return a copy? - return self._pipeline - - @property - def available_frames(self) -> list[str]: - """ - List of all the frame names in this WCS in their order in the pipeline - """ - return [step.frame.name for step in self._pipeline] - def get_frame( self, frame: str | CoordinateFrameProtocol ) -> CoordinateFrameProtocol: @@ -225,7 +501,7 @@ def _wrap_step( The step to wrap in a Step object and check. replace_index : int or None The index of the step to replace in the pipeline, this ensures that - we can inplace replace a step using the same frame as the one being + we can in place replace a step using the same frame as the one being replaced. This frame will be removed from the frames to check against If None (default), do not remove any frames for checking. """ @@ -293,41 +569,6 @@ def unit(self) -> Unit | None: """The unit of the coordinates in the output coordinate system.""" return self._pipeline[-1].frame.unit if self._pipeline else None - @staticmethod - def _frame_name(frame: str | CoordinateFrameProtocol) -> str: - """ - Return the name of the frame. - - Parameters - ---------- - frame : str or `~gwcs.coordinate_frames.CoordinateFrameProtocol` - Name of the frame or the frame object. - - Returns - ------- - Name of the frame. - """ - return frame.name if isinstance(frame, CoordinateFrameProtocol) else frame - - def _frame_index(self, frame: str | CoordinateFrameProtocol) -> int: - """ - Return the index of the given frame in the pipeline. - - Parameters - ---------- - frame : str or `~gwcs.coordinate_frames.CoordinateFrameProtocol` - Name of the frame or the frame object. - - Returns - ------- - Index of the frame in the pipeline. - """ - try: - return self.available_frames.index(self._frame_name(frame)) - except ValueError as err: - msg = f"Frame {self._frame_name(frame)} is not in the available frames" - raise CoordinateFrameError(msg) from err - def _get_step(self, frame: str | CoordinateFrameProtocol) -> IndexedStep: """ Get the index and step corresponding to the given frame. @@ -336,123 +577,6 @@ def _get_step(self, frame: str | CoordinateFrameProtocol) -> IndexedStep: return IndexedStep(index, self._pipeline[index]) - @property - def forward_transform(self) -> Model: - """ - Return the forward transform of the pipeline. - """ - transform = combine_transforms([step.transform for step in self._pipeline[:-1]]) - - if self.bounding_box is not None: - # Currently compound models do not attempt to combine individual model - # bounding boxes. Get the forward transform and assign the bounding_box - # to it before evaluating it. The order Model.bounding_box is reversed. - transform.bounding_box = self.bounding_box - - return transform - - @property - def backward_transform(self) -> Model: - """ - Return the total backward transform if available - from output to input - coordinate system. - - Raises - ------ - NotImplementedError : - An analytical inverse does not exist. - - """ - try: - backward = self.forward_transform.inverse - except NotImplementedError: - # Try to get the backward transform by combining the inverses of the - # individual steps. This will work even if the forward transform does - # not have an analytical inverse, as long as each step does. - backward = combine_transforms( - [step.inverse for step in self._pipeline[:-1][::-1]] - ) - - try: - _ = backward.inverse - except NotImplementedError: # means "hasattr" won't work - backward.inverse = self.forward_transform - return backward - - def pipeline_between( - self, - from_frame: str | CoordinateFrameProtocol, - to_frame: str | CoordinateFrameProtocol, - ) -> DirectionalPipeline[Self] | None: - """ - Return a pipeline between the two given frames. - - Parameters - ---------- - from_frame : str or `~gwcs.coordinate_frames.CoordinateFrame` - Initial coordinate frame name of object. - to_frame : str or `~gwcs.coordinate_frames.CoordinateFrame` - End coordinate frame name or object. - - Returns - ------- - pipeline : `~gwcs.wcs.Pipeline` - A NamedTuple (pipeline, forward) where pipeline is a Pipeline object - the pipeline between the two frames and forward is a boolean indicating - if the pipeline is forward (from_frame to to_frame) or - backward (to_frame to from_frame). - """ - from_index = self._frame_index(from_frame) - to_index = self._frame_index(to_frame) - - # Moving backwards over the pipeline - if to_index < from_index: - pipeline = self._pipeline[to_index:from_index] - pipeline.append(Step(self._pipeline[from_index].frame, None)) - - return DirectionalPipeline(wcs=type(self)(pipeline), forward=False) - - if from_index < to_index: - pipeline = self._pipeline[from_index:to_index] - pipeline.append(Step(self._pipeline[to_index].frame, None)) - - return DirectionalPipeline(wcs=type(self)(pipeline), forward=True) - - return None # from and to are the same frame, so no pipeline needed - - def get_transform( - self, - from_frame: str | CoordinateFrameProtocol, - to_frame: str | CoordinateFrameProtocol, - ) -> Mdl: - """ - Return a transform between two coordinate frames. - - Parameters - ---------- - from_frame : str or `~gwcs.coordinate_frames.CoordinateFrameProtocol` - Initial coordinate frame name of object. - to_frame : str or `~gwcs.coordinate_frames.CoordinateFrameProtocol` - End coordinate frame name or object. - - Returns - ------- - transform : `~astropy.modeling.Model` - Transform between two frames. - """ - direction = self.pipeline_between(from_frame, to_frame) - - # from and to are the same frame, so no transform needed - if direction is None: - return None - - # If the pipeline is forward, return the forward transform - if direction.forward: - return direction.wcs.forward_transform - - # Otherwise it is backward, so return the backward transform - return direction.wcs.backward_transform - def set_transform( self, from_frame: str | CoordinateFrameProtocol, @@ -603,9 +727,9 @@ def bounding_box(self) -> ModelBoundingBox | CompoundBoundingBox | None: if transform is None: return None - try: - bounding_box = transform.bounding_box - except NotImplementedError: + bounding_box = super().bounding_box + + if bounding_box is None: return None if ( diff --git a/gwcs/wcs/_step.py b/gwcs/wcs/_step.py index 2cb28790..a8191b3a 100644 --- a/gwcs/wcs/_step.py +++ b/gwcs/wcs/_step.py @@ -1,5 +1,7 @@ import sys import warnings +from copy import copy +from inspect import getattr_static from typing import NamedTuple, Self, TypeAlias, Union from astropy.modeling.core import Model @@ -10,6 +12,7 @@ CoordinateFrameProtocol, EmptyFrame, ) +from gwcs.coordinate_frames._base import _is_high_level, _LegacyCoordinateFrameProtocol __all__ = [ "IndexedStep", @@ -31,11 +34,55 @@ def _is_coordinate_frame(frame: str | CoordinateFrameProtocol) -> bool: return isinstance(frame, CoordinateFrameProtocol) + + def _is_legacy_coordinate_frame( + frame: str | CoordinateFrameProtocol | _LegacyCoordinateFrameProtocol, + ) -> bool: + return isinstance(frame, _LegacyCoordinateFrameProtocol) and not isinstance( + frame, CoordinateFrameProtocol + ) else: def _is_coordinate_frame(frame: str | CoordinateFrameProtocol) -> bool: return isinstance(frame, BaseCoordinateFrame | CoordinateFrame | EmptyFrame) + def _has_legacy_coordinate_frame_interface(frame: object) -> bool: + """ + Return `True` if ``frame`` looks like a legacy coordinate frame object. + + This supports duck-typed frames implementing the historical coordinate + frame API without ``is_high_level``. + """ + + required_members = ( + "naxes", + "name", + "unit", + "axes_names", + "axes_order", + "reference_frame", + "axes_type", + "axis_physical_types", + "world_axis_object_classes", + "world_axis_object_components", + "add_units", + "remove_units", + "to_high_level_coordinates", + "from_high_level_coordinates", + ) + + return all( + getattr_static(frame, member, None) is not None + for member in required_members + ) + + def _is_legacy_coordinate_frame( + frame: str | CoordinateFrameProtocol, + ) -> bool: + return _has_legacy_coordinate_frame_interface(frame) and not hasattr( + frame, "is_high_level" + ) + class Step: """ @@ -58,7 +105,7 @@ def __init__( # This is correct type-wise, but the Python 3.11 bugfix causes a MyPy error self.frame = ( frame - if _is_coordinate_frame(frame) + if _is_coordinate_frame(frame) or _is_legacy_coordinate_frame(frame) else EmptyFrame.from_transform(frame, transform) # type: ignore[assignment, arg-type] ) self.transform = transform @@ -69,7 +116,18 @@ def frame(self) -> CoordinateFrameProtocol: @frame.setter def frame(self, val: CoordinateFrameProtocol) -> None: - if not _is_coordinate_frame(val): + if is_legacy := _is_legacy_coordinate_frame(val): + msg = ( + "Coordinate frames that do not implement `is_high_level` are " + "deprecated. Please update your coordinate frame to add " + "`is_high_level`." + ) + warnings.warn(msg, DeprecationWarning, stacklevel=2) + # Copy the value to avoid mutating the original object. + val = copy(val) + val.is_high_level = lambda *args: _is_high_level(val, *args) # type: ignore[method-assign] + + if not (_is_coordinate_frame(val) or is_legacy): msg = '"frame" should be an instance of CoordinateFrameProtocol.' raise TypeError(msg) diff --git a/gwcs/wcs/_wcs.py b/gwcs/wcs/_wcs.py index 6aa5186e..4c2243b8 100644 --- a/gwcs/wcs/_wcs.py +++ b/gwcs/wcs/_wcs.py @@ -2,6 +2,7 @@ from __future__ import annotations import functools +import inspect import itertools import sys import warnings @@ -37,10 +38,10 @@ LowLevelInput, get_ctype_from_ucd, ) -from gwcs.utils import _compute_lon_pole, to_index +from gwcs.utils import _compute_lon_pole, correct_1d_output, to_index from ._exception import NoConvergence -from ._pipeline import ForwardTransform, Pipeline +from ._pipeline import BasePipeline, ForwardTransform, Pipeline from ._step import Step, StepTuple from ._utils import ( fit_2D_poly, @@ -89,12 +90,75 @@ def __init__(self, axis, frame, world_axis_order, cunit, ctype, input_axes): self.input_axes = input_axes +@functools.cache +def _func_accepts_correct_1d(func) -> bool: + """ + Determine whether a frame's coordinate-conversion method accepts the + ``correct_1d`` keyword. Legacy frames predate this argument. + """ + try: + sig = inspect.signature(func) + except (TypeError, ValueError): + return False + + return "correct_1d" in sig.parameters or any( + p.kind is inspect.Parameter.VAR_KEYWORD for p in sig.parameters.values() + ) + + +def _accepts_correct_1d(method) -> bool: + """ + Helper to pass to the _func_accepts_correct_1d function so that the cache + is not grown every time this function is called. + -> This is because the method is a bound method. + """ + return _func_accepts_correct_1d(getattr(method, "__func__", method)) + + +def _from_high_level_coordinates(frame, *args, correct_1d=True): + """ + Helper to support legacy frames whose ``from_high_level_coordinates`` does + not implement the ``correct_1d`` argument. + """ + method = frame.from_high_level_coordinates + if _accepts_correct_1d(method): + return method(*args, correct_1d=correct_1d) + + return correct_1d_output(method, pass_correct_1d=False)( + *args, correct_1d=correct_1d + ) + + +def _to_high_level_coordinates(frame, *args, correct_1d=True): + """ + Helper to support legacy frames whose ``to_high_level_coordinates`` does + not implement the ``correct_1d`` argument. + """ + method = frame.to_high_level_coordinates + if _accepts_correct_1d(method): + return method(*args, correct_1d=correct_1d) + + return correct_1d_output(method, pass_correct_1d=False)( + *args, correct_1d=correct_1d + ) + + class _UnitHandler: """ Class to ensure that the input-output type consistency for the GWCS native API Parameters ---------- + inputs : tuple of low level inputs + The inputs to be passed to the transform. + transform : Model + The transform to be applied to the inputs. + frame : CoordinateFrameProtocol + The frame to use for adding/removing units and converting to high level + objects. + + Attributes + ---------- args : The input arguments to be passed to the transform. _is_high_level : bool @@ -112,14 +176,21 @@ class _UnitHandler: def __init__( self, - inputs: tuple[LowLevelInput, ...], + inputs: tuple[LowLevelInput, ...] | LowLevelInput, transform: Model, frame: CoordinateFrameProtocol, ): # Handle the case where a HLO is passed into the Native API self._is_high_level = frame.is_high_level(*inputs) if self._is_high_level: - inputs = frame.from_high_level_coordinates(*inputs, correct_1d=False) + inputs = _from_high_level_coordinates(frame, *inputs, correct_1d=False) + + # Legacy frames will probably have stripped the 1d to a value rather than + # keeping it as a tuple + if not isinstance(inputs, list | tuple): + inputs = (inputs,) + + inputs = tuple(inputs) # Determine if the inputs are quantities input_is_quantity = any(isinstance(a, u.Quantity) for a in inputs) @@ -133,20 +204,28 @@ def __init__( or transform.input_units_equivalencies is not None ) - elif isinstance(transform, Identity): + # Note that upstream changes in astropy are planned to make parameterless + # models that have no parameters be able to set uses_quantity on their + # own. It is an ongoing discussion in astropy about whether Identity + # should be considered to use quantities or not (as there are other + # models like Scale and Shift that can be setup to effectively be an + # Identity and they make it clear when they use quantities or not). + # For now, we will assume that Identity does not use quantities, as that + # is the most common expected behavior of downstream users of GWCS. + elif isinstance(transform, Identity) or ( + transform is not None and len(transform.parameters) == 0 + ): transform_uses_quantity = False else: - transform_uses_quantity = not ( - transform is None or not transform.uses_quantity - ) + transform_uses_quantity = transform is not None and transform.uses_quantity # Using input_is_quantity and transform_uses_quantity, we determine if we # need to add or remove units from the inputs to make them compatible # with the transform and we determine if we need to add or remove units # from the outputs to make them consistent with the input's type. - # Case 1: No units are involved + # Case 1: No units are involved in anything # input -> no change # (we assume the user has passed numerical values # in the correct units) @@ -157,16 +236,16 @@ def __init__( self.args = inputs self._add_units = False - # Case 2: Units are involved and the transform supports them + # Case 2: The inputs have units and the transform supports units # input -> add units # (we add -- meaning convert -- units to the input values) # output -> add units - # (we convert the output to the correct units if necessary and + # (we add -- meaning convert -- units to the output values) elif input_is_quantity and transform_uses_quantity: self.args = frame.add_units(inputs) self._add_units = True - # Case 3: no input units, but transform needs units + # Case 3: The input does not have units, but the transform supports units # input -> add units # (we assume that the user has passed correct numerical values # and formally add units to them for the transform to work @@ -178,13 +257,13 @@ def __init__( self.args = frame.add_units(inputs) self._add_units = False - # Case 4: input has units, but transform does not use them + # Case 4: The input has units, but the transform does not support units # input -> remove units # (we convert the input to the correct units if necessary and # then remove them giving numerical values as input) # output -> add units - # (we convert the output to the correct units if necessary and - # then add them giving numerical values as output) + # (we assume that the output of the transform is in the correct + # units and then formally attach units to the output values) # This is the only other case: # input_is_quantity and not transform_uses_quantity else: @@ -196,7 +275,12 @@ def _handle_output_type( ) -> tuple[LowLevelInput, ...]: # Return a high level object if the input was a high level object if self._is_high_level: - return frame.to_high_level_coordinates(*outputs, correct_1d=False) # type: ignore[no-any-return] + values = _to_high_level_coordinates(frame, *outputs, correct_1d=False) # type: ignore[no-any-return] + + if not isinstance(values, list | tuple): + values = (values,) + + return tuple(values) # If the input had units return units if self._add_units: @@ -271,7 +355,7 @@ def __init__( @overload def __init__( self, - forward_transform: list[Step | StepTuple], + forward_transform: list[Step | StepTuple] | BasePipeline, input_frame: None = None, output_frame: None = None, name: str | None = None, @@ -471,8 +555,8 @@ def outside_footprint(self, world_arrays): not_numerical = False if self.output_frame.is_high_level(*world_arrays): not_numerical = True - world_arrays = self.output_frame.from_high_level_coordinates( - *world_arrays, correct_1d=False + world_arrays = _from_high_level_coordinates( + self.output_frame, *world_arrays, correct_1d=False ) for axtyp in axes_types: @@ -518,8 +602,8 @@ def outside_footprint(self, world_arrays): coord[outside] = np.nan world_arrays[idim] = coord if not_numerical: - world_arrays = self.output_frame.to_high_level_coordinates( - *world_arrays, correct_1d=False + world_arrays = _to_high_level_coordinates( + self.output_frame, *world_arrays, correct_1d=False ) return world_arrays diff --git a/gwcs/wcstools.py b/gwcs/wcstools.py index 1b11edd7..fbc7900f 100644 --- a/gwcs/wcstools.py +++ b/gwcs/wcstools.py @@ -64,7 +64,8 @@ def wcs_from_fiducial( The bounding box over which the WCS is valid. It is a tuple of tuples of size 2 where each tuple represents a range of (low, high) values. The ``bounding_box`` is in the - order of the axes, `~gwcs.coordinate_frames.CoordinateFrameProtocol.axes_order`. + order of the axes, + `~gwcs.coordinate_frames._LegacyCoordinateFrameProtocol.axes_order`. For two inputs and axes_order(0, 1) the bounding box is ((xlow, xhigh), (ylow, yhigh)). input_frame : ~gwcs.coordinate_frames.CoordinateFrameProtocol`