Coverage for python / lsst / images / _transforms / _projection.py: 52%
140 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-05 23:14 -0800
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 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 ._frames import Frame, SkyFrame
30from ._transform import Transform, TransformSerializationModel, _ast_apply
32# This pre-python-3.12 declaration is needed by Sphinx (probably the
33# autodoc-typehints plugin.
34F = TypeVar("F", bound=Frame)
37@final
38class Projection[F: Frame]:
39 """A transform from pixel coordinates to sky coordinates.
41 Parameters
42 ----------
43 pixel_to_sky
44 A low-level transform that maps pixel coordinates to sky coordinates.
45 fits_approximation
46 An approximation to ``pixel_to_sky`` that is guaranteed to have a
47 `~Transform.as_fits_wcs` method that does not return `None`. This
48 should not be provided if ``pixel_to_sky`` is itself representable
49 as a FITS WCS.
51 Notes
52 -----
53 `Transform` is conceptually immutable (the internal AST Mapping should
54 never be modified in-place after construction), and hence does not need to
55 be copied when any object that holds it is copied.
56 """
58 def __init__(
59 self, pixel_to_sky: Transform[F, SkyFrame], fits_approximation: Transform[F, SkyFrame] | None = None
60 ):
61 self._pixel_to_sky = pixel_to_sky
62 if pixel_to_sky.in_frame.unit != u.pix:
63 raise ValueError("Transform is not a mapping from pixel coordinates.")
64 if pixel_to_sky.out_frame != SkyFrame.ICRS:
65 raise ValueError("Transform is not a mapping to ICRS.")
66 self._fits_approximation = fits_approximation
68 @staticmethod
69 def from_fits_wcs(
70 fits_wcs: astropy.wcs.WCS,
71 pixel_frame: F,
72 pixel_bounds: Bounds | None = None,
73 x0: int = 0,
74 y0: int = 0,
75 ) -> Projection[F]:
76 """Construct a transform from a FITS WCS.
78 Parameters
79 ----------
80 fits_wcs
81 FITS WCS to convert.
82 pixel_frame
83 Coordinate frame for the pixel grid.
84 pixel_bounds
85 The region that bounds valid pixels for this transform.
86 x0
87 Logical coordinate of the first column in the array this WCS
88 relates to world coordinates.
89 y0
90 Logical coordinate of the first column in the array this WCS
91 relates to world coordinates.
93 Notes
94 -----
95 The ``x0`` and ``y0`` parameters reflect the fact that for FITS, the
96 first row and column are always labeled ``(1, 1)``, while in Astropy
97 and most other Python libraries, they are ``(0, 0)``. The `types` in
98 this package (e.g. `Image`, `Mask`) allow them to be any pair of
99 integers.
101 See Also
102 --------
103 Transform.from_fits_wcs
104 """
105 return Projection(
106 Transform.from_fits_wcs(
107 fits_wcs, pixel_frame, SkyFrame.ICRS, in_bounds=pixel_bounds, x0=x0, y0=y0
108 )
109 )
111 @property
112 def pixel_frame(self) -> F:
113 """Coordinate frame for the pixel grid."""
114 return self._pixel_to_sky.in_frame
116 @property
117 def sky_frame(self) -> SkyFrame:
118 """Coordinate frame for the sky."""
119 return self._pixel_to_sky.out_frame
121 @property
122 def pixel_bounds(self) -> Bounds | None:
123 """The region that bounds valid pixel points (`Bounds` | `None`)."""
124 return self._pixel_to_sky.in_bounds
126 @property
127 def pixel_to_sky_transform(self) -> Transform[F, SkyFrame]:
128 """Low-level transform from pixel to sky coordinates (`Transform`)."""
129 return self._pixel_to_sky
131 @property
132 def sky_to_pixel_transform(self) -> Transform[SkyFrame, F]:
133 """Low-level transform from sky to pixel coordinates (`Transform`)."""
134 return self._pixel_to_sky.inverted()
136 @property
137 def fits_approximation(self) -> Projection[F] | None:
138 """An approximation to this projection that is guaranteed to have an
139 `as_fits_wcs` method that does not return `None`.
140 """
141 return Projection(self._fits_approximation) if self._fits_approximation is not None else None
143 def show(self, simplified: bool = False, comments: bool = False) -> str:
144 """Return the AST native representation of the transform.
146 Parameters
147 ----------
148 simplified
149 Whether to ask AST to simplify the mapping before showing it.
150 This will make it much more likely that two equivalent transforms
151 have the same `show` result.
152 comments
153 Whether to include descriptive comments.
154 """
155 return self._pixel_to_sky.show(simplified=simplified, comments=comments)
157 def pixel_to_sky[T: np.ndarray | float](self, *, x: T, y: T) -> SkyCoord:
158 """Transform one or more pixel points to sky coordinates.
160 Parameters
161 ----------
162 x : `numpy.ndarray` | `float`
163 ``x`` values of the pixel points to transform.
164 y : `numpy.ndarray` | `float`
165 ``y`` values of the pixel points to transform.
167 Returns
168 -------
169 astropy.coordinates.SkyCoord
170 Transformed sky coordinates.
171 """
172 sky_rad = self._pixel_to_sky.apply_forward(x=x, y=y)
173 return SkyCoord(ra=sky_rad.x, dec=sky_rad.y, unit=u.rad)
175 def sky_to_pixel(self, sky: SkyCoord) -> XY[np.ndarray | float]:
176 """Transform one or more sky coordinates to pixels
178 Parameters
179 ----------
180 sky
181 Sky coordinates to transform.
183 Returns
184 -------
185 `XY` [`numpy.ndarray` | `float`]
186 Transformed pixel coordinates.
187 """
188 if sky.frame.name != "icrs":
189 sky = sky.transform_to("icrs")
190 ra: Longitude = sky.ra
191 dec: Latitude = sky.dec
192 return self._pixel_to_sky.apply_inverse(
193 x=ra.to_value(u.rad),
194 y=dec.to_value(u.rad),
195 )
197 def as_astropy(self, bbox: Box | None = None) -> ProjectionAstropyView:
198 """Return an `astropy.wcs` view of this `Projection`.
200 Parameters
201 ----------
202 bbox
203 Bounding box of the array the view will describe. This
204 projection object is assumed to work on the same coordinate system
205 in which ``bbox`` is defined, while the Astropy view will consider
206 the first row and column in that box to be ``(0, 0)``.
208 Notes
209 -----
210 This returns an object that satisfies the
211 `astropy.wcs.wcsapi.BaseHighLevelWCS` and
212 `astropy.wcs.wcsapi.BaseLowLevelWCS` interfaces while evaluating the
213 underlying `Projection` itself. It is *not* an `astropy.wcs.WCS`
214 instance, which is a type that also satisfies those interfaces but
215 only supports FITS WCS representations (see `as_fits_wcs`).
216 """
217 return ProjectionAstropyView(self._pixel_to_sky._ast_mapping, bbox)
219 def as_fits_wcs(self, bbox: Box, allow_approximation: bool = False) -> astropy.wcs.WCS | None:
220 """Return a FITS WCS representation of this projection, if possible.
222 Parameters
223 ----------
224 bbox
225 Bounding box of the array the FITS WCS will describe. This
226 transform object is assumed to work on the same coordinate system
227 in which ``bbox`` is defined, while the FITS WCS will consider the
228 first row and column in that box to be ``(0, 0)`` (in Astropy
229 interfaces) or ``(1, 1)`` (in the FITS representation itself).
230 allow_approximation
231 If `True` and this `Projection` holds a FITS approximation to
232 itself, return that approximation.
233 """
234 if allow_approximation and self._fits_approximation:
235 return self._fits_approximation.as_fits_wcs(bbox)
236 return self._pixel_to_sky.as_fits_wcs(bbox)
238 def serialize[P: pydantic.BaseModel](
239 self, archive: OutputArchive[P], *, use_frame_sets: bool = False
240 ) -> ProjectionSerializationModel[P]:
241 """Serialize a projection to an archive.
243 Parameters
244 ----------
245 archive
246 Archive to serialize to.
247 use_frame_sets
248 If `True`, decompose the underlying transform and try to reference
249 component mappings that were already serialized into a `FrameSet`
250 in the archive. The FITS approximation transform is never
251 decomposed.
253 Returns
254 -------
255 `ProjectionSerializationModel`
256 Serialized form of the projection.
257 """
258 pixel_to_sky = archive.serialize_direct(
259 "pixel_to_sky", functools.partial(self._pixel_to_sky.serialize, use_frame_sets=use_frame_sets)
260 )
261 fits_approximation = (
262 archive.serialize_direct("fits_approximation", self._fits_approximation.serialize)
263 if self._fits_approximation is not None
264 else None
265 )
266 return ProjectionSerializationModel(pixel_to_sky=pixel_to_sky, fits_approximation=fits_approximation)
268 @staticmethod
269 def deserialize[P: pydantic.BaseModel](
270 model: ProjectionSerializationModel[P], archive: InputArchive[P]
271 ) -> Projection[Any]:
272 """Deserialize a projection from an archive.
274 Parameters
275 ----------
276 model
277 Seralized form of the projection.
278 archive
279 Archive to read from.
280 """
281 pixel_to_sky = Transform.deserialize(model.pixel_to_sky, archive)
282 fits_approximation = (
283 Transform.deserialize(model.fits_approximation, archive)
284 if model.fits_approximation is not None
285 else None
286 )
287 return Projection(pixel_to_sky, fits_approximation=fits_approximation)
289 @staticmethod
290 def _get_archive_tree_type[P: pydantic.BaseModel](
291 pointer_type: type[P],
292 ) -> type[ProjectionSerializationModel[P]]:
293 """Return the serialization model type for this object for an archive
294 type that uses the given pointer type.
295 """
296 return ProjectionSerializationModel[pointer_type] # type: ignore
298 @staticmethod
299 def from_legacy(sky_wcs: Any, pixel_frame: F, pixel_bounds: Bounds | None = None) -> Projection[F]:
300 """Construct a transform from a legacy `lsst.afw.geom.SkyWcs`.
302 Parameters
303 ----------
304 sky_wcs : `lsst.afw.geom.SkyWcs`
305 Legacy WCS object.
306 pixel_frame
307 Coordinate frame for the pixel grid.
308 pixel_bounds
309 The region that bounds valid pixels for this transform.
310 """
311 fits_approximation: Transform[F, SkyFrame] | None = None
312 if (legacy_fits_approximation := sky_wcs.getFitsApproximation()) is not None:
313 fits_approximation = Transform.from_legacy(
314 legacy_fits_approximation.getTransform(), pixel_frame, SkyFrame.ICRS, in_bounds=pixel_bounds
315 )
316 return Projection(
317 Transform.from_legacy(sky_wcs.getTransform(), pixel_frame, SkyFrame.ICRS, in_bounds=pixel_bounds),
318 fits_approximation=fits_approximation,
319 )
322class ProjectionAstropyView(BaseLowLevelWCS, HighLevelWCSMixin):
323 """An Astropy-interface view of a `Projection`.
325 Notes
326 -----
327 The constructor of this classe is considered a private implementation
328 detail; use `Projection.as_astropy` instead.
330 This object satisfies the `astropy.wcs.wcsapi.BaseHighLevelWCS` and
331 `astropy.wcs.wcsapi.BaseLowLevelWCS` interfaces while evaluating the
332 underlying `Projection` itself. It is *not* an `astropy.wcs.WCS`
333 subclass, which is a type that also satisfies those interfaces but
334 only supports FITS WCS representations (see `Projection.as_fits_wcs`).
335 """
337 def __init__(self, ast_pixel_to_sky: Any, bbox: Box | None):
338 self._bbox = bbox
339 if bbox is not None:
340 import astshim
342 ast_pixel_to_sky = astshim.ShiftMap(list(bbox.start.xy)).then(ast_pixel_to_sky)
343 self._ast_pixel_to_sky = ast_pixel_to_sky
345 @property
346 def low_level_wcs(self) -> Self:
347 return self
349 @property
350 def array_shape(self) -> YX[int] | None:
351 return self._bbox.shape if self._bbox is not None else None
353 @property
354 def axis_correlation_matrix(self) -> np.ndarray:
355 return np.array([[True, True], [True, True]])
357 @property
358 def pixel_axis_names(self) -> XY[str]:
359 return XY("x", "y")
361 @property
362 def pixel_bounds(self) -> XY[tuple[int, int]] | None:
363 if self._bbox is None:
364 return None
365 return XY((self._bbox.x.min, self._bbox.x.max), (self._bbox.y.min, self._bbox.y.max))
367 @property
368 def pixel_n_dim(self) -> int:
369 return 2
371 @property
372 def pixel_shape(self) -> XY[int] | None:
373 array_shape = self.array_shape
374 return array_shape.xy if array_shape is not None else None
376 @property
377 def serialized_classes(self) -> bool:
378 return False
380 @property
381 def world_axis_names(self) -> tuple[str, str]:
382 return ("ra", "dec")
384 @property
385 def world_axis_object_classes(self) -> dict[str, tuple[type[SkyCoord], tuple[()], dict[str, Any]]]:
386 return {"celestial": (SkyCoord, (), {"frame": ICRS, "unit": (u.rad, u.rad)})}
388 @property
389 def world_axis_object_components(self) -> list[tuple[str, int, str]]:
390 return [("celestial", 0, "spherical.lon.radian"), ("celestial", 1, "spherical.lat.radian")]
392 @property
393 def world_axis_physical_types(self) -> tuple[str, str]:
394 return ("pos.eq.ra", "pos.eq.dec")
396 @property
397 def world_axis_units(self) -> tuple[str, str]:
398 return ("rad", "rad")
400 @property
401 def world_n_dim(self) -> int:
402 return 2
404 def pixel_to_world_values(self, x: np.ndarray, y: np.ndarray) -> XY[np.ndarray]:
405 return _ast_apply(self._ast_pixel_to_sky.applyForward, x=x, y=y)
407 def world_to_pixel_values(self, ra: np.ndarray, dec: np.ndarray) -> XY[np.ndarray]:
408 return _ast_apply(self._ast_pixel_to_sky.applyInverse, x=ra, y=dec)
411class ProjectionSerializationModel[P: pydantic.BaseModel](ArchiveTree):
412 """Serialization model for projetions."""
414 pixel_to_sky: TransformSerializationModel[P] = pydantic.Field(
415 description="The transform that maps pixel coordinates to the sky."
416 )
417 fits_approximation: TransformSerializationModel[P] | None = pydantic.Field(
418 default=None,
419 description=(
420 "An approximation of the pixel-to-sky transform that is exactly representable as a FITS WCS."
421 ),
422 exclude_if=is_none,
423 )