Coverage for python/lsst/images/_transforms/_transform.py: 77%
278 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-16 15:24 -0700
1# This file is part of lsst-images.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# Use of this source code is governed by a 3-clause BSD-style
10# license that can be found in the LICENSE file.
12from __future__ import annotations
14__all__ = (
15 "Transform",
16 "TransformCompositionError",
17 "TransformSerializationModel",
18)
20import textwrap
21from collections.abc import Iterable
22from typing import TYPE_CHECKING, Any, ClassVar, TypeVar, assert_type, cast, final, overload
24import astropy.io.fits.header
25import astropy.units as u
26import numpy as np
27import numpy.typing as npt
28import pydantic
30from .._concrete_bounds import BoundsSerializationModel
31from .._geom import XY, YX, Bounds, Box
32from ..serialization import ArchiveReadError, ArchiveTree, InputArchive, InvalidParameterError, OutputArchive
33from . import _ast as astshim
34from ._frames import Frame, SerializableFrame, SkyFrame
36if TYPE_CHECKING:
37 try:
38 from lsst.afw.geom import TransformPoint2ToPoint2 as LegacyTransform
39 except ImportError:
40 type LegacyTransform = Any # type: ignore[no-redef]
42# These pre-python-3.12 declaration are needed by Sphinx (probably the
43# autodoc-typehints plugin.
44I = TypeVar("I", bound=Frame) # noqa: E741
45O = TypeVar("O", bound=Frame) # noqa: E741
46P = TypeVar("P", bound=pydantic.BaseModel)
49class TransformCompositionError(RuntimeError):
50 """Exception raised when two transforms cannot be composed."""
53@final
54class Transform[I: Frame, O: Frame]:
55 """A transform that maps two coordinate frames.
57 Parameters
58 ----------
59 in_frame
60 Input coordinate frame.
61 out_frame
62 Output coordinate frame.
63 ast_mapping
64 AST mapping that implements the transform.
65 in_bounds
66 Bounds of the input frame, defaulting to the input frame's
67 bounding box.
68 out_bounds
69 Bounds of the output frame, defaulting to the output frame's
70 bounding box.
71 components
72 Component transforms that this transform was composed from.
74 Notes
75 -----
76 The `Transform` class constructor is considered a private implementation
77 detail. Instead of using this, various factory methods are available:
79 - `from_fits_wcs` constructs a transform from a FITS WCS, as represented
80 `astropy.wcs.WCS`;
81 - `then` composes two transforms;
82 - `identity` constructs a trivial transform that does nothing;
83 - `affine` contructs an affine transform from a 2x2 or 3x3 matrix;
84 - `inverted` returns the inverse of a transform;
85 - `from_legacy` converts an `lsst.afw.geom.Transform` instance.
87 When applied to celestial coordinate systems, ``x=ra`` and ``y=dec``.
88 `SkyProjection` provides a more natural interface for pixel-to-sky
89 transforms.
91 `Transform` is conceptually immutable (the internal AST Mapping should
92 never be modified in-place after construction), and hence does not need to
93 be copied when any object that holds it is copied.
94 """
96 def __init__(
97 self,
98 in_frame: I,
99 out_frame: O,
100 ast_mapping: astshim.Mapping,
101 in_bounds: Bounds | None = None,
102 out_bounds: Bounds | None = None,
103 components: Iterable[Transform[Any, Any]] = (),
104 ) -> None:
105 self._in_frame = in_frame
106 self._out_frame = out_frame
107 self._ast_mapping = ast_mapping
108 self._in_bounds = in_bounds or getattr(in_frame, "bbox", None)
109 self._out_bounds = out_bounds or getattr(out_frame, "bbox", None)
110 self._components = list(components)
112 def __eq__(self, other: Any) -> bool:
113 if self is other:
114 # Short circuit for case where you are quickly checking
115 # that the image WCS and variance WCS are the same object.
116 return True
117 if not isinstance(other, Transform):
118 return NotImplemented
119 if self._ast_mapping != other._ast_mapping:
120 return False
121 if self._in_bounds != other._in_bounds:
122 return False
123 if self._out_bounds != other._out_bounds:
124 return False
125 if self._in_frame != other._in_frame:
126 return False
127 if self._out_frame != other._out_frame:
128 return False
129 if self._components != other._components:
130 return False
131 return True
133 @staticmethod
134 def from_fits_wcs(
135 fits_wcs: astropy.wcs.WCS,
136 in_frame: I,
137 out_frame: O,
138 in_bounds: Bounds | None = None,
139 out_bounds: Bounds | None = None,
140 x0: int = 0,
141 y0: int = 0,
142 ) -> Transform[I, O]:
143 """Construct a transform from a FITS WCS.
145 Parameters
146 ----------
147 fits_wcs
148 FITS WCS to convert.
149 in_frame
150 Coordinate frame for input points to the forward transform.
151 out_frame
152 Coordinate frame for output points from the forward transform.
153 in_bounds
154 The region that bounds valid input points.
155 out_bounds
156 The region that bounds valid output points.
157 x0
158 Logical coordinate of the first column in the array this WCS
159 relates to world coordinates.
160 y0
161 Logical coordinate of the first column in the array this WCS
162 relates to world coordinates.
164 Notes
165 -----
166 The ``x0`` and ``y0`` parameters reflect the fact that for FITS, the
167 first row and column are always labeled ``(1, 1)``, while in Astropy
168 and most other Python libraries, they are ``(0, 0)``. The `types` in
169 this package (e.g. `Image`, `Mask`) allow them to be any pair of
170 integers.
172 See Also
173 --------
174 SkyProjection.from_fits_wcs
175 """
176 ast_stream = astshim.StringStream(fits_wcs.to_header_string(relax=True))
177 ast_fits_chan = astshim.FitsChan(ast_stream, "Encoding=FITS-WCS, SipReplace=0, IWC=1")
178 ast_frame_set = ast_fits_chan.read()
179 _prepend_ast_shift(ast_frame_set, x=x0 - 1.0, y=y0 - 1.0, ast_domain="PIXEL")
180 return Transform(
181 in_frame,
182 out_frame,
183 ast_frame_set,
184 in_bounds=in_bounds,
185 out_bounds=out_bounds,
186 )
188 @staticmethod
189 def identity(frame: I) -> Transform[I, I]:
190 """Construct a trivial transform that maps a frame to itelf.
192 Parameters
193 ----------
194 frame
195 Frame used for both input and output points.
196 """
197 return Transform(frame, frame, astshim.UnitMap(2))
199 @staticmethod
200 def affine(in_frame: I, out_frame: O, matrix: np.ndarray) -> Transform[I, O]:
201 """Construct an affine transform from a matrix.
203 Parameters
204 ----------
205 in_frame
206 Coordinate frame for input points to the forward transform.
207 out_frame
208 Coordinate frame for output points from the forward transform.
209 matrix
210 Matrix of coefficients, either a 2x2 linear transform or a 3x3
211 augmented affine transform, with a shift embedded in the third
212 column and ``[0, 0, 1]`` the third row.
213 """
214 if matrix.shape == (2, 2):
215 return Transform(in_frame, out_frame, astshim.MatrixMap(matrix.copy()))
216 elif matrix.shape == (3, 3): 216 ↛ 223line 216 didn't jump to line 223 because the condition on line 216 was always true
217 linear = astshim.MatrixMap(matrix[:2, :2].copy())
218 shift = astshim.ShiftMap(matrix[:2, 2])
219 if not np.array_equal(matrix[2, :], np.array([0.0, 0.0, 1.0])): 219 ↛ 220line 219 didn't jump to line 220 because the condition on line 219 was never true
220 raise ValueError("3x3 affine transform array must have [0, 0, 1] in its last row.")
221 return Transform(in_frame, out_frame, linear.then(shift))
222 else:
223 raise ValueError("Affine transform array must be 2x2 or 3x3.")
225 @property
226 def in_frame(self) -> I:
227 """Coordinate frame for input points."""
228 return self._in_frame
230 @property
231 def out_frame(self) -> O:
232 """Coordinate frame for output points."""
233 return self._out_frame
235 @property
236 def in_bounds(self) -> Bounds | None:
237 """The region that bounds valid input points (`Bounds` | `None`)."""
238 return self._in_bounds
240 @property
241 def out_bounds(self) -> Bounds | None:
242 """The region that bounds valid output points (`Bounds` | `None`)."""
243 return self._out_bounds
245 def show(self, simplified: bool = False, comments: bool = False) -> str:
246 """Return the AST native representation of the transform.
248 Parameters
249 ----------
250 simplified
251 Whether to ask AST to simplify the mapping before showing it.
252 This will make it much more likely that two equivalent transforms
253 have the same `show` result. If the internal mapping is actually
254 a frame set (as needed to round-trip legacy
255 `lsst.afw.geom.SkyWcs` objects), this will also just show the
256 mapping with no frame set information.
257 comments
258 Whether to include descriptive comments.
259 """
260 ast_mapping = self._ast_mapping
261 if simplified: 261 ↛ 262line 261 didn't jump to line 262 because the condition on line 261 was never true
262 if isinstance(ast_mapping, astshim.FrameSet):
263 ast_mapping = ast_mapping.getMapping()
264 ast_mapping = ast_mapping.simplified()
265 return ast_mapping.show(comments)
267 @overload
268 def apply_forward(self, point: XY[int | float] | YX[int | float], /) -> XY[float]: ... 268 ↛ exitline 268 didn't return from function 'apply_forward' because
270 @overload
271 def apply_forward(self, point: XY[npt.ArrayLike] | YX[npt.ArrayLike], /) -> XY[np.ndarray]: ... 271 ↛ exitline 271 didn't return from function 'apply_forward' because
273 @overload
274 def apply_forward(self, /, *, x: int | float, y: int | float) -> XY[float]: ... 274 ↛ exitline 274 didn't return from function 'apply_forward' because
276 @overload
277 def apply_forward(self, /, *, x: npt.ArrayLike, y: npt.ArrayLike) -> XY[np.ndarray]: ... 277 ↛ exitline 277 didn't return from function 'apply_forward' because
279 def apply_forward(
280 self, point: XY[Any] | YX[Any] | None = None, /, *, x: Any = None, y: Any = None
281 ) -> XY[float] | XY[np.ndarray]:
282 """Apply the forward transform to one or more points.
284 Parameters
285 ----------
286 point
287 An `XY` or `YX` coordinate pair to transform. Mutually exclusive
288 with ``x`` and ``y``.
289 x : `float` | array-like
290 ``x`` values of the points to transform, as a scalar or any
291 array-like. Results are broadcast against ``y``.
292 Mutually exclusive with ``point``.
293 y : `float` | array-like
294 ``y`` values of the points to transform, as a scalar or any
295 array-like. Results are broadcast against ``x``.
296 Mutually exclusive with ``point``.
298 Returns
299 -------
300 `XY` [`float` | `numpy.ndarray`]
301 The transformed point or points. A scalar input pair returns
302 `XY` of `float`; array-like inputs return `XY` of
303 `numpy.ndarray` with the broadcast shape of ``x`` and ``y``.
304 """
305 match point:
306 case None:
307 if x is None or y is None:
308 raise TypeError("Pass either a point or both x= and y= to 'apply_forward'.")
309 case XY() | YX(): 309 ↛ 313line 309 didn't jump to line 313 because the pattern on line 309 always matched
310 if x is not None or y is not None:
311 raise TypeError("'apply_forward' point argument is mutually exclusive with x= and y=.")
312 x, y = point.x, point.y
313 case _:
314 raise TypeError(f"Unexpected positional argument type: {type(point)!r}.")
315 return _standardize_xy(
316 _ast_apply(
317 self._ast_mapping.applyForward,
318 x=self._in_frame.standardize_x(x),
319 y=self._in_frame.standardize_y(y),
320 ),
321 self._out_frame,
322 )
324 @overload
325 def apply_inverse(self, point: XY[int | float] | YX[int | float], /) -> XY[float]: ... 325 ↛ exitline 325 didn't return from function 'apply_inverse' because
327 @overload
328 def apply_inverse(self, point: XY[npt.ArrayLike] | YX[npt.ArrayLike], /) -> XY[np.ndarray]: ... 328 ↛ exitline 328 didn't return from function 'apply_inverse' because
330 @overload
331 def apply_inverse(self, /, *, x: int | float, y: int | float) -> XY[float]: ... 331 ↛ exitline 331 didn't return from function 'apply_inverse' because
333 @overload
334 def apply_inverse(self, /, *, x: npt.ArrayLike, y: npt.ArrayLike) -> XY[np.ndarray]: ... 334 ↛ exitline 334 didn't return from function 'apply_inverse' because
336 def apply_inverse(
337 self, point: XY[Any] | YX[Any] | None = None, /, *, x: Any = None, y: Any = None
338 ) -> XY[float] | XY[np.ndarray]:
339 """Apply the inverse transform to one or more points.
341 Parameters
342 ----------
343 point
344 An `XY` or `YX` coordinate pair to transform. Mutually exclusive
345 with ``x`` and ``y``.
346 x : `float` | array-like
347 ``x`` values of the points to transform, as a scalar or any
348 array-like. Results are broadcast against ``y``.
349 Mutually exclusive with ``point``.
350 y : `float` | array-like
351 ``y`` values of the points to transform, as a scalar or any
352 array-like. Results are broadcast against ``x``.
353 Mutually exclusive with ``point``.
355 Returns
356 -------
357 `XY` [`float` | `numpy.ndarray`]
358 The transformed point or points. A scalar input pair returns
359 `XY` of `float`; array-like inputs return `XY` of
360 `numpy.ndarray` with the broadcast shape of ``x`` and ``y``.
361 """
362 match point:
363 case None:
364 if x is None or y is None: 364 ↛ 365line 364 didn't jump to line 365 because the condition on line 364 was never true
365 raise TypeError("Pass either a point or both x= and y= to 'apply_inverse'.")
366 case XY() | YX(): 366 ↛ 370line 366 didn't jump to line 370 because the pattern on line 366 always matched
367 if x is not None or y is not None: 367 ↛ 368line 367 didn't jump to line 368 because the condition on line 367 was never true
368 raise TypeError("'apply_inverse' point argument is mutually exclusive with x= and y=.")
369 x, y = point.x, point.y
370 case _:
371 raise TypeError(f"Unexpected positional argument type: {type(point)!r}.")
372 return _standardize_xy(
373 _ast_apply(
374 self._ast_mapping.applyInverse,
375 x=self._out_frame.standardize_x(x),
376 y=self._out_frame.standardize_y(y),
377 ),
378 self._in_frame,
379 )
381 @overload
382 def apply_forward_q(self, point: XY[u.Quantity] | YX[u.Quantity], /) -> XY[u.Quantity]: ... 382 ↛ exitline 382 didn't return from function 'apply_forward_q' because
384 @overload
385 def apply_forward_q(self, /, *, x: u.Quantity, y: u.Quantity) -> XY[u.Quantity]: ... 385 ↛ exitline 385 didn't return from function 'apply_forward_q' because
387 def apply_forward_q(
388 self, point: XY[u.Quantity] | YX[u.Quantity] | None = None, /, *, x: Any = None, y: Any = None
389 ) -> XY[u.Quantity]:
390 """Apply the forward transform to one or more unit-aware points.
392 Parameters
393 ----------
394 point
395 An `XY` or `YX` coordinate pair of `~astropy.units.Quantity` to
396 transform. Mutually exclusive with ``x`` and ``y``.
397 x
398 ``x`` values of the points to transform.
399 Mutually exclusive with ``point``.
400 y
401 ``y`` values of the points to transform.
402 Mutually exclusive with ``point``.
404 Returns
405 -------
406 `XY` [`astropy.units.Quantity`]
407 The transformed point or points.
408 """
409 match point:
410 case None:
411 if x is None or y is None: 411 ↛ 412line 411 didn't jump to line 412 because the condition on line 411 was never true
412 raise TypeError("Pass either a point or both x= and y= to 'apply_forward_q'.")
413 case XY() | YX(): 413 ↛ 417line 413 didn't jump to line 417 because the pattern on line 413 always matched
414 if x is not None or y is not None: 414 ↛ 415line 414 didn't jump to line 415 because the condition on line 414 was never true
415 raise TypeError("'apply_forward_q' point argument is mutually exclusive with x= and y=.")
416 x, y = point.x, point.y
417 case _:
418 raise TypeError(f"Unexpected positional argument type: {type(point)!r}.")
419 xy = self.apply_forward(x=x.to_value(self._in_frame.unit), y=y.to_value(self._in_frame.unit))
420 return XY(xy.x * self._out_frame.unit, xy.y * self._out_frame.unit)
422 @overload
423 def apply_inverse_q(self, point: XY[u.Quantity] | YX[u.Quantity], /) -> XY[u.Quantity]: ... 423 ↛ exitline 423 didn't return from function 'apply_inverse_q' because
425 @overload
426 def apply_inverse_q(self, /, *, x: u.Quantity, y: u.Quantity) -> XY[u.Quantity]: ... 426 ↛ exitline 426 didn't return from function 'apply_inverse_q' because
428 def apply_inverse_q(
429 self, point: XY[u.Quantity] | YX[u.Quantity] | None = None, /, *, x: Any = None, y: Any = None
430 ) -> XY[u.Quantity]:
431 """Apply the inverse transform to one or more unit-aware points.
433 Parameters
434 ----------
435 point
436 An `XY` or `YX` coordinate pair of `~astropy.units.Quantity` to
437 transform. Mutually exclusive with ``x`` and ``y``.
438 x
439 ``x`` values of the points to transform.
440 Mutually exclusive with ``point``.
441 y
442 ``y`` values of the points to transform.
443 Mutually exclusive with ``point``.
445 Returns
446 -------
447 `XY` [`astropy.units.Quantity`]
448 The transformed point or points.
449 """
450 match point:
451 case None:
452 if x is None or y is None: 452 ↛ 453line 452 didn't jump to line 453 because the condition on line 452 was never true
453 raise TypeError("Pass either a point or both x= and y= to 'apply_inverse_q'.")
454 case XY() | YX(): 454 ↛ 458line 454 didn't jump to line 458 because the pattern on line 454 always matched
455 if x is not None or y is not None: 455 ↛ 456line 455 didn't jump to line 456 because the condition on line 455 was never true
456 raise TypeError("'apply_inverse_q' point argument is mutually exclusive with x= and y=.")
457 x, y = point.x, point.y
458 case _:
459 raise TypeError(f"Unexpected positional argument type: {type(point)!r}.")
460 xy = self.apply_inverse(x=x.to_value(self._out_frame.unit), y=y.to_value(self._out_frame.unit))
461 return XY(xy.x * self._in_frame.unit, xy.y * self._in_frame.unit)
463 def decompose(self) -> list[Transform[Any, Any]]:
464 """Deconstruct a composed transform into its constituent parts.
466 Notes
467 -----
468 Most transforms will just return a single-element list holding
469 ``self``. Identity transform will return an empty list, and
470 transforms composed with `then` will return the original transforms.
471 Transforms constructed by `FrameSet` may or may not be decomposable.
472 """
473 if not self._components: 473 ↛ 479line 473 didn't jump to line 479 because the condition on line 473 was always true
474 if self.in_frame == self._out_frame: 474 ↛ 477line 474 didn't jump to line 477 because the condition on line 474 was always true
475 return []
476 else:
477 return [self]
478 else:
479 return list(self._components)
481 def inverted(self) -> Transform[O, I]:
482 """Return the inverse of this transform."""
483 return Transform[O, I](
484 self._out_frame,
485 self._in_frame,
486 self._ast_mapping.inverted(),
487 in_bounds=self.out_bounds,
488 out_bounds=self.in_bounds,
489 components=[t.inverted() for t in reversed(self._components)],
490 )
492 def then[F: Frame](self, next: Transform[O, F], remember_components: bool = True) -> Transform[I, F]:
493 """Compose two transforms into another.
495 Parameters
496 ----------
497 next
498 Another transform to apply after ``self``.
499 remember_components
500 If `True`, the returned composed transform will remember ``self``
501 and ``other`` so they can be returned by `decompose`.
502 """
503 if self._out_frame != next._in_frame:
504 raise TransformCompositionError(
505 "Cannot compose transforms that do not share a common intermediate frame: "
506 f"{self._out_frame} != {next._in_frame}."
507 )
508 components = self.decompose() + next.decompose() if remember_components else ()
509 return Transform(
510 self._in_frame,
511 next._out_frame,
512 self._ast_mapping.then(next._ast_mapping),
513 in_bounds=self.in_bounds,
514 out_bounds=next.out_bounds,
515 components=components,
516 )
518 def as_fits_wcs(self, bbox: Box) -> astropy.wcs.WCS | None:
519 """Return a FITS WCS representation of this transform, if possible.
521 Parameters
522 ----------
523 bbox
524 Bounding box of the array the FITS WCS will describe. This
525 transform object is assumed to work on the same coordinate system
526 in which ``bbox`` is defined, while the FITS WCS will consider the
527 first row and column in that box to be ``(0, 0)`` (in Astropy
528 interfaces) or ``(1, 1)`` (in the FITS representation itself).
530 Notes
531 -----
532 This method assumes the transform maps pixel coordinates to world
533 coordinates.
535 Not all transforms can be represented exactly; when a FITS
536 represention is not possible, `None` is returned. When the returned
537 WCS is not `None`, it will have the same functional form, but it may
538 not evaluate identically due to small implementation differences in
539 the order of floating-point operations.
540 """
541 ast_frame_set = self._get_ast_frame_set()
542 _prepend_ast_shift(ast_frame_set, x=1.0 - bbox.x.start, y=1.0 - bbox.y.start, ast_domain="GRID")
543 ast_stream = astshim.StringStream()
544 ast_fits_chan = astshim.FitsChan(
545 ast_stream, "Encoding=FITS-WCS, CDMatrix=1, FitsAxisOrder=<copy>, FitsTol=0.0001"
546 )
547 ast_fits_chan.setFitsI("NAXIS1", bbox.x.size)
548 ast_fits_chan.setFitsI("NAXIS2", bbox.y.size)
549 n_writes = ast_fits_chan.write(ast_frame_set)
550 if not n_writes:
551 return None
552 header = astropy.io.fits.Header(astropy.io.fits.Card.fromstring(c) for c in ast_fits_chan)
553 return astropy.wcs.WCS(header)
555 def serialize[P: pydantic.BaseModel](
556 self, archive: OutputArchive[P], *, use_frame_sets: bool = False
557 ) -> TransformSerializationModel[P]:
558 """Serialize a transform to an archive.
560 Parameters
561 ----------
562 archive
563 Archive to serialize to.
564 use_frame_sets
565 If `True`, decompose the transform and try to reference component
566 mappings that were already serialized into a `FrameSet` in the
567 archive. Note that if multiple transforms exist between a pair of
568 frames (e.g. a `SkyProjection` and its FITS approximation), this
569 may cause the wrong one to be saved. When this option is used, the
570 frame set must be saved before the transform, and it must be
571 deserialized before the transform as well.
573 Returns
574 -------
575 `TransformSerializationModel`
576 Serialized form of the transform.
577 """
578 model = TransformSerializationModel[P]()
579 if use_frame_sets: 579 ↛ 580line 579 didn't jump to line 580 because the condition on line 579 was never true
580 for link in self.decompose():
581 model.frames.append(link.in_frame.serialize())
582 model.bounds.append(link.in_bounds.serialize() if link.in_bounds is not None else None)
583 for frame_set, pointer in archive.iter_frame_sets():
584 if link.in_frame in frame_set and link.out_frame in frame_set:
585 model.mappings.append(pointer)
586 break
587 else:
588 model.mappings.append(MappingSerializationModel(ast=link._ast_mapping.show()))
589 else:
590 model.frames.append(self.in_frame.serialize())
591 model.bounds.append(self.in_bounds.serialize() if self.in_bounds is not None else None)
592 model.mappings.append(MappingSerializationModel(ast=self._ast_mapping.show()))
593 model.frames.append(self.out_frame.serialize())
594 model.bounds.append(self.out_bounds.serialize() if self.out_bounds is not None else None)
595 return model
597 @staticmethod
598 def _get_archive_tree_type[P: pydantic.BaseModel](
599 pointer_type: type[P],
600 ) -> type[TransformSerializationModel[P]]:
601 """Return the serialization model type for this object for an archive
602 type that uses the given pointer type.
603 """
604 return TransformSerializationModel[pointer_type] # type: ignore
606 @staticmethod
607 def from_legacy(
608 legacy: LegacyTransform,
609 in_frame: I,
610 out_frame: O,
611 in_bounds: Bounds | None = None,
612 out_bounds: Bounds | None = None,
613 ) -> Transform[I, O]:
614 """Construct a transform from a legacy `lsst.afw.geom.Transform`.
616 Parameters
617 ----------
618 legacy : `lsst.afw.geom.Transform`
619 Legacy transform object.
620 in_frame
621 Coordinate frame for input points to the forward transform.
622 out_frame
623 Coordinate frame for output points from the forward transform.
624 in_bounds
625 The region that bounds valid input points.
626 out_bounds
627 The region that bounds valid output points.
628 """
629 return Transform(
630 in_frame,
631 out_frame,
632 legacy.getMapping(),
633 in_bounds=in_bounds,
634 out_bounds=out_bounds,
635 )
637 def to_legacy(self) -> LegacyTransform:
638 """Convert to a legacy `lsst.afw.geom.TransformPoint2ToPoint2`
639 instance.
640 """
641 from lsst.afw.geom import TransformPoint2ToPoint2 as LegacyTransform
643 return LegacyTransform(self._ast_mapping, False)
645 def _get_ast_frame_set(self) -> Any:
646 ast_frame_set = astshim.FrameSet(_make_ast_frame(self._in_frame))
647 ast_frame_set.addFrame(astshim.FrameSet.BASE, self._ast_mapping, _make_ast_frame(self._out_frame))
648 return ast_frame_set
651def _ast_apply(method: Any, *, x: Any, y: Any) -> XY[float] | XY[np.ndarray]:
652 # TODO: add bounds argument and check inputs
653 xa = np.asarray(x)
654 ya = np.asarray(y)
655 broadcast_shape = np.broadcast(xa, ya).shape
656 scalar = not broadcast_shape
657 xb, yb = np.broadcast_arrays(xa, ya)
658 xy_in = np.vstack([xb.ravel(), yb.ravel()]).astype(np.float64)
659 xy_out = method(xy_in)
660 if scalar:
661 return XY(float(xy_out[0, 0]), float(xy_out[1, 0]))
662 return XY(xy_out[0].reshape(broadcast_shape), xy_out[1].reshape(broadcast_shape))
665def _prepend_ast_shift(ast_frame_set: Any, x: float, y: float, ast_domain: str) -> None:
666 ast_output_frame_id = ast_frame_set.current
667 ast_frame_set.addFrame(
668 astshim.FrameSet.BASE,
669 astshim.ShiftMap([x, y]),
670 astshim.Frame(2, f"Domain={ast_domain}"),
671 )
672 ast_frame_set.base = ast_frame_set.current
673 ast_frame_set.current = ast_output_frame_id
676def _make_ast_frame(frame: Frame) -> Any:
677 if frame is SkyFrame.ICRS:
678 return astshim.SkyFrame("")
679 ast_frame = astshim.Frame(2, f"Ident={frame._ast_ident}")
680 if frame.unit is not None: 680 ↛ 684line 680 didn't jump to line 684 because the condition on line 680 was always true
681 fits_unit = frame.unit.to_string(format="fits")
682 ast_frame.setUnit(1, fits_unit)
683 ast_frame.setUnit(2, fits_unit)
684 ast_frame.setLabel(1, "x")
685 ast_frame.setLabel(2, "y")
686 return ast_frame
689def _standardize_xy(xy: XY[Any], frame: Frame) -> XY[Any]:
690 return XY(x=frame.standardize_x(xy.x), y=frame.standardize_y(xy.y))
693class MappingSerializationModel(pydantic.BaseModel):
694 """Serialization model for an AST Mapping."""
696 ast: str = pydantic.Field(description="A serialized Starlink AST Mapping, using the AST native encoding.")
699class TransformSerializationModel[P: pydantic.BaseModel](ArchiveTree):
700 """Serialization model for coordinate transforms."""
702 SCHEMA_NAME: ClassVar[str] = "transform"
703 SCHEMA_VERSION: ClassVar[str] = "1.0.0"
704 MIN_READ_VERSION: ClassVar[int] = 1
705 PUBLIC_TYPE: ClassVar[type] = Transform
707 frames: list[SerializableFrame] = pydantic.Field(
708 default_factory=list,
709 description=textwrap.dedent(
710 """
711 List of frames that this transform passes through.
713 All transforms include at least two frames (the endpoints). Others
714 intermediate frames may be included to facilitate data-sharing
715 between transforms.
716 """
717 ),
718 )
720 bounds: list[BoundsSerializationModel | None] = pydantic.Field(
721 default_factory=list,
722 description=textwrap.dedent(
723 """
724 List of the bounds of the ``frames`` for this transform.
726 This always has the same number of elements as ``frames``.
727 """
728 ),
729 )
731 mappings: list[P | MappingSerializationModel] = pydantic.Field(
732 default_factory=list,
733 description=textwrap.dedent(
734 """
735 The actual mappings between frames, or archive pointers to
736 serialized FrameSet objects from which they can be obtained.
738 This always has one fewer element than ``frames``.
739 """
740 ),
741 )
743 def deserialize(self, archive: InputArchive[P], **kwargs: Any) -> Transform[Any, Any]:
744 """Deserialize a transform from an archive.
746 Parameters
747 ----------
748 archive
749 Archive to read from.
750 **kwargs
751 Unsupported keyword arguments are accepted only to provide better
752 error messages (raising `serialization.InvalidParameterError`).
753 """
754 if kwargs: 754 ↛ 755line 754 didn't jump to line 755 because the condition on line 754 was never true
755 raise InvalidParameterError(f"Unrecognized parameters for Transform: {set(kwargs.keys())}.")
756 if len(self.frames) != len(self.bounds): 756 ↛ 757line 756 didn't jump to line 757 because the condition on line 756 was never true
757 raise ArchiveReadError(
758 f"Inconsistent lengths for 'frames' ({len(self.frames)}) and 'bounds' ({len(self.bounds)})."
759 )
760 if len(self.frames) != len(self.mappings) + 1: 760 ↛ 761line 760 didn't jump to line 761 because the condition on line 760 was never true
761 raise ArchiveReadError(
762 f"Inconsistent lengths for 'frames' ({len(self.frames)}) and "
763 f"'mappings' ({len(self.mappings)}; should be one less)."
764 )
765 # We can't just compose onto an identity Transform if we want to
766 # preserve the FrameSet-ness of any of these mappings.
767 transform: Transform | None = None
768 for n, mapping in enumerate(self.mappings):
769 match mapping:
770 case MappingSerializationModel(ast=serialized_mapping): 770 ↛ 781line 770 didn't jump to line 781 because the pattern on line 770 always matched
771 ast_mapping = astshim.Mapping.fromString(serialized_mapping)
772 in_bounds = self.bounds[n]
773 out_bounds = self.bounds[n + 1]
774 new_transform = Transform(
775 self.frames[n].deserialize(),
776 self.frames[n + 1].deserialize(),
777 ast_mapping,
778 in_bounds.deserialize() if in_bounds is not None else None,
779 out_bounds.deserialize() if out_bounds is not None else None,
780 )
781 case reference:
782 frame_set = archive.get_frame_set(reference)
783 new_transform = frame_set[self.frames[n].deserialize(), self.frames[n + 1].deserialize()]
784 if transform is None: 784 ↛ 787line 784 didn't jump to line 787 because the condition on line 784 was always true
785 transform = new_transform
786 else:
787 transform = transform.then(new_transform)
788 if transform is None: 788 ↛ 789line 788 didn't jump to line 789 because the condition on line 788 was never true
789 transform = Transform.identity(self.frames[0].deserialize())
790 return transform
793if TYPE_CHECKING:
795 def _test_types() -> None:
796 t = cast(Transform, None)
797 arr = np.zeros(3)
799 # Scalar inputs → XY[float]
800 assert_type(t.apply_forward(x=1.0, y=2.0), XY[float])
801 assert_type(t.apply_inverse(x=1.0, y=2.0), XY[float])
803 # Array inputs → XY[np.ndarray]
804 assert_type(t.apply_forward(x=arr, y=arr), XY[np.ndarray])
805 assert_type(t.apply_inverse(x=arr, y=arr), XY[np.ndarray])
807 # Array-like (list) inputs → XY[np.ndarray]
808 assert_type(t.apply_forward(x=[1.0, 2.0], y=[3.0, 4.0]), XY[np.ndarray])
809 assert_type(t.apply_inverse(x=[1.0, 2.0], y=[3.0, 4.0]), XY[np.ndarray])