Coverage for python/lsst/images/_transforms/_sky_projection.py: 82%
198 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__ = ("SkyProjection", "SkyProjectionAstropyView", "SkyProjectionSerializationModel")
16import functools
17from typing import TYPE_CHECKING, Any, ClassVar, Self, TypeVar, assert_type, cast, final, overload
19import astropy.units as u
20import astropy.wcs
21import numpy as np
22import numpy.typing as npt
23import pydantic
24from astropy.coordinates import ICRS, Latitude, Longitude, SkyCoord
25from astropy.wcs.wcsapi import BaseLowLevelWCS, HighLevelWCSMixin
27from .._geom import XY, YX, Bounds, Box
28from ..serialization import ArchiveTree, InputArchive, InvalidParameterError, OutputArchive
29from ..utils import is_none
30from . import _ast as astshim
31from ._frames import Frame, SkyFrame
32from ._transform import Transform, TransformSerializationModel, _ast_apply
34if TYPE_CHECKING:
35 try:
36 from lsst.afw.geom import SkyWcs as LegacySkyWcs
37 except ImportError:
38 type LegacySkyWcs = Any # type: ignore[no-redef]
41# This pre-python-3.12 declaration is needed by Sphinx (probably the
42# autodoc-typehints plugin.
43F = TypeVar("F", bound=Frame)
44P = TypeVar("P", bound=pydantic.BaseModel)
47def _set_ast_skyframe_system(frame: astshim.SkyFrame, system: str) -> None:
48 """Set an AST SkyFrame coordinate system across supported wrappers."""
49 if hasattr(frame, "_impl"): 49 ↛ 50line 49 didn't jump to line 50 because the condition on line 49 was never true
50 frame._impl.System = system
51 else:
52 setattr(frame, "system", system)
55@final
56class SkyProjection[F: Frame]:
57 """A transform from pixel coordinates to sky coordinates.
59 Parameters
60 ----------
61 pixel_to_sky
62 A low-level transform that maps pixel coordinates to sky coordinates.
63 fits_approximation
64 An approximation to ``pixel_to_sky`` that is guaranteed to have a
65 `~Transform.as_fits_wcs` method that does not return `None`. This
66 should not be provided if ``pixel_to_sky`` is itself representable
67 as a FITS WCS.
69 Notes
70 -----
71 `Transform` is conceptually immutable (the internal AST Mapping should
72 never be modified in-place after construction), and hence does not need to
73 be copied when any object that holds it is copied.
74 """
76 def __init__(
77 self, pixel_to_sky: Transform[F, SkyFrame], fits_approximation: Transform[F, SkyFrame] | None = None
78 ) -> None:
79 self._pixel_to_sky = pixel_to_sky
80 if pixel_to_sky.in_frame.unit != u.pix:
81 raise ValueError("Transform is not a mapping from pixel coordinates.")
82 if pixel_to_sky.out_frame != SkyFrame.ICRS:
83 raise ValueError("Transform is not a mapping to ICRS.")
84 self._fits_approximation = fits_approximation
86 def __eq__(self, other: Any) -> bool:
87 if self is other:
88 return True
89 if not isinstance(other, SkyProjection):
90 return NotImplemented
91 # Even though two approximations could be different and yet consistent
92 # with the primary mapping (for example using different tolerances
93 # on construction) we require them to be equal to declare that the
94 # two objects are equal.
95 if self._fits_approximation != other._fits_approximation:
96 return False
97 return self._pixel_to_sky == other._pixel_to_sky
99 @staticmethod
100 def from_fits_wcs(
101 fits_wcs: astropy.wcs.WCS,
102 pixel_frame: F,
103 pixel_bounds: Bounds | None = None,
104 x0: int = 0,
105 y0: int = 0,
106 ) -> SkyProjection[F]:
107 """Construct a transform from a FITS WCS.
109 Parameters
110 ----------
111 fits_wcs
112 FITS WCS to convert.
113 pixel_frame
114 Coordinate frame for the pixel grid.
115 pixel_bounds
116 The region that bounds valid pixels for this transform.
117 x0
118 Logical coordinate of the first column in the array this WCS
119 relates to world coordinates.
120 y0
121 Logical coordinate of the first column in the array this WCS
122 relates to world coordinates.
124 Notes
125 -----
126 The ``x0`` and ``y0`` parameters reflect the fact that for FITS, the
127 first row and column are always labeled ``(1, 1)``, while in Astropy
128 and most other Python libraries, they are ``(0, 0)``. The `types` in
129 this package (e.g. `Image`, `Mask`) allow them to be any pair of
130 integers.
132 See Also
133 --------
134 Transform.from_fits_wcs
135 """
136 return SkyProjection(
137 Transform.from_fits_wcs(
138 fits_wcs, pixel_frame, SkyFrame.ICRS, in_bounds=pixel_bounds, x0=x0, y0=y0
139 )
140 )
142 @staticmethod
143 def from_ast_frame_set(
144 ast_frame_set: astshim.FrameSet,
145 pixel_frame: F,
146 pixel_bounds: Bounds | None = None,
147 ) -> SkyProjection[F]:
148 """Construct a sky projection from an AST FrameSet.
150 The current frame of the FrameSet must be an AST SkyFrame. Its
151 coordinate system is forced to ICRS (AST adjusts the mapping
152 automatically) so the resulting Projection is always in ICRS
153 regardless of the original sky system.
155 Parameters
156 ----------
157 ast_frame_set
158 An AST FrameSet whose base frame is pixel coordinates and
159 whose current frame is a SkyFrame (in any supported sky
160 coordinate system).
161 pixel_frame
162 Coordinate frame for the pixel grid.
163 pixel_bounds
164 The region that bounds valid pixels for this transform.
166 Raises
167 ------
168 ValueError
169 If the current frame of the FrameSet is not a SkyFrame.
170 """
171 current_frame = ast_frame_set.getFrame(ast_frame_set.current, copy=False)
172 if not isinstance(current_frame, astshim.SkyFrame): 172 ↛ 173line 172 didn't jump to line 173 because the condition on line 172 was never true
173 raise ValueError(
174 "The current frame of the AST FrameSet is not a SkyFrame "
175 f"(got {type(current_frame).__name__})."
176 )
177 _set_ast_skyframe_system(current_frame, "ICRS")
178 return SkyProjection(Transform(pixel_frame, SkyFrame.ICRS, ast_frame_set, in_bounds=pixel_bounds))
180 @property
181 def pixel_frame(self) -> F:
182 """Coordinate frame for the pixel grid."""
183 return self._pixel_to_sky.in_frame
185 @property
186 def sky_frame(self) -> SkyFrame:
187 """Coordinate frame for the sky."""
188 return self._pixel_to_sky.out_frame
190 @property
191 def pixel_bounds(self) -> Bounds | None:
192 """The region that bounds valid pixel points (`Bounds` | `None`)."""
193 return self._pixel_to_sky.in_bounds
195 @property
196 def pixel_to_sky_transform(self) -> Transform[F, SkyFrame]:
197 """Low-level transform from pixel to sky coordinates (`Transform`)."""
198 return self._pixel_to_sky
200 @property
201 def sky_to_pixel_transform(self) -> Transform[SkyFrame, F]:
202 """Low-level transform from sky to pixel coordinates (`Transform`)."""
203 return self._pixel_to_sky.inverted()
205 @property
206 def fits_approximation(self) -> SkyProjection[F] | None:
207 """An approximation to this projection that is guaranteed to have an
208 `as_fits_wcs` method that does not return `None`.
209 """
210 return SkyProjection(self._fits_approximation) if self._fits_approximation is not None else None
212 def show(self, simplified: bool = False, comments: bool = False) -> str:
213 """Return the AST native representation of the transform.
215 Parameters
216 ----------
217 simplified
218 Whether to ask AST to simplify the mapping before showing it.
219 This will make it much more likely that two equivalent transforms
220 have the same `show` result.
221 comments
222 Whether to include descriptive comments.
223 """
224 return self._pixel_to_sky.show(simplified=simplified, comments=comments)
226 @overload
227 def pixel_to_sky(self, point: XY[int | float] | YX[int | float], /) -> SkyCoord: ... 227 ↛ exitline 227 didn't return from function 'pixel_to_sky' because
229 @overload
230 def pixel_to_sky(self, point: XY[npt.ArrayLike] | YX[npt.ArrayLike], /) -> SkyCoord: ... 230 ↛ exitline 230 didn't return from function 'pixel_to_sky' because
232 @overload
233 def pixel_to_sky( 233 ↛ exitline 233 didn't return from function 'pixel_to_sky' because
234 self, /, *, x: int | float | npt.ArrayLike, y: int | float | npt.ArrayLike
235 ) -> SkyCoord: ...
237 def pixel_to_sky(
238 self,
239 point: XY[Any] | YX[Any] | None = None,
240 /,
241 *,
242 x: Any = None,
243 y: Any = None,
244 ) -> SkyCoord:
245 """Transform one or more pixel points to sky coordinates.
247 Parameters
248 ----------
249 point
250 An `XY` or `YX` coordinate pair of pixel positions to transform.
251 Mutually exclusive with ``x`` and ``y``.
252 x : `float` | array-like
253 ``x`` values of the pixel points to transform, as a scalar or
254 any array-like. Results are broadcast against ``y``.
255 Mutually exclusive with ``point``.
256 y : `float` | array-like
257 ``y`` values of the pixel points to transform, as a scalar or
258 any array-like. Results are broadcast against ``x``.
259 Mutually exclusive with ``point``.
261 Returns
262 -------
263 astropy.coordinates.SkyCoord
264 Transformed sky coordinates, with the broadcast shape of ``x``
265 and ``y``.
266 """
267 match point:
268 case None:
269 if x is None or y is None:
270 raise TypeError("Pass either a point or both x= and y= to 'pixel_to_sky'.")
271 case XY() | YX(): 271 ↛ 275line 271 didn't jump to line 275 because the pattern on line 271 always matched
272 if x is not None or y is not None:
273 raise TypeError("'pixel_to_sky' point argument is mutually exclusive with x= and y=.")
274 x, y = point.x, point.y
275 case _:
276 raise TypeError(f"Unexpected positional argument type: {type(point)!r}.")
277 sky_rad = self._pixel_to_sky.apply_forward(x=x, y=y)
278 return SkyCoord(ra=sky_rad.x, dec=sky_rad.y, unit=u.rad)
280 def sky_to_pixel(self, sky: SkyCoord) -> XY[Any]:
281 """Transform one or more sky coordinates to pixels.
283 Parameters
284 ----------
285 sky
286 Sky coordinates to transform. Any shape is supported; the
287 result has the same shape as ``sky``.
289 Returns
290 -------
291 `XY` [`numpy.ndarray` | `float`]
292 Transformed pixel coordinates with the same shape as ``sky``.
293 """
294 if sky.frame.name != "icrs": 294 ↛ 295line 294 didn't jump to line 295 because the condition on line 294 was never true
295 sky = sky.transform_to("icrs")
296 ra: Longitude = sky.ra
297 dec: Latitude = sky.dec
298 # cast works around a mypy false positive specific to generic
299 # NamedTuple classes: returning XY[Any] from an overloaded function
300 # when called with Any-typed arguments produces "Incompatible return
301 # value type (got XY[Any], expected XY[Any])".
302 return cast(
303 XY[Any],
304 self._pixel_to_sky.apply_inverse(
305 x=ra.to_value(u.rad),
306 y=dec.to_value(u.rad),
307 ),
308 )
310 def as_astropy(self, bbox: Box | None = None) -> SkyProjectionAstropyView:
311 """Return an `astropy.wcs` view of this `SkyProjection`.
313 Parameters
314 ----------
315 bbox
316 Bounding box of the array the view will describe. This
317 projection object is assumed to work on the same coordinate system
318 in which ``bbox`` is defined, while the Astropy view will consider
319 the first row and column in that box to be ``(0, 0)``.
321 Notes
322 -----
323 This returns an object that satisfies the
324 `astropy.wcs.wcsapi.BaseHighLevelWCS` and
325 `astropy.wcs.wcsapi.BaseLowLevelWCS` interfaces while evaluating the
326 underlying `SkyProjection` itself. It is *not* an `astropy.wcs.WCS`
327 instance, which is a type that also satisfies those interfaces but
328 only supports FITS WCS representations (see `as_fits_wcs`).
329 """
330 return SkyProjectionAstropyView(self._pixel_to_sky._ast_mapping, bbox)
332 def as_fits_wcs(self, bbox: Box, allow_approximation: bool = False) -> astropy.wcs.WCS | None:
333 """Return a FITS WCS representation of this projection, if possible.
335 Parameters
336 ----------
337 bbox
338 Bounding box of the array the FITS WCS will describe. This
339 transform object is assumed to work on the same coordinate system
340 in which ``bbox`` is defined, while the FITS WCS will consider the
341 first row and column in that box to be ``(0, 0)`` (in Astropy
342 interfaces) or ``(1, 1)`` (in the FITS representation itself).
343 allow_approximation
344 If `True` and this `SkyProjection` holds a FITS approximation to
345 itself, return that approximation.
346 """
347 if allow_approximation and self._fits_approximation: 347 ↛ 348line 347 didn't jump to line 348 because the condition on line 347 was never true
348 return self._fits_approximation.as_fits_wcs(bbox)
349 return self._pixel_to_sky.as_fits_wcs(bbox)
351 def serialize[P: pydantic.BaseModel](
352 self, archive: OutputArchive[P], *, use_frame_sets: bool = False
353 ) -> SkyProjectionSerializationModel[P]:
354 """Serialize a projection to an archive.
356 Parameters
357 ----------
358 archive
359 Archive to serialize to.
360 use_frame_sets
361 If `True`, decompose the underlying transform and try to reference
362 component mappings that were already serialized into a `FrameSet`
363 in the archive. The FITS approximation transform is never
364 decomposed.
366 Returns
367 -------
368 `SkyProjectionSerializationModel`
369 Serialized form of the projection.
370 """
371 pixel_to_sky = archive.serialize_direct(
372 "pixel_to_sky", functools.partial(self._pixel_to_sky.serialize, use_frame_sets=use_frame_sets)
373 )
374 fits_approximation = (
375 archive.serialize_direct("fits_approximation", self._fits_approximation.serialize)
376 if self._fits_approximation is not None
377 else None
378 )
379 return SkyProjectionSerializationModel(
380 pixel_to_sky=pixel_to_sky, fits_approximation=fits_approximation
381 )
383 @staticmethod
384 def _get_archive_tree_type[P: pydantic.BaseModel](
385 pointer_type: type[P],
386 ) -> type[SkyProjectionSerializationModel[P]]:
387 """Return the serialization model type for this object for an archive
388 type that uses the given pointer type.
389 """
390 return SkyProjectionSerializationModel[pointer_type] # type: ignore
392 @staticmethod
393 def from_legacy(
394 sky_wcs: LegacySkyWcs, pixel_frame: F, pixel_bounds: Bounds | None = None
395 ) -> SkyProjection[F]:
396 """Construct a transform from a legacy `lsst.afw.geom.SkyWcs`.
398 Parameters
399 ----------
400 sky_wcs : `lsst.afw.geom.SkyWcs`
401 Legacy WCS object.
402 pixel_frame
403 Coordinate frame for the pixel grid.
404 pixel_bounds
405 The region that bounds valid pixels for this transform.
406 """
407 fits_approximation: Transform[F, SkyFrame] | None = None
408 if (legacy_fits_approximation := sky_wcs.getFitsApproximation()) is not None:
409 fits_approximation = Transform(
410 pixel_frame,
411 SkyFrame.ICRS,
412 legacy_fits_approximation.getFrameDict(),
413 pixel_bounds,
414 )
415 return SkyProjection(
416 Transform(pixel_frame, SkyFrame.ICRS, sky_wcs.getFrameDict(), pixel_bounds),
417 fits_approximation=fits_approximation,
418 )
420 def to_legacy(self) -> LegacySkyWcs:
421 """Convert to a legacy `lsst.afw.geom.SkyWcs` instance."""
422 from lsst.afw.geom import SkyWcs as LegacySkyWcs
424 try:
425 ast_mapping = astshim.FrameDict(self._pixel_to_sky._ast_mapping)
426 except TypeError as err:
427 err.add_note(
428 "Only Projections created by from_legacy and from_fits_wcs "
429 "are guaranteed to be convertible to SkyWcs."
430 )
431 raise
432 # SkyWcs requires the pixel frame's domain to be PIXELS, while AST
433 # (and hence NDF and from_fits_wcs) uses PIXEL, so rename it in the
434 # FrameDict (a deep copy, so this projection itself is unchanged).
435 if ast_mapping.hasDomain("PIXEL") and not ast_mapping.hasDomain("PIXELS"): 435 ↛ 440line 435 didn't jump to line 440 because the condition on line 435 was always true
436 saved_current = ast_mapping.current
437 ast_mapping.setCurrent("PIXEL")
438 ast_mapping.setDomain("PIXELS")
439 ast_mapping.current = saved_current
440 legacy_wcs = LegacySkyWcs(ast_mapping)
441 if self.fits_approximation is not None: 441 ↛ 442line 441 didn't jump to line 442 because the condition on line 441 was never true
442 legacy_wcs = legacy_wcs.copyWithFitsApproximation(self.fits_approximation.to_legacy())
443 return legacy_wcs
446class SkyProjectionAstropyView(BaseLowLevelWCS, HighLevelWCSMixin):
447 """An Astropy-interface view of a `SkyProjection`.
449 Parameters
450 ----------
451 ast_pixel_to_sky
452 AST mapping from pixel coordinates to sky coordinates.
453 bbox
454 Bounding box of the projection, or `None` if unbounded.
456 Notes
457 -----
458 The constructor of this classe is considered a private implementation
459 detail; use `SkyProjection.as_astropy` instead.
461 This object satisfies the `astropy.wcs.wcsapi.BaseHighLevelWCS` and
462 `astropy.wcs.wcsapi.BaseLowLevelWCS` interfaces while evaluating the
463 underlying `SkyProjection` itself. It is *not* an `astropy.wcs.WCS`
464 subclass, which is a type that also satisfies those interfaces but
465 only supports FITS WCS representations (see `SkyProjection.as_fits_wcs`).
466 """
468 def __init__(self, ast_pixel_to_sky: astshim.Mapping, bbox: Box | None) -> None:
469 self._bbox = bbox
470 if bbox is not None:
471 ast_pixel_to_sky = astshim.ShiftMap(list(bbox.start.xy)).then(ast_pixel_to_sky)
472 self._ast_pixel_to_sky = ast_pixel_to_sky
474 @property
475 def low_level_wcs(self) -> Self:
476 return self
478 @property
479 def array_shape(self) -> YX[int] | None:
480 return self._bbox.shape if self._bbox is not None else None
482 @property
483 def axis_correlation_matrix(self) -> np.ndarray:
484 return np.array([[True, True], [True, True]])
486 @property
487 def pixel_axis_names(self) -> XY[str]:
488 return XY("x", "y")
490 @property
491 def pixel_bounds(self) -> XY[tuple[int, int]] | None:
492 if self._bbox is None:
493 return None
494 return XY((self._bbox.x.min, self._bbox.x.max), (self._bbox.y.min, self._bbox.y.max))
496 @property
497 def pixel_n_dim(self) -> int:
498 return 2
500 @property
501 def pixel_shape(self) -> XY[int] | None:
502 array_shape = self.array_shape
503 return array_shape.xy if array_shape is not None else None
505 @property
506 def serialized_classes(self) -> bool:
507 return False
509 @property
510 def world_axis_names(self) -> tuple[str, str]:
511 return ("ra", "dec")
513 @property
514 def world_axis_object_classes(self) -> dict[str, tuple[type[SkyCoord], tuple[()], dict[str, Any]]]:
515 return {"celestial": (SkyCoord, (), {"frame": ICRS, "unit": (u.rad, u.rad)})}
517 @property
518 def world_axis_object_components(self) -> list[tuple[str, int, str]]:
519 return [("celestial", 0, "spherical.lon.radian"), ("celestial", 1, "spherical.lat.radian")]
521 @property
522 def world_axis_physical_types(self) -> tuple[str, str]:
523 return ("pos.eq.ra", "pos.eq.dec")
525 @property
526 def world_axis_units(self) -> tuple[str, str]:
527 return ("rad", "rad")
529 @property
530 def world_n_dim(self) -> int:
531 return 2
533 def pixel_to_world_values(self, x: np.ndarray, y: np.ndarray) -> XY[Any]:
534 return _ast_apply(self._ast_pixel_to_sky.applyForward, x=x, y=y)
536 def world_to_pixel_values(self, ra: np.ndarray, dec: np.ndarray) -> XY[Any]:
537 return _ast_apply(self._ast_pixel_to_sky.applyInverse, x=ra, y=dec)
540class SkyProjectionSerializationModel[P: pydantic.BaseModel](ArchiveTree):
541 """Serialization model for projetions."""
543 SCHEMA_NAME: ClassVar[str] = "sky_projection"
544 SCHEMA_VERSION: ClassVar[str] = "1.0.0"
545 MIN_READ_VERSION: ClassVar[int] = 1
546 PUBLIC_TYPE: ClassVar[type] = SkyProjection
548 pixel_to_sky: TransformSerializationModel[P] = pydantic.Field(
549 description="The transform that maps pixel coordinates to the sky."
550 )
551 fits_approximation: TransformSerializationModel[P] | None = pydantic.Field(
552 default=None,
553 description=(
554 "An approximation of the pixel-to-sky transform that is exactly representable as a FITS WCS."
555 ),
556 exclude_if=is_none,
557 )
559 def deserialize(self, archive: InputArchive[P], **kwargs: Any) -> SkyProjection[Any]:
560 """Deserialize a projection from an archive.
562 Parameters
563 ----------
564 archive
565 Archive to read from.
566 **kwargs
567 Unsupported keyword arguments are accepted only to provide better
568 error messages (raising `serialization.InvalidParameterError`).
569 """
570 if kwargs: 570 ↛ 571line 570 didn't jump to line 571 because the condition on line 570 was never true
571 raise InvalidParameterError(f"Unrecognized parameters for Projection: {set(kwargs.keys())}.")
572 pixel_to_sky = self.pixel_to_sky.deserialize(archive)
573 fits_approximation = (
574 self.fits_approximation.deserialize(archive) if self.fits_approximation is not None else None
575 )
576 return SkyProjection(pixel_to_sky, fits_approximation=fits_approximation)
579if TYPE_CHECKING:
581 def _test_types() -> None:
582 sp = cast(SkyProjection, None)
583 sky = cast(SkyCoord, None)
585 # sky_to_pixel returns XY[Any] so that callers need no cast regardless
586 # of whether sky is scalar or array-shaped.
587 assert_type(sp.sky_to_pixel(sky), XY[Any])