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"), 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/__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 284ff68f..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,7 +290,8 @@ def remove_units( for array in self.add_units(arrays) ) - def to_high_level_coordinates(self, *values): + @correct_1d_output + def to_high_level_coordinates(self, *values, correct_1d=True): """ Convert "values" to high level coordinate objects described by this frame. @@ -308,12 +317,10 @@ def to_high_level_coordinates(self, *values): 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 len(high_level) == 1: - high_level = high_level[0] - return high_level + return values_to_high_level_objects(*values, low_level_wcs=self) - def from_high_level_coordinates(self, *high_level_coords): + @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. @@ -331,10 +338,76 @@ def from_high_level_coordinates(self, *high_level_coords): 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 len(values) == 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 e46f47a7..614102d7 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,12 @@ def test_coordinate_frame_api(): pixel = wcs.world_to_pixel(world) assert isinstance(pixel, float) - with pytest.raises(TypeError): - _ = wcs.invert(world) + # 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_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..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, @@ -646,14 +647,53 @@ def add_units(self, arrays): def remove_units(self, arrays): return arrays - def to_high_level_coordinates(self, *inputs): + def is_high_level(self, *inputs): + return False + + 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 21e37152..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) @@ -1788,10 +1852,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 +2003,298 @@ 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) + 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) - # 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) + + # 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): diff --git a/gwcs/utils.py b/gwcs/utils.py index 40a55f2d..345aa4ba 100644 --- a/gwcs/utils.py +++ b/gwcs/utils.py @@ -495,34 +495,35 @@ def create_projection_transform(projcode): return projklass(**projparams) -def is_high_level(*args, low_level_wcs): +def correct_1d_output(func, pass_correct_1d: bool = True): """ - Determine if args matches the high level classes as defined by - ``low_level_wcs``. + Decorator to correct the output of a function to be 1D if the input is 1D. """ - 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] + @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 all(types_are_high_level): - return True + if correct_1d and isinstance(value, (tuple, list)) and len(value) == 1: + return value[0] - 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 value + + return wrapper - return False + +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 60045198..782e1a5c 100644 --- a/gwcs/wcs/__init__.py +++ b/gwcs/wcs/__init__.py @@ -1,10 +1,14 @@ from ._exception import GwcsBoundingBoxWarning, NoConvergence +from ._pipeline import BasePipeline, DirectionalWCS, Pipeline from ._step import Step from ._wcs import WCS __all__ = [ "WCS", + "BasePipeline", + "DirectionalWCS", "GwcsBoundingBoxWarning", "NoConvergence", + "Pipeline", "Step", ] diff --git a/gwcs/wcs/_pipeline.py b/gwcs/wcs/_pipeline.py index 3020dfb5..17eaf83e 100644 --- a/gwcs/wcs/_pipeline.py +++ b/gwcs/wcs/_pipeline.py @@ -1,7 +1,9 @@ from __future__ import annotations import warnings -from typing import TypeAlias, Union, overload +from collections.abc import Sequence +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 @@ -13,19 +15,321 @@ from ._exception import GwcsBoundingBoxWarning, GwcsFrameExistsError from ._step import IndexedStep, Mdl, Step, StepTuple -__all__ = ["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, list[Step | StepTuple]] # noqa: UP007 +ForwardTransform: TypeAlias = Union[Model, Sequence[Step | StepTuple] | BasePipeline] # noqa: UP007 + + +_T = TypeVar("_T", bound=BasePipeline) + + +@dataclass(frozen=True, slots=True) +class DirectionalWCS(Generic[_T]): + """ + 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 @@ -40,7 +344,7 @@ def __init__( @overload def __init__( self, - forward_transform: list[Step | StepTuple], + forward_transform: Sequence[Step | StepTuple] | BasePipeline, *, input_frame: None = None, output_frame: None = None, @@ -53,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, @@ -103,7 +410,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 +444,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 " @@ -158,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: @@ -208,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. """ @@ -276,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. @@ -319,47 +577,6 @@ def _get_step(self, frame: str | CoordinateFrameProtocol) -> IndexedStep: return IndexedStep(index, self._pipeline[index]) - 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. - """ - from_index = self._frame_index(from_frame) - to_index = self._frame_index(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: - return None - - # Moving forwards over the pipeline - else: - transforms = [ - step.transform for step in self._pipeline[from_index:to_index] - ] - - return combine_transforms(transforms) - def set_transform( self, from_frame: str | CoordinateFrameProtocol, @@ -505,15 +722,14 @@ 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 - try: - bounding_box = transform.bounding_box - except NotImplementedError: + bounding_box = super().bounding_box + + if bounding_box is None: return None if ( @@ -556,8 +772,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 +792,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 +807,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/_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 d6006edf..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 @@ -14,6 +15,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, @@ -23,10 +25,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,10 +38,10 @@ LowLevelInput, get_ctype_from_ucd, ) -from gwcs.utils import _compute_lon_pole, is_high_level, 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, @@ -92,6 +90,238 @@ 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 + 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, ...] | 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 = _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) + + # 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 + ) + + # 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 = 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 in anything + # 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: 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 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: 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 + # 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: 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 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: + 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: + 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: + 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. @@ -125,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, @@ -160,34 +390,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, @@ -214,93 +426,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 or transform.parameters.size - ): - 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, @@ -394,40 +519,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) - - input_is_quantity, transform_uses_quantity = self._units_are_present( - args, transform - ) + inputs = self.outside_footprint(inputs) - 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 +544,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] @@ -456,11 +553,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 = _from_high_level_coordinates( + self.output_frame, *world_arrays, correct_1d=False ) + for axtyp in axes_types: ind = np.asarray(np.asarray(self.output_frame.axes_type) == axtyp) @@ -504,8 +602,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 = _to_high_level_coordinates( + self.output_frame, *world_arrays, correct_1d=False ) return world_arrays @@ -1193,47 +1291,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: 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`