Coverage for python / lsst / images / _transforms / _projection.py: 49%
151 statements
« prev ^ index » next coverage.py v7.14.0, created at 2026-05-14 08:07 +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__ = ("Projection", "ProjectionAstropyView", "ProjectionSerializationModel")
16import functools
17from typing import TYPE_CHECKING, Any, 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, 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)
46@final
47class Projection[F: Frame]:
48 """A transform from pixel coordinates to sky coordinates.
50 Parameters
51 ----------
52 pixel_to_sky
53 A low-level transform that maps pixel coordinates to sky coordinates.
54 fits_approximation
55 An approximation to ``pixel_to_sky`` that is guaranteed to have a
56 `~Transform.as_fits_wcs` method that does not return `None`. This
57 should not be provided if ``pixel_to_sky`` is itself representable
58 as a FITS WCS.
60 Notes
61 -----
62 `Transform` is conceptually immutable (the internal AST Mapping should
63 never be modified in-place after construction), and hence does not need to
64 be copied when any object that holds it is copied.
65 """
67 def __init__(
68 self, pixel_to_sky: Transform[F, SkyFrame], fits_approximation: Transform[F, SkyFrame] | None = None
69 ):
70 self._pixel_to_sky = pixel_to_sky
71 if pixel_to_sky.in_frame.unit != u.pix:
72 raise ValueError("Transform is not a mapping from pixel coordinates.")
73 if pixel_to_sky.out_frame != SkyFrame.ICRS:
74 raise ValueError("Transform is not a mapping to ICRS.")
75 self._fits_approximation = fits_approximation
77 @staticmethod
78 def from_fits_wcs(
79 fits_wcs: astropy.wcs.WCS,
80 pixel_frame: F,
81 pixel_bounds: Bounds | None = None,
82 x0: int = 0,
83 y0: int = 0,
84 ) -> Projection[F]:
85 """Construct a transform from a FITS WCS.
87 Parameters
88 ----------
89 fits_wcs
90 FITS WCS to convert.
91 pixel_frame
92 Coordinate frame for the pixel grid.
93 pixel_bounds
94 The region that bounds valid pixels for this transform.
95 x0
96 Logical coordinate of the first column in the array this WCS
97 relates to world coordinates.
98 y0
99 Logical coordinate of the first column in the array this WCS
100 relates to world coordinates.
102 Notes
103 -----
104 The ``x0`` and ``y0`` parameters reflect the fact that for FITS, the
105 first row and column are always labeled ``(1, 1)``, while in Astropy
106 and most other Python libraries, they are ``(0, 0)``. The `types` in
107 this package (e.g. `Image`, `Mask`) allow them to be any pair of
108 integers.
110 See Also
111 --------
112 Transform.from_fits_wcs
113 """
114 return Projection(
115 Transform.from_fits_wcs(
116 fits_wcs, pixel_frame, SkyFrame.ICRS, in_bounds=pixel_bounds, x0=x0, y0=y0
117 )
118 )
120 @property
121 def pixel_frame(self) -> F:
122 """Coordinate frame for the pixel grid."""
123 return self._pixel_to_sky.in_frame
125 @property
126 def sky_frame(self) -> SkyFrame:
127 """Coordinate frame for the sky."""
128 return self._pixel_to_sky.out_frame
130 @property
131 def pixel_bounds(self) -> Bounds | None:
132 """The region that bounds valid pixel points (`Bounds` | `None`)."""
133 return self._pixel_to_sky.in_bounds
135 @property
136 def pixel_to_sky_transform(self) -> Transform[F, SkyFrame]:
137 """Low-level transform from pixel to sky coordinates (`Transform`)."""
138 return self._pixel_to_sky
140 @property
141 def sky_to_pixel_transform(self) -> Transform[SkyFrame, F]:
142 """Low-level transform from sky to pixel coordinates (`Transform`)."""
143 return self._pixel_to_sky.inverted()
145 @property
146 def fits_approximation(self) -> Projection[F] | None:
147 """An approximation to this projection that is guaranteed to have an
148 `as_fits_wcs` method that does not return `None`.
149 """
150 return Projection(self._fits_approximation) if self._fits_approximation is not None else None
152 def show(self, simplified: bool = False, comments: bool = False) -> str:
153 """Return the AST native representation of the transform.
155 Parameters
156 ----------
157 simplified
158 Whether to ask AST to simplify the mapping before showing it.
159 This will make it much more likely that two equivalent transforms
160 have the same `show` result.
161 comments
162 Whether to include descriptive comments.
163 """
164 return self._pixel_to_sky.show(simplified=simplified, comments=comments)
166 def pixel_to_sky[T: np.ndarray | float](self, *, x: T, y: T) -> SkyCoord:
167 """Transform one or more pixel points to sky coordinates.
169 Parameters
170 ----------
171 x : `numpy.ndarray` | `float`
172 ``x`` values of the pixel points to transform.
173 y : `numpy.ndarray` | `float`
174 ``y`` values of the pixel points to transform.
176 Returns
177 -------
178 astropy.coordinates.SkyCoord
179 Transformed sky coordinates.
180 """
181 sky_rad = self._pixel_to_sky.apply_forward(x=x, y=y)
182 return SkyCoord(ra=sky_rad.x, dec=sky_rad.y, unit=u.rad)
184 def sky_to_pixel(self, sky: SkyCoord) -> XY[np.ndarray | float]:
185 """Transform one or more sky coordinates to pixels
187 Parameters
188 ----------
189 sky
190 Sky coordinates to transform.
192 Returns
193 -------
194 `XY` [`numpy.ndarray` | `float`]
195 Transformed pixel coordinates.
196 """
197 if sky.frame.name != "icrs":
198 sky = sky.transform_to("icrs")
199 ra: Longitude = sky.ra
200 dec: Latitude = sky.dec
201 return self._pixel_to_sky.apply_inverse(
202 x=ra.to_value(u.rad),
203 y=dec.to_value(u.rad),
204 )
206 def as_astropy(self, bbox: Box | None = None) -> ProjectionAstropyView:
207 """Return an `astropy.wcs` view of this `Projection`.
209 Parameters
210 ----------
211 bbox
212 Bounding box of the array the view will describe. This
213 projection object is assumed to work on the same coordinate system
214 in which ``bbox`` is defined, while the Astropy view will consider
215 the first row and column in that box to be ``(0, 0)``.
217 Notes
218 -----
219 This returns an object that satisfies the
220 `astropy.wcs.wcsapi.BaseHighLevelWCS` and
221 `astropy.wcs.wcsapi.BaseLowLevelWCS` interfaces while evaluating the
222 underlying `Projection` itself. It is *not* an `astropy.wcs.WCS`
223 instance, which is a type that also satisfies those interfaces but
224 only supports FITS WCS representations (see `as_fits_wcs`).
225 """
226 return ProjectionAstropyView(self._pixel_to_sky._ast_mapping, bbox)
228 def as_fits_wcs(self, bbox: Box, allow_approximation: bool = False) -> astropy.wcs.WCS | None:
229 """Return a FITS WCS representation of this projection, if possible.
231 Parameters
232 ----------
233 bbox
234 Bounding box of the array the FITS WCS will describe. This
235 transform object is assumed to work on the same coordinate system
236 in which ``bbox`` is defined, while the FITS WCS will consider the
237 first row and column in that box to be ``(0, 0)`` (in Astropy
238 interfaces) or ``(1, 1)`` (in the FITS representation itself).
239 allow_approximation
240 If `True` and this `Projection` holds a FITS approximation to
241 itself, return that approximation.
242 """
243 if allow_approximation and self._fits_approximation:
244 return self._fits_approximation.as_fits_wcs(bbox)
245 return self._pixel_to_sky.as_fits_wcs(bbox)
247 def serialize[P: pydantic.BaseModel](
248 self, archive: OutputArchive[P], *, use_frame_sets: bool = False
249 ) -> ProjectionSerializationModel[P]:
250 """Serialize a projection to an archive.
252 Parameters
253 ----------
254 archive
255 Archive to serialize to.
256 use_frame_sets
257 If `True`, decompose the underlying transform and try to reference
258 component mappings that were already serialized into a `FrameSet`
259 in the archive. The FITS approximation transform is never
260 decomposed.
262 Returns
263 -------
264 `ProjectionSerializationModel`
265 Serialized form of the projection.
266 """
267 pixel_to_sky = archive.serialize_direct(
268 "pixel_to_sky", functools.partial(self._pixel_to_sky.serialize, use_frame_sets=use_frame_sets)
269 )
270 fits_approximation = (
271 archive.serialize_direct("fits_approximation", self._fits_approximation.serialize)
272 if self._fits_approximation is not None
273 else None
274 )
275 return ProjectionSerializationModel(pixel_to_sky=pixel_to_sky, fits_approximation=fits_approximation)
277 @staticmethod
278 def _get_archive_tree_type[P: pydantic.BaseModel](
279 pointer_type: type[P],
280 ) -> type[ProjectionSerializationModel[P]]:
281 """Return the serialization model type for this object for an archive
282 type that uses the given pointer type.
283 """
284 return ProjectionSerializationModel[pointer_type] # type: ignore
286 @staticmethod
287 def from_legacy(
288 sky_wcs: LegacySkyWcs, pixel_frame: F, pixel_bounds: Bounds | None = None
289 ) -> Projection[F]:
290 """Construct a transform from a legacy `lsst.afw.geom.SkyWcs`.
292 Parameters
293 ----------
294 sky_wcs : `lsst.afw.geom.SkyWcs`
295 Legacy WCS object.
296 pixel_frame
297 Coordinate frame for the pixel grid.
298 pixel_bounds
299 The region that bounds valid pixels for this transform.
300 """
301 fits_approximation: Transform[F, SkyFrame] | None = None
302 if (legacy_fits_approximation := sky_wcs.getFitsApproximation()) is not None:
303 fits_approximation = Transform(
304 pixel_frame,
305 SkyFrame.ICRS,
306 legacy_fits_approximation.getFrameDict(),
307 pixel_bounds,
308 )
309 return Projection(
310 Transform(pixel_frame, SkyFrame.ICRS, sky_wcs.getFrameDict(), pixel_bounds),
311 fits_approximation=fits_approximation,
312 )
314 def to_legacy(self) -> LegacySkyWcs:
315 """Convert to a legacy `lsst.afw.geom.SkyWcs` instance."""
316 from lsst.afw.geom import SkyWcs as LegacySkyWcs
318 try:
319 ast_mapping = astshim.FrameDict(self._pixel_to_sky._ast_mapping)
320 except TypeError as err:
321 err.add_note(
322 "Only Projections created by from_legacy and from_fits_wcs "
323 "are guaranteed to be convertible to SkyWcs."
324 )
325 raise
326 legacy_wcs = LegacySkyWcs(ast_mapping)
327 if self.fits_approximation is not None:
328 legacy_wcs = legacy_wcs.copyWithFitsApproximation(self.fits_approximation.to_legacy())
329 return legacy_wcs
332class ProjectionAstropyView(BaseLowLevelWCS, HighLevelWCSMixin):
333 """An Astropy-interface view of a `Projection`.
335 Notes
336 -----
337 The constructor of this classe is considered a private implementation
338 detail; use `Projection.as_astropy` instead.
340 This object satisfies the `astropy.wcs.wcsapi.BaseHighLevelWCS` and
341 `astropy.wcs.wcsapi.BaseLowLevelWCS` interfaces while evaluating the
342 underlying `Projection` itself. It is *not* an `astropy.wcs.WCS`
343 subclass, which is a type that also satisfies those interfaces but
344 only supports FITS WCS representations (see `Projection.as_fits_wcs`).
345 """
347 def __init__(self, ast_pixel_to_sky: astshim.Mapping, bbox: Box | None):
348 self._bbox = bbox
349 if bbox is not None:
350 ast_pixel_to_sky = astshim.ShiftMap(list(bbox.start.xy)).then(ast_pixel_to_sky)
351 self._ast_pixel_to_sky = ast_pixel_to_sky
353 @property
354 def low_level_wcs(self) -> Self:
355 return self
357 @property
358 def array_shape(self) -> YX[int] | None:
359 return self._bbox.shape if self._bbox is not None else None
361 @property
362 def axis_correlation_matrix(self) -> np.ndarray:
363 return np.array([[True, True], [True, True]])
365 @property
366 def pixel_axis_names(self) -> XY[str]:
367 return XY("x", "y")
369 @property
370 def pixel_bounds(self) -> XY[tuple[int, int]] | None:
371 if self._bbox is None:
372 return None
373 return XY((self._bbox.x.min, self._bbox.x.max), (self._bbox.y.min, self._bbox.y.max))
375 @property
376 def pixel_n_dim(self) -> int:
377 return 2
379 @property
380 def pixel_shape(self) -> XY[int] | None:
381 array_shape = self.array_shape
382 return array_shape.xy if array_shape is not None else None
384 @property
385 def serialized_classes(self) -> bool:
386 return False
388 @property
389 def world_axis_names(self) -> tuple[str, str]:
390 return ("ra", "dec")
392 @property
393 def world_axis_object_classes(self) -> dict[str, tuple[type[SkyCoord], tuple[()], dict[str, Any]]]:
394 return {"celestial": (SkyCoord, (), {"frame": ICRS, "unit": (u.rad, u.rad)})}
396 @property
397 def world_axis_object_components(self) -> list[tuple[str, int, str]]:
398 return [("celestial", 0, "spherical.lon.radian"), ("celestial", 1, "spherical.lat.radian")]
400 @property
401 def world_axis_physical_types(self) -> tuple[str, str]:
402 return ("pos.eq.ra", "pos.eq.dec")
404 @property
405 def world_axis_units(self) -> tuple[str, str]:
406 return ("rad", "rad")
408 @property
409 def world_n_dim(self) -> int:
410 return 2
412 def pixel_to_world_values(self, x: np.ndarray, y: np.ndarray) -> XY[np.ndarray]:
413 return _ast_apply(self._ast_pixel_to_sky.applyForward, x=x, y=y)
415 def world_to_pixel_values(self, ra: np.ndarray, dec: np.ndarray) -> XY[np.ndarray]:
416 return _ast_apply(self._ast_pixel_to_sky.applyInverse, x=ra, y=dec)
419class ProjectionSerializationModel[P: pydantic.BaseModel](ArchiveTree):
420 """Serialization model for projetions."""
422 pixel_to_sky: TransformSerializationModel[P] = pydantic.Field(
423 description="The transform that maps pixel coordinates to the sky."
424 )
425 fits_approximation: TransformSerializationModel[P] | None = pydantic.Field(
426 default=None,
427 description=(
428 "An approximation of the pixel-to-sky transform that is exactly representable as a FITS WCS."
429 ),
430 exclude_if=is_none,
431 )
433 def deserialize(self, archive: InputArchive[P]) -> Projection[Any]:
434 """Deserialize a projection from an archive.
436 Parameters
437 ----------
438 archive
439 Archive to read from.
440 """
441 pixel_to_sky = self.pixel_to_sky.deserialize(archive)
442 fits_approximation = (
443 self.fits_approximation.deserialize(archive) if self.fits_approximation is not None else None
444 )
445 return Projection(pixel_to_sky, fits_approximation=fits_approximation)