Coverage for python/lsst/images/_transforms/_sky_projection.py: 75%
193 statements
« prev ^ index » next coverage.py v7.14.3, created at 2026-07-15 02:03 -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 legacy_wcs = LegacySkyWcs(ast_mapping)
433 if self.fits_approximation is not None:
434 legacy_wcs = legacy_wcs.copyWithFitsApproximation(self.fits_approximation.to_legacy())
435 return legacy_wcs
438class SkyProjectionAstropyView(BaseLowLevelWCS, HighLevelWCSMixin):
439 """An Astropy-interface view of a `SkyProjection`.
441 Parameters
442 ----------
443 ast_pixel_to_sky
444 AST mapping from pixel coordinates to sky coordinates.
445 bbox
446 Bounding box of the projection, or `None` if unbounded.
448 Notes
449 -----
450 The constructor of this classe is considered a private implementation
451 detail; use `SkyProjection.as_astropy` instead.
453 This object satisfies the `astropy.wcs.wcsapi.BaseHighLevelWCS` and
454 `astropy.wcs.wcsapi.BaseLowLevelWCS` interfaces while evaluating the
455 underlying `SkyProjection` itself. It is *not* an `astropy.wcs.WCS`
456 subclass, which is a type that also satisfies those interfaces but
457 only supports FITS WCS representations (see `SkyProjection.as_fits_wcs`).
458 """
460 def __init__(self, ast_pixel_to_sky: astshim.Mapping, bbox: Box | None) -> None:
461 self._bbox = bbox
462 if bbox is not None:
463 ast_pixel_to_sky = astshim.ShiftMap(list(bbox.start.xy)).then(ast_pixel_to_sky)
464 self._ast_pixel_to_sky = ast_pixel_to_sky
466 @property
467 def low_level_wcs(self) -> Self:
468 return self
470 @property
471 def array_shape(self) -> YX[int] | None:
472 return self._bbox.shape if self._bbox is not None else None
474 @property
475 def axis_correlation_matrix(self) -> np.ndarray:
476 return np.array([[True, True], [True, True]])
478 @property
479 def pixel_axis_names(self) -> XY[str]:
480 return XY("x", "y")
482 @property
483 def pixel_bounds(self) -> XY[tuple[int, int]] | None:
484 if self._bbox is None:
485 return None
486 return XY((self._bbox.x.min, self._bbox.x.max), (self._bbox.y.min, self._bbox.y.max))
488 @property
489 def pixel_n_dim(self) -> int:
490 return 2
492 @property
493 def pixel_shape(self) -> XY[int] | None:
494 array_shape = self.array_shape
495 return array_shape.xy if array_shape is not None else None
497 @property
498 def serialized_classes(self) -> bool:
499 return False
501 @property
502 def world_axis_names(self) -> tuple[str, str]:
503 return ("ra", "dec")
505 @property
506 def world_axis_object_classes(self) -> dict[str, tuple[type[SkyCoord], tuple[()], dict[str, Any]]]:
507 return {"celestial": (SkyCoord, (), {"frame": ICRS, "unit": (u.rad, u.rad)})}
509 @property
510 def world_axis_object_components(self) -> list[tuple[str, int, str]]:
511 return [("celestial", 0, "spherical.lon.radian"), ("celestial", 1, "spherical.lat.radian")]
513 @property
514 def world_axis_physical_types(self) -> tuple[str, str]:
515 return ("pos.eq.ra", "pos.eq.dec")
517 @property
518 def world_axis_units(self) -> tuple[str, str]:
519 return ("rad", "rad")
521 @property
522 def world_n_dim(self) -> int:
523 return 2
525 def pixel_to_world_values(self, x: np.ndarray, y: np.ndarray) -> XY[Any]:
526 return _ast_apply(self._ast_pixel_to_sky.applyForward, x=x, y=y)
528 def world_to_pixel_values(self, ra: np.ndarray, dec: np.ndarray) -> XY[Any]:
529 return _ast_apply(self._ast_pixel_to_sky.applyInverse, x=ra, y=dec)
532class SkyProjectionSerializationModel[P: pydantic.BaseModel](ArchiveTree):
533 """Serialization model for projetions."""
535 SCHEMA_NAME: ClassVar[str] = "sky_projection"
536 SCHEMA_VERSION: ClassVar[str] = "1.0.0"
537 MIN_READ_VERSION: ClassVar[int] = 1
538 PUBLIC_TYPE: ClassVar[type] = SkyProjection
540 pixel_to_sky: TransformSerializationModel[P] = pydantic.Field(
541 description="The transform that maps pixel coordinates to the sky."
542 )
543 fits_approximation: TransformSerializationModel[P] | None = pydantic.Field(
544 default=None,
545 description=(
546 "An approximation of the pixel-to-sky transform that is exactly representable as a FITS WCS."
547 ),
548 exclude_if=is_none,
549 )
551 def deserialize(self, archive: InputArchive[P], **kwargs: Any) -> SkyProjection[Any]:
552 """Deserialize a projection from an archive.
554 Parameters
555 ----------
556 archive
557 Archive to read from.
558 **kwargs
559 Unsupported keyword arguments are accepted only to provide better
560 error messages (raising `serialization.InvalidParameterError`).
561 """
562 if kwargs: 562 ↛ 563line 562 didn't jump to line 563 because the condition on line 562 was never true
563 raise InvalidParameterError(f"Unrecognized parameters for Projection: {set(kwargs.keys())}.")
564 pixel_to_sky = self.pixel_to_sky.deserialize(archive)
565 fits_approximation = (
566 self.fits_approximation.deserialize(archive) if self.fits_approximation is not None else None
567 )
568 return SkyProjection(pixel_to_sky, fits_approximation=fits_approximation)
571if TYPE_CHECKING:
573 def _test_types() -> None:
574 sp = cast(SkyProjection, None)
575 sky = cast(SkyCoord, None)
577 # sky_to_pixel returns XY[Any] so that callers need no cast regardless
578 # of whether sky is scalar or array-shaped.
579 assert_type(sp.sky_to_pixel(sky), XY[Any])