Coverage for python/lsst/images/_transforms/_sky_projection.py: 69%
176 statements
« prev ^ index » next coverage.py v7.14.3, created at 2026-07-02 08:55 +0000
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, final
19import astropy.units as u
20import astropy.wcs
21import numpy as np
22import pydantic
23from astropy.coordinates import ICRS, Latitude, Longitude, SkyCoord
24from astropy.wcs.wcsapi import BaseLowLevelWCS, HighLevelWCSMixin
26from .._geom import XY, YX, Bounds, Box
27from ..serialization import ArchiveTree, InputArchive, InvalidParameterError, OutputArchive
28from ..utils import is_none
29from . import _ast as astshim
30from ._frames import Frame, SkyFrame
31from ._transform import Transform, TransformSerializationModel, _ast_apply
33if TYPE_CHECKING:
34 try:
35 from lsst.afw.geom import SkyWcs as LegacySkyWcs
36 except ImportError:
37 type LegacySkyWcs = Any # type: ignore[no-redef]
40# This pre-python-3.12 declaration is needed by Sphinx (probably the
41# autodoc-typehints plugin.
42F = TypeVar("F", bound=Frame)
43P = TypeVar("P", bound=pydantic.BaseModel)
46def _set_ast_skyframe_system(frame: astshim.SkyFrame, system: str) -> None:
47 """Set an AST SkyFrame coordinate system across supported wrappers."""
48 if hasattr(frame, "_impl"): 48 ↛ 49line 48 didn't jump to line 49 because the condition on line 48 was never true
49 frame._impl.System = system
50 else:
51 setattr(frame, "system", system)
54@final
55class SkyProjection[F: Frame]:
56 """A transform from pixel coordinates to sky coordinates.
58 Parameters
59 ----------
60 pixel_to_sky
61 A low-level transform that maps pixel coordinates to sky coordinates.
62 fits_approximation
63 An approximation to ``pixel_to_sky`` that is guaranteed to have a
64 `~Transform.as_fits_wcs` method that does not return `None`. This
65 should not be provided if ``pixel_to_sky`` is itself representable
66 as a FITS WCS.
68 Notes
69 -----
70 `Transform` is conceptually immutable (the internal AST Mapping should
71 never be modified in-place after construction), and hence does not need to
72 be copied when any object that holds it is copied.
73 """
75 def __init__(
76 self, pixel_to_sky: Transform[F, SkyFrame], fits_approximation: Transform[F, SkyFrame] | None = None
77 ) -> None:
78 self._pixel_to_sky = pixel_to_sky
79 if pixel_to_sky.in_frame.unit != u.pix:
80 raise ValueError("Transform is not a mapping from pixel coordinates.")
81 if pixel_to_sky.out_frame != SkyFrame.ICRS:
82 raise ValueError("Transform is not a mapping to ICRS.")
83 self._fits_approximation = fits_approximation
85 def __eq__(self, other: Any) -> bool:
86 if self is other:
87 return True
88 if not isinstance(other, SkyProjection):
89 return NotImplemented
90 # Even though two approximations could be different and yet consistent
91 # with the primary mapping (for example using different tolerances
92 # on construction) we require them to be equal to declare that the
93 # two objects are equal.
94 if self._fits_approximation != other._fits_approximation:
95 return False
96 return self._pixel_to_sky == other._pixel_to_sky
98 @staticmethod
99 def from_fits_wcs(
100 fits_wcs: astropy.wcs.WCS,
101 pixel_frame: F,
102 pixel_bounds: Bounds | None = None,
103 x0: int = 0,
104 y0: int = 0,
105 ) -> SkyProjection[F]:
106 """Construct a transform from a FITS WCS.
108 Parameters
109 ----------
110 fits_wcs
111 FITS WCS to convert.
112 pixel_frame
113 Coordinate frame for the pixel grid.
114 pixel_bounds
115 The region that bounds valid pixels for this transform.
116 x0
117 Logical coordinate of the first column in the array this WCS
118 relates to world coordinates.
119 y0
120 Logical coordinate of the first column in the array this WCS
121 relates to world coordinates.
123 Notes
124 -----
125 The ``x0`` and ``y0`` parameters reflect the fact that for FITS, the
126 first row and column are always labeled ``(1, 1)``, while in Astropy
127 and most other Python libraries, they are ``(0, 0)``. The `types` in
128 this package (e.g. `Image`, `Mask`) allow them to be any pair of
129 integers.
131 See Also
132 --------
133 Transform.from_fits_wcs
134 """
135 return SkyProjection(
136 Transform.from_fits_wcs(
137 fits_wcs, pixel_frame, SkyFrame.ICRS, in_bounds=pixel_bounds, x0=x0, y0=y0
138 )
139 )
141 @staticmethod
142 def from_ast_frame_set(
143 ast_frame_set: astshim.FrameSet,
144 pixel_frame: F,
145 pixel_bounds: Bounds | None = None,
146 ) -> SkyProjection[F]:
147 """Construct a sky projection from an AST FrameSet.
149 The current frame of the FrameSet must be an AST SkyFrame. Its
150 coordinate system is forced to ICRS (AST adjusts the mapping
151 automatically) so the resulting Projection is always in ICRS
152 regardless of the original sky system.
154 Parameters
155 ----------
156 ast_frame_set
157 An AST FrameSet whose base frame is pixel coordinates and
158 whose current frame is a SkyFrame (in any supported sky
159 coordinate system).
160 pixel_frame
161 Coordinate frame for the pixel grid.
162 pixel_bounds
163 The region that bounds valid pixels for this transform.
165 Raises
166 ------
167 ValueError
168 If the current frame of the FrameSet is not a SkyFrame.
169 """
170 current_frame = ast_frame_set.getFrame(ast_frame_set.current, copy=False)
171 if not isinstance(current_frame, astshim.SkyFrame): 171 ↛ 172line 171 didn't jump to line 172 because the condition on line 171 was never true
172 raise ValueError(
173 "The current frame of the AST FrameSet is not a SkyFrame "
174 f"(got {type(current_frame).__name__})."
175 )
176 _set_ast_skyframe_system(current_frame, "ICRS")
177 return SkyProjection(Transform(pixel_frame, SkyFrame.ICRS, ast_frame_set, in_bounds=pixel_bounds))
179 @property
180 def pixel_frame(self) -> F:
181 """Coordinate frame for the pixel grid."""
182 return self._pixel_to_sky.in_frame
184 @property
185 def sky_frame(self) -> SkyFrame:
186 """Coordinate frame for the sky."""
187 return self._pixel_to_sky.out_frame
189 @property
190 def pixel_bounds(self) -> Bounds | None:
191 """The region that bounds valid pixel points (`Bounds` | `None`)."""
192 return self._pixel_to_sky.in_bounds
194 @property
195 def pixel_to_sky_transform(self) -> Transform[F, SkyFrame]:
196 """Low-level transform from pixel to sky coordinates (`Transform`)."""
197 return self._pixel_to_sky
199 @property
200 def sky_to_pixel_transform(self) -> Transform[SkyFrame, F]:
201 """Low-level transform from sky to pixel coordinates (`Transform`)."""
202 return self._pixel_to_sky.inverted()
204 @property
205 def fits_approximation(self) -> SkyProjection[F] | None:
206 """An approximation to this projection that is guaranteed to have an
207 `as_fits_wcs` method that does not return `None`.
208 """
209 return SkyProjection(self._fits_approximation) if self._fits_approximation is not None else None
211 def show(self, simplified: bool = False, comments: bool = False) -> str:
212 """Return the AST native representation of the transform.
214 Parameters
215 ----------
216 simplified
217 Whether to ask AST to simplify the mapping before showing it.
218 This will make it much more likely that two equivalent transforms
219 have the same `show` result.
220 comments
221 Whether to include descriptive comments.
222 """
223 return self._pixel_to_sky.show(simplified=simplified, comments=comments)
225 def pixel_to_sky[T: np.ndarray | float](self, *, x: T, y: T) -> SkyCoord:
226 """Transform one or more pixel points to sky coordinates.
228 Parameters
229 ----------
230 x : `numpy.ndarray` | `float`
231 ``x`` values of the pixel points to transform.
232 y : `numpy.ndarray` | `float`
233 ``y`` values of the pixel points to transform.
235 Returns
236 -------
237 astropy.coordinates.SkyCoord
238 Transformed sky coordinates.
239 """
240 sky_rad = self._pixel_to_sky.apply_forward(x=x, y=y)
241 return SkyCoord(ra=sky_rad.x, dec=sky_rad.y, unit=u.rad)
243 def sky_to_pixel(self, sky: SkyCoord) -> XY[np.ndarray | float]:
244 """Transform one or more sky coordinates to pixels.
246 Parameters
247 ----------
248 sky
249 Sky coordinates to transform.
251 Returns
252 -------
253 `XY` [`numpy.ndarray` | `float`]
254 Transformed pixel coordinates.
255 """
256 if sky.frame.name != "icrs":
257 sky = sky.transform_to("icrs")
258 ra: Longitude = sky.ra
259 dec: Latitude = sky.dec
260 return self._pixel_to_sky.apply_inverse(
261 x=ra.to_value(u.rad),
262 y=dec.to_value(u.rad),
263 )
265 def as_astropy(self, bbox: Box | None = None) -> SkyProjectionAstropyView:
266 """Return an `astropy.wcs` view of this `SkyProjection`.
268 Parameters
269 ----------
270 bbox
271 Bounding box of the array the view will describe. This
272 projection object is assumed to work on the same coordinate system
273 in which ``bbox`` is defined, while the Astropy view will consider
274 the first row and column in that box to be ``(0, 0)``.
276 Notes
277 -----
278 This returns an object that satisfies the
279 `astropy.wcs.wcsapi.BaseHighLevelWCS` and
280 `astropy.wcs.wcsapi.BaseLowLevelWCS` interfaces while evaluating the
281 underlying `SkyProjection` itself. It is *not* an `astropy.wcs.WCS`
282 instance, which is a type that also satisfies those interfaces but
283 only supports FITS WCS representations (see `as_fits_wcs`).
284 """
285 return SkyProjectionAstropyView(self._pixel_to_sky._ast_mapping, bbox)
287 def as_fits_wcs(self, bbox: Box, allow_approximation: bool = False) -> astropy.wcs.WCS | None:
288 """Return a FITS WCS representation of this projection, if possible.
290 Parameters
291 ----------
292 bbox
293 Bounding box of the array the FITS WCS will describe. This
294 transform object is assumed to work on the same coordinate system
295 in which ``bbox`` is defined, while the FITS WCS will consider the
296 first row and column in that box to be ``(0, 0)`` (in Astropy
297 interfaces) or ``(1, 1)`` (in the FITS representation itself).
298 allow_approximation
299 If `True` and this `SkyProjection` holds a FITS approximation to
300 itself, return that approximation.
301 """
302 if allow_approximation and self._fits_approximation: 302 ↛ 303line 302 didn't jump to line 303 because the condition on line 302 was never true
303 return self._fits_approximation.as_fits_wcs(bbox)
304 return self._pixel_to_sky.as_fits_wcs(bbox)
306 def serialize[P: pydantic.BaseModel](
307 self, archive: OutputArchive[P], *, use_frame_sets: bool = False
308 ) -> SkyProjectionSerializationModel[P]:
309 """Serialize a projection to an archive.
311 Parameters
312 ----------
313 archive
314 Archive to serialize to.
315 use_frame_sets
316 If `True`, decompose the underlying transform and try to reference
317 component mappings that were already serialized into a `FrameSet`
318 in the archive. The FITS approximation transform is never
319 decomposed.
321 Returns
322 -------
323 `SkyProjectionSerializationModel`
324 Serialized form of the projection.
325 """
326 pixel_to_sky = archive.serialize_direct(
327 "pixel_to_sky", functools.partial(self._pixel_to_sky.serialize, use_frame_sets=use_frame_sets)
328 )
329 fits_approximation = (
330 archive.serialize_direct("fits_approximation", self._fits_approximation.serialize)
331 if self._fits_approximation is not None
332 else None
333 )
334 return SkyProjectionSerializationModel(
335 pixel_to_sky=pixel_to_sky, fits_approximation=fits_approximation
336 )
338 @staticmethod
339 def _get_archive_tree_type[P: pydantic.BaseModel](
340 pointer_type: type[P],
341 ) -> type[SkyProjectionSerializationModel[P]]:
342 """Return the serialization model type for this object for an archive
343 type that uses the given pointer type.
344 """
345 return SkyProjectionSerializationModel[pointer_type] # type: ignore
347 @staticmethod
348 def from_legacy(
349 sky_wcs: LegacySkyWcs, pixel_frame: F, pixel_bounds: Bounds | None = None
350 ) -> SkyProjection[F]:
351 """Construct a transform from a legacy `lsst.afw.geom.SkyWcs`.
353 Parameters
354 ----------
355 sky_wcs : `lsst.afw.geom.SkyWcs`
356 Legacy WCS object.
357 pixel_frame
358 Coordinate frame for the pixel grid.
359 pixel_bounds
360 The region that bounds valid pixels for this transform.
361 """
362 fits_approximation: Transform[F, SkyFrame] | None = None
363 if (legacy_fits_approximation := sky_wcs.getFitsApproximation()) is not None:
364 fits_approximation = Transform(
365 pixel_frame,
366 SkyFrame.ICRS,
367 legacy_fits_approximation.getFrameDict(),
368 pixel_bounds,
369 )
370 return SkyProjection(
371 Transform(pixel_frame, SkyFrame.ICRS, sky_wcs.getFrameDict(), pixel_bounds),
372 fits_approximation=fits_approximation,
373 )
375 def to_legacy(self) -> LegacySkyWcs:
376 """Convert to a legacy `lsst.afw.geom.SkyWcs` instance."""
377 from lsst.afw.geom import SkyWcs as LegacySkyWcs
379 try:
380 ast_mapping = astshim.FrameDict(self._pixel_to_sky._ast_mapping)
381 except TypeError as err:
382 err.add_note(
383 "Only Projections created by from_legacy and from_fits_wcs "
384 "are guaranteed to be convertible to SkyWcs."
385 )
386 raise
387 legacy_wcs = LegacySkyWcs(ast_mapping)
388 if self.fits_approximation is not None:
389 legacy_wcs = legacy_wcs.copyWithFitsApproximation(self.fits_approximation.to_legacy())
390 return legacy_wcs
393class SkyProjectionAstropyView(BaseLowLevelWCS, HighLevelWCSMixin):
394 """An Astropy-interface view of a `SkyProjection`.
396 Parameters
397 ----------
398 ast_pixel_to_sky
399 AST mapping from pixel coordinates to sky coordinates.
400 bbox
401 Bounding box of the projection, or `None` if unbounded.
403 Notes
404 -----
405 The constructor of this classe is considered a private implementation
406 detail; use `SkyProjection.as_astropy` instead.
408 This object satisfies the `astropy.wcs.wcsapi.BaseHighLevelWCS` and
409 `astropy.wcs.wcsapi.BaseLowLevelWCS` interfaces while evaluating the
410 underlying `SkyProjection` itself. It is *not* an `astropy.wcs.WCS`
411 subclass, which is a type that also satisfies those interfaces but
412 only supports FITS WCS representations (see `SkyProjection.as_fits_wcs`).
413 """
415 def __init__(self, ast_pixel_to_sky: astshim.Mapping, bbox: Box | None) -> None:
416 self._bbox = bbox
417 if bbox is not None: 417 ↛ 419line 417 didn't jump to line 419 because the condition on line 417 was always true
418 ast_pixel_to_sky = astshim.ShiftMap(list(bbox.start.xy)).then(ast_pixel_to_sky)
419 self._ast_pixel_to_sky = ast_pixel_to_sky
421 @property
422 def low_level_wcs(self) -> Self:
423 return self
425 @property
426 def array_shape(self) -> YX[int] | None:
427 return self._bbox.shape if self._bbox is not None else None
429 @property
430 def axis_correlation_matrix(self) -> np.ndarray:
431 return np.array([[True, True], [True, True]])
433 @property
434 def pixel_axis_names(self) -> XY[str]:
435 return XY("x", "y")
437 @property
438 def pixel_bounds(self) -> XY[tuple[int, int]] | None:
439 if self._bbox is None:
440 return None
441 return XY((self._bbox.x.min, self._bbox.x.max), (self._bbox.y.min, self._bbox.y.max))
443 @property
444 def pixel_n_dim(self) -> int:
445 return 2
447 @property
448 def pixel_shape(self) -> XY[int] | None:
449 array_shape = self.array_shape
450 return array_shape.xy if array_shape is not None else None
452 @property
453 def serialized_classes(self) -> bool:
454 return False
456 @property
457 def world_axis_names(self) -> tuple[str, str]:
458 return ("ra", "dec")
460 @property
461 def world_axis_object_classes(self) -> dict[str, tuple[type[SkyCoord], tuple[()], dict[str, Any]]]:
462 return {"celestial": (SkyCoord, (), {"frame": ICRS, "unit": (u.rad, u.rad)})}
464 @property
465 def world_axis_object_components(self) -> list[tuple[str, int, str]]:
466 return [("celestial", 0, "spherical.lon.radian"), ("celestial", 1, "spherical.lat.radian")]
468 @property
469 def world_axis_physical_types(self) -> tuple[str, str]:
470 return ("pos.eq.ra", "pos.eq.dec")
472 @property
473 def world_axis_units(self) -> tuple[str, str]:
474 return ("rad", "rad")
476 @property
477 def world_n_dim(self) -> int:
478 return 2
480 def pixel_to_world_values(self, x: np.ndarray, y: np.ndarray) -> XY[np.ndarray]:
481 return _ast_apply(self._ast_pixel_to_sky.applyForward, x=x, y=y)
483 def world_to_pixel_values(self, ra: np.ndarray, dec: np.ndarray) -> XY[np.ndarray]:
484 return _ast_apply(self._ast_pixel_to_sky.applyInverse, x=ra, y=dec)
487class SkyProjectionSerializationModel[P: pydantic.BaseModel](ArchiveTree):
488 """Serialization model for projetions."""
490 SCHEMA_NAME: ClassVar[str] = "sky_projection"
491 SCHEMA_VERSION: ClassVar[str] = "1.0.0"
492 MIN_READ_VERSION: ClassVar[int] = 1
493 PUBLIC_TYPE: ClassVar[type] = SkyProjection
495 pixel_to_sky: TransformSerializationModel[P] = pydantic.Field(
496 description="The transform that maps pixel coordinates to the sky."
497 )
498 fits_approximation: TransformSerializationModel[P] | None = pydantic.Field(
499 default=None,
500 description=(
501 "An approximation of the pixel-to-sky transform that is exactly representable as a FITS WCS."
502 ),
503 exclude_if=is_none,
504 )
506 def deserialize(self, archive: InputArchive[P], **kwargs: Any) -> SkyProjection[Any]:
507 """Deserialize a projection from an archive.
509 Parameters
510 ----------
511 archive
512 Archive to read from.
513 **kwargs
514 Unsupported keyword arguments are accepted only to provide better
515 error messages (raising `serialization.InvalidParameterError`).
516 """
517 if kwargs: 517 ↛ 518line 517 didn't jump to line 518 because the condition on line 517 was never true
518 raise InvalidParameterError(f"Unrecognized parameters for Projection: {set(kwargs.keys())}.")
519 pixel_to_sky = self.pixel_to_sky.deserialize(archive)
520 fits_approximation = (
521 self.fits_approximation.deserialize(archive) if self.fits_approximation is not None else None
522 )
523 return SkyProjection(pixel_to_sky, fits_approximation=fits_approximation)