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