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