Coverage for python/lsst/images/_image.py: 82%
228 statements
« prev ^ index » next coverage.py v7.14.3, created at 2026-06-29 22:39 +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__ = ("Image", "ImageSerializationModel")
16from collections.abc import Callable, Sequence
17from contextlib import ExitStack
18from types import EllipsisType
19from typing import TYPE_CHECKING, Any, ClassVar, final
21import astropy.io.fits
22import astropy.units
23import astropy.wcs
24import numpy as np
25import numpy.typing as npt
26import pydantic
28from lsst.resources import ResourcePath, ResourcePathExpression
30from . import fits
31from ._generalized_image import GeneralizedImage
32from ._geom import YX, Box
33from ._transforms import Frame, GeneralFrame, SkyProjection, SkyProjectionSerializationModel
34from .serialization import (
35 ArchiveTree,
36 ArrayReferenceModel,
37 ArrayReferenceQuantityModel,
38 InlineArrayModel,
39 InlineArrayQuantityModel,
40 InputArchive,
41 InvalidParameterError,
42 MetadataValue,
43 OutputArchive,
44 no_header_updates,
45)
46from .utils import is_none
48if TYPE_CHECKING:
49 try:
50 from lsst.afw.image import Image as LegacyImage
51 except ImportError:
52 type LegacyImage = Any # type: ignore[no-redef]
55DEFAULT_PIXEL_FRAME = GeneralFrame(unit=astropy.units.pix)
56"""The pixel-grid `Frame` assumed when reconstructing a `SkyProjection` from a
57FITS header that does not otherwise identify its pixel frame, consistent with
58the FITS standard's notion of a plain pixel axis.
59"""
62@final
63class Image(GeneralizedImage):
64 """A 2-d array that may be augmented with units and a nonzero origin.
66 Parameters
67 ----------
68 array_or_fill
69 Array or fill value for the image. If a fill value, ``bbox`` or
70 ``shape`` must be provided.
71 bbox
72 Bounding box for the image.
73 yx0
74 Logical coordinates of the first pixel in the array, ordered ``y``,
75 ``x`` (unless an `XY` instance is passed). Ignored if
76 ``bbox`` is provided. Defaults to zeros.
77 shape
78 Leading dimensions of the array, ordered ``y``, ``x`` (unless an `XY`
79 instance is passed). Only needed if ``array_or_fill`` is not an
80 array and ``bbox`` is not provided. Like the bbox, this does not
81 include the last dimension of the array.
82 dtype
83 Pixel data type override.
84 unit
85 Units for the image's pixel values.
86 sky_projection
87 Projection that maps the pixel grid to the sky.
88 metadata
89 Arbitrary flexible metadata to associate with the image.
91 Notes
92 -----
93 Indexing the `array` attribute of an `Image` does not take into account its
94 ``yx0`` offset, but accessing a subimage by indexing an `Image` with a
95 `Box` does, and the `bbox` of the subimage is set to match its location
96 within the original image.
98 Indexed assignment to a subimage requires consistency between the
99 coordinate systems and units of both operands, but it will automatically
100 select a subimage of the right-hand side and convert compatible units when
101 possible. In other words::
103 a[box] = b
105 is a shortcut for
107 a[box].quantity = b[box].quantity
109 An ellipsis (``...``) can be used instead of a `Box` to assign to the full
110 image.
111 """
113 def __init__(
114 self,
115 array_or_fill: np.ndarray | int | float = 0,
116 /,
117 *,
118 bbox: Box | None = None,
119 yx0: Sequence[int] | None = None,
120 shape: Sequence[int] | None = None,
121 dtype: npt.DTypeLike | None = None,
122 unit: astropy.units.UnitBase | None = None,
123 sky_projection: SkyProjection[Any] | None = None,
124 metadata: dict[str, MetadataValue] | None = None,
125 ) -> None:
126 super().__init__(metadata)
127 if isinstance(array_or_fill, np.ndarray):
128 if dtype is not None:
129 array = np.array(array_or_fill, dtype=dtype, copy=None)
130 else:
131 array = array_or_fill
132 if bbox is None:
133 bbox = Box.from_shape(array.shape, start=yx0)
134 elif bbox.shape != array.shape:
135 raise ValueError(
136 f"Explicit bbox shape {bbox.shape} does not match array with shape {array.shape}."
137 )
138 if shape is not None and shape != array.shape:
139 raise ValueError(f"Explicit shape {shape} does not match array with shape {array.shape}.")
140 else:
141 if bbox is None:
142 if shape is None:
143 raise TypeError("No bbox, shape, or array provided.")
144 bbox = Box.from_shape(shape, start=yx0)
145 elif shape is not None and shape != bbox.shape:
146 raise ValueError(f"Explicit shape {shape} does not match bbox shape {bbox.shape}.")
147 array = np.full(bbox.shape, array_or_fill, dtype=dtype)
148 self._array: np.ndarray = array
149 self._bbox: Box = bbox
150 self._unit = unit
151 self._sky_projection = sky_projection
153 @property
154 def array(self) -> np.ndarray:
155 """The low-level array (`numpy.ndarray`).
157 Assigning to this attribute modifies the existing array in place; the
158 bounding box and underlying data pointer are never changed.
159 """
160 return self._array
162 @array.setter
163 def array(self, value: np.ndarray | int | float) -> None:
164 self._array[...] = value
166 @property
167 def quantity(self) -> astropy.units.Quantity:
168 """The low-level array with units (`astropy.units.Quantity`).
170 Assigning to this attribute modifies the existing array in place; the
171 bounding box and underlying data pointer are never changed.
172 """
173 return astropy.units.Quantity(self._array, self._unit, copy=False)
175 @quantity.setter
176 def quantity(self, value: astropy.units.Quantity) -> None:
177 self.quantity[...] = value
179 @property
180 def bbox(self) -> Box:
181 """Bounding box for the image (`Box`)."""
182 return self._bbox
184 @property
185 def unit(self) -> astropy.units.UnitBase | None:
186 """Units for the image's pixel values (`astropy.units.Unit` or
187 `None`).
188 """
189 return self._unit
191 @property
192 def sky_projection(self) -> SkyProjection[Any] | None:
193 """The projection that maps this image's pixel grid to the sky
194 (`SkyProjection` | `None`).
196 Notes
197 -----
198 The pixel coordinates used by this projection account for the bounding
199 box ``start``; they are not just array indices.
200 """
201 return self._sky_projection
203 def __getitem__(self, bbox: Box | EllipsisType) -> Image:
204 if bbox is ...:
205 return self
206 super().__getitem__(bbox)
207 indices = bbox.slice_within(self._bbox)
208 return self._transfer_metadata(
209 Image(self._array[indices], bbox=bbox, unit=self._unit, sky_projection=self._sky_projection),
210 bbox=bbox,
211 )
213 def __setitem__(self, bbox: Box | EllipsisType, value: Image) -> None:
214 self[bbox].quantity[...] = value.quantity
216 def __str__(self) -> str:
217 return f"Image({self.bbox!s}, {self.array.dtype.type.__name__})"
219 def __repr__(self) -> str:
220 return f"Image(..., bbox={self.bbox!r}, dtype={self.array.dtype!r})"
222 def __eq__(self, other: object) -> bool:
223 if not isinstance(other, Image):
224 return NotImplemented
225 return (
226 self._bbox == other._bbox
227 and self._unit == other._unit
228 and np.array_equal(self._array, other._array, equal_nan=True)
229 )
231 def copy(self) -> Image:
232 return self._transfer_metadata(
233 Image(self._array.copy(), bbox=self._bbox, unit=self._unit, sky_projection=self._sky_projection),
234 copy=True,
235 )
237 def view(
238 self,
239 *,
240 unit: astropy.units.UnitBase | None | EllipsisType = ...,
241 sky_projection: SkyProjection | None | EllipsisType = ...,
242 yx0: Sequence[int] | EllipsisType = ...,
243 ) -> Image:
244 """Make a view of the image, with optional updates.
246 Parameters
247 ----------
248 unit
249 Units for the view's pixel values. Defaults to the units of this
250 image.
251 sky_projection
252 Projection that maps the pixel grid to the sky. Defaults to the
253 projection of this image.
254 yx0
255 Logical coordinates of the first pixel, ordered ``y``, ``x``.
256 Defaults to the ``start`` of this image's bounding box.
257 """
258 if unit is ...: 258 ↛ 260line 258 didn't jump to line 260 because the condition on line 258 was always true
259 unit = self._unit
260 if sky_projection is ...: 260 ↛ 261line 260 didn't jump to line 261 because the condition on line 260 was never true
261 sky_projection = self._sky_projection
262 if yx0 is ...: 262 ↛ 264line 262 didn't jump to line 264 because the condition on line 262 was always true
263 yx0 = self._bbox.start
264 return self._transfer_metadata(Image(self._array, yx0=yx0, unit=unit, sky_projection=sky_projection))
266 def serialize[P: pydantic.BaseModel](
267 self,
268 archive: OutputArchive[P],
269 *,
270 update_header: Callable[[astropy.io.fits.Header], None] = no_header_updates,
271 save_projection: bool = True,
272 add_offset_wcs: str | None = "A",
273 tile_shape: tuple[int, ...] | None = None,
274 options_name: str | None = None,
275 ) -> ImageSerializationModel[P]:
276 """Serialize the image to an output archive.
278 Parameters
279 ----------
280 archive
281 Archive to write to.
282 update_header
283 A callback that will be given the FITS header for the HDU
284 containing this image in order to add keys to it. This callback
285 may be provided but will not be called if the output format is not
286 FITS.
287 save_projection
288 If `True`, save the `SkyProjection` attached to the image, if there
289 is one. This does not affect whether a FITS WCS corresponding to
290 the projection is written (it always is, if available, and if
291 ``add_offset_wcs`` is not ``" "``).
292 add_offset_wcs
293 A FITS WCS single-character suffix to use when adding a linear
294 WCS that maps the FITS array to the logical pixel coordinates
295 defined by ``bbox.start``. Set to `None` to not write this WCS.
296 If this is set to ``" "``, it will prevent the `SkyProjection` from
297 being saved as a FITS WCS.
298 tile_shape
299 The recommended shape of each tile, if the archive will save
300 the array in distinct tiles for faster subarray retrieval.
301 This is a hint; archives are not required to use this value.
302 options_name
303 Use this name to look up archive options.
304 """
306 def _update_header(header: astropy.io.fits.Header) -> None:
307 update_header(header)
308 if self.unit is not None:
309 try:
310 header["BUNIT"] = self.unit.to_string(format="fits")
311 except ValueError:
312 # Units not supported by FITS; write it anyway because
313 # the accepted units are just a recommendation in the
314 # standard.
315 header["BUNIT"] = self.unit.to_string()
316 if self.sky_projection is not None and add_offset_wcs != " ":
317 if self.fits_wcs:
318 header.update(self.fits_wcs.to_header(relax=True))
319 if add_offset_wcs is not None: 319 ↛ exitline 319 didn't return from function '_update_header' because the condition on line 319 was always true
320 fits.add_offset_wcs(header, x=self.bbox.x.start, y=self.bbox.y.start, key=add_offset_wcs)
322 array_model = archive.add_array(
323 self.array, update_header=_update_header, tile_shape=tile_shape, options_name=options_name
324 )
325 serialized_projection: SkyProjectionSerializationModel[P] | None = None
326 if save_projection and self.sky_projection is not None:
327 serialized_projection = archive.serialize_direct("sky_projection", self.sky_projection.serialize)
328 data = array_model if self.unit is None else array_model.with_units(self.unit)
329 return ImageSerializationModel.model_construct(
330 data=data,
331 yx0=list(self.bbox.start),
332 sky_projection=serialized_projection,
333 metadata=self.metadata,
334 )
336 @staticmethod
337 def _get_archive_tree_type[P: pydantic.BaseModel](
338 pointer_type: type[P],
339 ) -> type[ImageSerializationModel[P]]:
340 """Return the serialization model type for this object for an archive
341 type that uses the given pointer type.
342 """
343 return ImageSerializationModel[pointer_type] # type: ignore
345 _archive_default_name: ClassVar[str] = "image"
346 """The name this object should be serialized with when written as the
347 top-level object.
348 """
350 @staticmethod
351 def from_legacy(legacy: LegacyImage, unit: astropy.units.UnitBase | None = None) -> Image:
352 """Convert from an `lsst.afw.image.Image` instance.
354 Parameters
355 ----------
356 legacy
357 An `lsst.afw.image.Image` instance that will share pixel data with
358 the returned object.
359 unit
360 Units of the image.
361 """
362 return Image(legacy.array, yx0=YX(y=legacy.getY0(), x=legacy.getX0()), unit=unit)
364 def to_legacy(self, *, copy: bool | None = None) -> LegacyImage:
365 """Convert to an `lsst.afw.image.Image` instance.
367 Parameters
368 ----------
369 copy
370 If `True`, always copy the pixel data. If `False`, return a view,
371 and raise `TypeError` if the pixel data is read-only (this is not
372 supported by afw). If `None`, onyl if the pixel data is
373 read-only.
374 """
375 import lsst.afw.image
376 import lsst.geom
378 array = self._array
379 if copy: 379 ↛ 380line 379 didn't jump to line 380 because the condition on line 379 was never true
380 array = array.copy()
381 elif not self._array.flags.writeable: 381 ↛ 382line 381 didn't jump to line 382 because the condition on line 381 was never true
382 if copy is None:
383 array = array.copy()
384 else:
385 raise TypeError("Cannot create a legacy lsst.afw.image.Image view into a read-only array.")
387 return lsst.afw.image.Image(
388 array,
389 deep=False,
390 dtype=array.dtype.type,
391 xy0=lsst.geom.Point2I(self._bbox.x.min, self._bbox.y.min),
392 )
394 @classmethod
395 def from_hdu_list(
396 cls,
397 hdu_list: astropy.io.fits.HDUList,
398 *,
399 fits_wcs_frame: Frame | None = DEFAULT_PIXEL_FRAME,
400 ) -> Image:
401 """Reconstruct an `~lsst.images.Image` from a cut-down ``lsst.images``
402 HDU list.
404 This reads only the first two HDUs (the primary HDU and the image
405 HDU), as written for the image-only cut-outs produced by
406 ``dax_images_cutout``: a real ``lsst.images`` file with its JSON-tree,
407 index, and any nested-archive HDUs dropped.
409 Parameters
410 ----------
411 hdu_list
412 HDU list whose first HDU is the primary header and whose second
413 HDU holds the image pixels.
414 fits_wcs_frame
415 Pixel-grid `~lsst.images.Frame` for the
416 `~lsst.images.SkyProjection` reconstructed from the image HDU's
417 FITS WCS. Defaults to a plain pixel frame; pass `None` to skip
418 attaching a projection.
420 Returns
421 -------
422 `~lsst.images.Image`
423 The reconstructed image, ready to be re-serialized as a normal
424 ``lsst.images`` file.
426 Notes
427 -----
428 The headers of the consumed HDUs are modified in place (WCS and other
429 interpreted cards are stripped), as in `read_legacy`.
430 """
431 opaque_metadata = fits.FitsOpaqueMetadata()
432 opaque_metadata.add_cutdown_primary_header(hdu_list[0].header)
433 result = cls._read_legacy_hdu(
434 hdu_list[1], opaque_metadata, preserve_bintable=None, fits_wcs_frame=fits_wcs_frame
435 )
436 result._opaque_metadata = opaque_metadata
437 return result
439 @staticmethod
440 def read_legacy(
441 uri: ResourcePathExpression,
442 *,
443 preserve_quantization: bool = False,
444 ext: str | int = 1,
445 fits_wcs_frame: Frame | None = None,
446 ) -> Image:
447 """Read a FITS file written by `lsst.afw.image.Image.writeFits`.
449 Parameters
450 ----------
451 uri
452 URI or file name.
453 preserve_quantization
454 If `True`, ensure that writing the image back out again will
455 exactly preserve quantization-compressed pixel values. This causes
456 the arrays to be marked as read-only and stores the original binary
457 table data for those planes in memory. If the `Image` is copied,
458 the precompressed pixel values are not transferred to the copy.
459 ext
460 Name or index of the FITS HDU to read.
461 fits_wcs_frame
462 If not `None` and the HDU containing the image has a FITS WCS,
463 attach a `SkyProjection` to the returned image by converting that
464 WCS.
465 """
466 opaque_metadata = fits.FitsOpaqueMetadata()
467 with ExitStack() as exit_stack:
468 fs, fspath = ResourcePath(uri).to_fsspec()
469 stream = exit_stack.enter_context(fs.open(fspath))
470 hdu_list = exit_stack.enter_context(astropy.io.fits.open(stream))
471 opaque_metadata.extract_legacy_primary_header(hdu_list[0].header)
472 bintable_hdu: astropy.io.fits.BinTableHDU | None = None
473 if preserve_quantization:
474 bintable_stream = exit_stack.enter_context(fs.open(fspath))
475 bintable_hdu_list = exit_stack.enter_context(
476 astropy.io.fits.open(bintable_stream, disable_image_compression=True)
477 )
478 bintable_hdu = bintable_hdu_list[ext]
479 result = Image._read_legacy_hdu(
480 hdu_list[ext], opaque_metadata, preserve_bintable=bintable_hdu, fits_wcs_frame=fits_wcs_frame
481 )
482 result._opaque_metadata = opaque_metadata
483 return result
485 @staticmethod
486 def _read_legacy_hdu(
487 hdu: astropy.io.fits.ImageHDU | astropy.io.fits.CompImageHDU,
488 opaque_metadata: fits.FitsOpaqueMetadata,
489 *,
490 preserve_bintable: astropy.io.fits.BinTableHDU | None,
491 fits_wcs_frame: Frame | None = None,
492 ) -> Image:
493 unit: astropy.units.UnitBase | None = None
494 if (fits_unit := hdu.header.pop("BUNIT", None)) is not None:
495 try:
496 unit = astropy.units.Unit(fits_unit, format="fits")
497 except ValueError:
498 # Accept non-FITS units by assuming Astropy can still figure
499 # them out if we don't specify the format.
500 unit = astropy.units.Unit(fits_unit)
501 if opaque_metadata.get_instrumental_unit() == astropy.units.electron: 501 ↛ 503line 501 didn't jump to line 503 because the condition on line 501 was never true
502 # Fix incorrect BUNIT='adu' in LSST preliminary_visit_image.
503 if unit == astropy.units.adu:
504 unit = astropy.units.electron
505 if unit == astropy.units.adu**2:
506 unit = astropy.units.electron**2
507 yx0 = fits.read_yx0(hdu.header)
508 hdu.header.remove("LTV1", ignore_missing=True)
509 hdu.header.remove("LTV2", ignore_missing=True)
510 read_only: bool = False
511 if preserve_bintable is not None: 511 ↛ 512line 511 didn't jump to line 512 because the condition on line 511 was never true
512 opaque_metadata.precompressed[hdu.name] = fits.PrecompressedImage.from_bintable(preserve_bintable)
513 read_only = True
514 sky_projection: SkyProjection | None = None
515 if fits_wcs_frame is not None:
516 try:
517 fits_wcs = astropy.wcs.WCS(hdu.header)
518 except KeyError:
519 pass
520 else:
521 sky_projection = SkyProjection.from_fits_wcs(
522 fits_wcs, pixel_frame=fits_wcs_frame, x0=yx0.x, y0=yx0.y
523 )
524 image = Image(hdu.data, yx0=yx0, unit=unit, sky_projection=sky_projection)
525 if read_only: 525 ↛ 526line 525 didn't jump to line 526 because the condition on line 525 was never true
526 image._array.flags["WRITEABLE"] = False
527 fits.strip_wcs_cards(hdu.header)
528 hdu.header.strip()
529 hdu.header.remove("EXTTYPE", ignore_missing=True)
530 hdu.header.remove("INHERIT", ignore_missing=True)
531 hdu.header.remove("UZSCALE", ignore_missing=True)
532 opaque_metadata.add_header(hdu.header)
533 return image
536class ImageSerializationModel[P: pydantic.BaseModel](ArchiveTree):
537 """Pydantic model used to represent the serialized form of an `.Image`."""
539 SCHEMA_NAME: ClassVar[str] = "image"
540 SCHEMA_VERSION: ClassVar[str] = "1.0.0"
541 MIN_READ_VERSION: ClassVar[int] = 1
542 PUBLIC_TYPE: ClassVar[type] = Image
544 data: ArrayReferenceQuantityModel | ArrayReferenceModel | InlineArrayModel | InlineArrayQuantityModel = (
545 pydantic.Field(description="Reference to pixel data.")
546 )
547 yx0: list[int] = pydantic.Field(
548 description="Coordinate of the first pixels in the array, ordered (y, x)."
549 )
550 sky_projection: SkyProjectionSerializationModel[P] | None = pydantic.Field(
551 default=None,
552 exclude_if=is_none,
553 description="Projection that maps the logical pixel grid onto the sky.",
554 )
556 @property
557 def bbox(self) -> Box:
558 """The bounding box of the image."""
559 match self.data:
560 case ArrayReferenceQuantityModel() | InlineArrayQuantityModel():
561 shape = self.data.value.shape
562 case ArrayReferenceModel() | InlineArrayModel(): 562 ↛ 564line 562 didn't jump to line 564 because the pattern on line 562 always matched
563 shape = self.data.shape
564 return Box.from_shape(shape, self.yx0)
566 def deserialize(
567 self,
568 archive: InputArchive[Any],
569 *,
570 bbox: Box | None = None,
571 strip_header: Callable[[astropy.io.fits.Header], None] = no_header_updates,
572 **kwargs: Any,
573 ) -> Image:
574 """Deserialize an image from an input archive.
576 Parameters
577 ----------
578 archive
579 Archive to read from.
580 bbox
581 Bounding box of a subimage to read instead.
582 strip_header
583 A callable that strips out any FITS header cards added by the
584 ``update_header`` argument in the corresponding call to
585 `Image.serialize`.
586 **kwargs
587 Unsupported keyword arguments are accepted only to provide better
588 error messages (raising `serialization.InvalidParameterError`).
589 """
590 if kwargs: 590 ↛ 591line 590 didn't jump to line 591 because the condition on line 590 was never true
591 raise InvalidParameterError(f"Unrecognized parameters for Image: {set(kwargs.keys())}.")
592 array_model: ArrayReferenceModel | InlineArrayModel
593 unit: astropy.units.UnitBase | None = None
594 if isinstance(self.data, ArrayReferenceQuantityModel | InlineArrayQuantityModel):
595 array_model = self.data.value
596 unit = self.data.unit
597 else:
598 array_model = self.data
600 def _strip_header(header: astropy.io.fits.Header) -> None:
601 if unit is not None:
602 header.pop("BUNIT", None)
603 fits.strip_wcs_cards(header)
604 strip_header(header)
606 slices = bbox.slice_within(self.bbox) if bbox is not None else ...
607 array = archive.get_array(array_model, strip_header=_strip_header, slices=slices)
608 sky_projection = self.sky_projection.deserialize(archive) if self.sky_projection is not None else None
609 return Image(
610 array,
611 yx0=self.yx0 if bbox is None else bbox.start,
612 unit=unit,
613 sky_projection=sky_projection,
614 )._finish_deserialize(self)
616 def deserialize_component(self, component: str, archive: InputArchive[Any], **kwargs: Any) -> Any:
617 if kwargs:
618 raise InvalidParameterError(f"Unsupported parameters for Image components: {set(kwargs.keys())}.")
619 return super().deserialize_component(component, archive)