Coverage for python/lsst/images/_image.py: 82%
225 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-16 15:24 -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__ = ("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 bbox, indices = self._handle_getitem_args(bbox)
205 return self._transfer_metadata(
206 Image(self._array[indices], bbox=bbox, unit=self._unit, sky_projection=self._sky_projection),
207 bbox=bbox,
208 )
210 def __setitem__(self, bbox: Box | EllipsisType, value: Image) -> None:
211 self[bbox].quantity[...] = value.quantity
213 def __str__(self) -> str:
214 return f"Image({self.bbox!s}, {self.array.dtype.type.__name__})"
216 def __repr__(self) -> str:
217 return f"Image(..., bbox={self.bbox!r}, dtype={self.array.dtype!r})"
219 def __eq__(self, other: object) -> bool:
220 if not isinstance(other, Image):
221 return NotImplemented
222 return (
223 self._bbox == other._bbox
224 and self._unit == other._unit
225 and np.array_equal(self._array, other._array, equal_nan=True)
226 )
228 def copy(self) -> Image:
229 return self._transfer_metadata(
230 Image(self._array.copy(), bbox=self._bbox, unit=self._unit, sky_projection=self._sky_projection),
231 copy=True,
232 )
234 def view(
235 self,
236 *,
237 unit: astropy.units.UnitBase | None | EllipsisType = ...,
238 sky_projection: SkyProjection | None | EllipsisType = ...,
239 yx0: Sequence[int] | EllipsisType = ...,
240 ) -> Image:
241 """Make a view of the image, with optional updates.
243 Parameters
244 ----------
245 unit
246 Units for the view's pixel values. Defaults to the units of this
247 image.
248 sky_projection
249 Projection that maps the pixel grid to the sky. Defaults to the
250 projection of this image.
251 yx0
252 Logical coordinates of the first pixel, ordered ``y``, ``x``.
253 Defaults to the ``start`` of this image's bounding box.
254 """
255 if unit is ...: 255 ↛ 257line 255 didn't jump to line 257 because the condition on line 255 was always true
256 unit = self._unit
257 if sky_projection is ...: 257 ↛ 258line 257 didn't jump to line 258 because the condition on line 257 was never true
258 sky_projection = self._sky_projection
259 if yx0 is ...: 259 ↛ 261line 259 didn't jump to line 261 because the condition on line 259 was always true
260 yx0 = self._bbox.start
261 return self._transfer_metadata(Image(self._array, yx0=yx0, unit=unit, sky_projection=sky_projection))
263 def serialize[P: pydantic.BaseModel](
264 self,
265 archive: OutputArchive[P],
266 *,
267 update_header: Callable[[astropy.io.fits.Header], None] = no_header_updates,
268 save_projection: bool = True,
269 add_offset_wcs: str | None = "A",
270 tile_shape: tuple[int, ...] | None = None,
271 options_name: str | None = None,
272 ) -> ImageSerializationModel[P]:
273 """Serialize the image to an output archive.
275 Parameters
276 ----------
277 archive
278 Archive to write to.
279 update_header
280 A callback that will be given the FITS header for the HDU
281 containing this image in order to add keys to it. This callback
282 may be provided but will not be called if the output format is not
283 FITS.
284 save_projection
285 If `True`, save the `SkyProjection` attached to the image, if there
286 is one. This does not affect whether a FITS WCS corresponding to
287 the projection is written (it always is, if available, and if
288 ``add_offset_wcs`` is not ``" "``).
289 add_offset_wcs
290 A FITS WCS single-character suffix to use when adding a linear
291 WCS that maps the FITS array to the logical pixel coordinates
292 defined by ``bbox.start``. Set to `None` to not write this WCS.
293 If this is set to ``" "``, it will prevent the `SkyProjection` from
294 being saved as a FITS WCS.
295 tile_shape
296 The recommended shape of each tile, if the archive will save
297 the array in distinct tiles for faster subarray retrieval.
298 This is a hint; archives are not required to use this value.
299 options_name
300 Use this name to look up archive options.
301 """
303 def _update_header(header: astropy.io.fits.Header) -> None:
304 update_header(header)
305 if self.unit is not None:
306 try:
307 header["BUNIT"] = self.unit.to_string(format="fits")
308 except ValueError:
309 # Units not supported by FITS; write it anyway because
310 # the accepted units are just a recommendation in the
311 # standard.
312 header["BUNIT"] = self.unit.to_string()
313 if self.sky_projection is not None and add_offset_wcs != " ":
314 if self.fits_wcs:
315 header.update(self.fits_wcs.to_header(relax=True))
316 if add_offset_wcs is not None: 316 ↛ exitline 316 didn't return from function '_update_header' because the condition on line 316 was always true
317 fits.add_offset_wcs(header, x=self.bbox.x.start, y=self.bbox.y.start, key=add_offset_wcs)
319 array_model = archive.add_array(
320 self.array, update_header=_update_header, tile_shape=tile_shape, options_name=options_name
321 )
322 serialized_projection: SkyProjectionSerializationModel[P] | None = None
323 if save_projection and self.sky_projection is not None:
324 serialized_projection = archive.serialize_direct("sky_projection", self.sky_projection.serialize)
325 data = array_model if self.unit is None else array_model.with_units(self.unit)
326 return ImageSerializationModel.model_construct(
327 data=data,
328 yx0=list(self.bbox.start),
329 sky_projection=serialized_projection,
330 metadata=self.metadata,
331 )
333 @staticmethod
334 def _get_archive_tree_type[P: pydantic.BaseModel](
335 pointer_type: type[P],
336 ) -> type[ImageSerializationModel[P]]:
337 """Return the serialization model type for this object for an archive
338 type that uses the given pointer type.
339 """
340 return ImageSerializationModel[pointer_type] # type: ignore
342 _archive_default_name: ClassVar[str] = "image"
343 """The name this object should be serialized with when written as the
344 top-level object.
345 """
347 @staticmethod
348 def from_legacy(legacy: LegacyImage, unit: astropy.units.UnitBase | None = None) -> Image:
349 """Convert from an `lsst.afw.image.Image` instance.
351 Parameters
352 ----------
353 legacy
354 An `lsst.afw.image.Image` instance that will share pixel data with
355 the returned object.
356 unit
357 Units of the image.
358 """
359 return Image(legacy.array, yx0=YX(y=legacy.getY0(), x=legacy.getX0()), unit=unit)
361 def to_legacy(self, *, copy: bool | None = None) -> LegacyImage:
362 """Convert to an `lsst.afw.image.Image` instance.
364 Parameters
365 ----------
366 copy
367 If `True`, always copy the pixel data. If `False`, return a view,
368 and raise `TypeError` if the pixel data is read-only (this is not
369 supported by afw). If `None`, only copy if the pixel data is
370 read-only.
371 """
372 import lsst.afw.image
373 import lsst.geom
375 array = self._array
376 if copy: 376 ↛ 377line 376 didn't jump to line 377 because the condition on line 376 was never true
377 array = array.copy()
378 elif not self._array.flags.writeable: 378 ↛ 379line 378 didn't jump to line 379 because the condition on line 378 was never true
379 if copy is None:
380 array = array.copy()
381 else:
382 raise TypeError("Cannot create a legacy lsst.afw.image.Image view into a read-only array.")
384 return lsst.afw.image.Image(
385 array,
386 deep=False,
387 dtype=array.dtype.type,
388 xy0=lsst.geom.Point2I(self._bbox.x.min, self._bbox.y.min),
389 )
391 @classmethod
392 def from_hdu_list(
393 cls,
394 hdu_list: astropy.io.fits.HDUList,
395 *,
396 fits_wcs_frame: Frame | None = DEFAULT_PIXEL_FRAME,
397 ) -> Image:
398 """Reconstruct an `~lsst.images.Image` from a cut-down ``lsst.images``
399 HDU list.
401 This reads only the first two HDUs (the primary HDU and the image
402 HDU), as written for the image-only cut-outs produced by
403 ``dax_images_cutout``: a real ``lsst.images`` file with its JSON-tree,
404 index, and any nested-archive HDUs dropped.
406 Parameters
407 ----------
408 hdu_list
409 HDU list whose first HDU is the primary header and whose second
410 HDU holds the image pixels.
411 fits_wcs_frame
412 Pixel-grid `~lsst.images.Frame` for the
413 `~lsst.images.SkyProjection` reconstructed from the image HDU's
414 FITS WCS. Defaults to a plain pixel frame; pass `None` to skip
415 attaching a projection.
417 Returns
418 -------
419 `~lsst.images.Image`
420 The reconstructed image, ready to be re-serialized as a normal
421 ``lsst.images`` file.
423 Notes
424 -----
425 The headers of the consumed HDUs are modified in place (WCS and other
426 interpreted cards are stripped), as in `read_legacy`.
427 """
428 opaque_metadata = fits.FitsOpaqueMetadata()
429 opaque_metadata.add_cutdown_primary_header(hdu_list[0].header)
430 result = cls._read_legacy_hdu(
431 hdu_list[1], opaque_metadata, preserve_bintable=None, fits_wcs_frame=fits_wcs_frame
432 )
433 result._opaque_metadata = opaque_metadata
434 return result
436 @staticmethod
437 def read_legacy(
438 uri: ResourcePathExpression,
439 *,
440 preserve_quantization: bool = False,
441 ext: str | int = 1,
442 fits_wcs_frame: Frame | None = None,
443 ) -> Image:
444 """Read a FITS file written by `lsst.afw.image.Image.writeFits`.
446 Parameters
447 ----------
448 uri
449 URI or file name.
450 preserve_quantization
451 If `True`, ensure that writing the image back out again will
452 exactly preserve quantization-compressed pixel values. This causes
453 the arrays to be marked as read-only and stores the original binary
454 table data for those planes in memory. If the `Image` is copied,
455 the precompressed pixel values are not transferred to the copy.
456 ext
457 Name or index of the FITS HDU to read.
458 fits_wcs_frame
459 If not `None` and the HDU containing the image has a FITS WCS,
460 attach a `SkyProjection` to the returned image by converting that
461 WCS.
462 """
463 opaque_metadata = fits.FitsOpaqueMetadata()
464 with ExitStack() as exit_stack:
465 fs, fspath = ResourcePath(uri).to_fsspec()
466 stream = exit_stack.enter_context(fs.open(fspath))
467 hdu_list = exit_stack.enter_context(astropy.io.fits.open(stream))
468 opaque_metadata.extract_legacy_primary_header(hdu_list[0].header)
469 bintable_hdu: astropy.io.fits.BinTableHDU | None = None
470 if preserve_quantization:
471 bintable_stream = exit_stack.enter_context(fs.open(fspath))
472 bintable_hdu_list = exit_stack.enter_context(
473 astropy.io.fits.open(bintable_stream, disable_image_compression=True)
474 )
475 bintable_hdu = bintable_hdu_list[ext]
476 result = Image._read_legacy_hdu(
477 hdu_list[ext], opaque_metadata, preserve_bintable=bintable_hdu, fits_wcs_frame=fits_wcs_frame
478 )
479 result._opaque_metadata = opaque_metadata
480 return result
482 @staticmethod
483 def _read_legacy_hdu(
484 hdu: astropy.io.fits.ImageHDU | astropy.io.fits.CompImageHDU,
485 opaque_metadata: fits.FitsOpaqueMetadata,
486 *,
487 preserve_bintable: astropy.io.fits.BinTableHDU | None,
488 fits_wcs_frame: Frame | None = None,
489 ) -> Image:
490 unit: astropy.units.UnitBase | None = None
491 if (fits_unit := hdu.header.pop("BUNIT", None)) is not None:
492 try:
493 unit = astropy.units.Unit(fits_unit, format="fits")
494 except ValueError:
495 # Accept non-FITS units by assuming Astropy can still figure
496 # them out if we don't specify the format.
497 unit = astropy.units.Unit(fits_unit)
498 if opaque_metadata.get_instrumental_unit() == astropy.units.electron: 498 ↛ 500line 498 didn't jump to line 500 because the condition on line 498 was never true
499 # Fix incorrect BUNIT='adu' in LSST preliminary_visit_image.
500 if unit == astropy.units.adu:
501 unit = astropy.units.electron
502 if unit == astropy.units.adu**2:
503 unit = astropy.units.electron**2
504 yx0 = fits.read_yx0(hdu.header)
505 hdu.header.remove("LTV1", ignore_missing=True)
506 hdu.header.remove("LTV2", ignore_missing=True)
507 read_only: bool = False
508 if preserve_bintable is not None: 508 ↛ 509line 508 didn't jump to line 509 because the condition on line 508 was never true
509 opaque_metadata.precompressed[hdu.name] = fits.PrecompressedImage.from_bintable(preserve_bintable)
510 read_only = True
511 sky_projection: SkyProjection | None = None
512 if fits_wcs_frame is not None:
513 try:
514 fits_wcs = astropy.wcs.WCS(hdu.header)
515 except KeyError:
516 pass
517 else:
518 sky_projection = SkyProjection.from_fits_wcs(
519 fits_wcs, pixel_frame=fits_wcs_frame, x0=yx0.x, y0=yx0.y
520 )
521 image = Image(hdu.data, yx0=yx0, unit=unit, sky_projection=sky_projection)
522 if read_only: 522 ↛ 523line 522 didn't jump to line 523 because the condition on line 522 was never true
523 image._array.flags["WRITEABLE"] = False
524 fits.strip_wcs_cards(hdu.header)
525 hdu.header.strip()
526 hdu.header.remove("EXTTYPE", ignore_missing=True)
527 hdu.header.remove("INHERIT", ignore_missing=True)
528 hdu.header.remove("UZSCALE", ignore_missing=True)
529 opaque_metadata.add_header(hdu.header)
530 return image
533class ImageSerializationModel[P: pydantic.BaseModel](ArchiveTree):
534 """Pydantic model used to represent the serialized form of an `.Image`."""
536 SCHEMA_NAME: ClassVar[str] = "image"
537 SCHEMA_VERSION: ClassVar[str] = "1.0.0"
538 MIN_READ_VERSION: ClassVar[int] = 1
539 PUBLIC_TYPE: ClassVar[type] = Image
541 data: ArrayReferenceQuantityModel | ArrayReferenceModel | InlineArrayModel | InlineArrayQuantityModel = (
542 pydantic.Field(description="Reference to pixel data.")
543 )
544 yx0: list[int] = pydantic.Field(
545 description="Coordinate of the first pixels in the array, ordered (y, x)."
546 )
547 sky_projection: SkyProjectionSerializationModel[P] | None = pydantic.Field(
548 default=None,
549 exclude_if=is_none,
550 description="Projection that maps the logical pixel grid onto the sky.",
551 )
553 @property
554 def bbox(self) -> Box:
555 """The bounding box of the image."""
556 match self.data:
557 case ArrayReferenceQuantityModel() | InlineArrayQuantityModel():
558 shape = self.data.value.shape
559 case ArrayReferenceModel() | InlineArrayModel(): 559 ↛ 561line 559 didn't jump to line 561 because the pattern on line 559 always matched
560 shape = self.data.shape
561 return Box.from_shape(shape, self.yx0)
563 def deserialize(
564 self,
565 archive: InputArchive[Any],
566 *,
567 bbox: Box | None = None,
568 strip_header: Callable[[astropy.io.fits.Header], None] = no_header_updates,
569 **kwargs: Any,
570 ) -> Image:
571 """Deserialize an image from an input archive.
573 Parameters
574 ----------
575 archive
576 Archive to read from.
577 bbox
578 Bounding box of a subimage to read instead.
579 strip_header
580 A callable that strips out any FITS header cards added by the
581 ``update_header`` argument in the corresponding call to
582 `Image.serialize`.
583 **kwargs
584 Unsupported keyword arguments are accepted only to provide better
585 error messages (raising `serialization.InvalidParameterError`).
586 """
587 if kwargs: 587 ↛ 588line 587 didn't jump to line 588 because the condition on line 587 was never true
588 raise InvalidParameterError(f"Unrecognized parameters for Image: {set(kwargs.keys())}.")
589 array_model: ArrayReferenceModel | InlineArrayModel
590 unit: astropy.units.UnitBase | None = None
591 if isinstance(self.data, ArrayReferenceQuantityModel | InlineArrayQuantityModel):
592 array_model = self.data.value
593 unit = self.data.unit
594 else:
595 array_model = self.data
597 def _strip_header(header: astropy.io.fits.Header) -> None:
598 if unit is not None:
599 header.pop("BUNIT", None)
600 fits.strip_wcs_cards(header)
601 strip_header(header)
603 slices = bbox.slice_within(self.bbox) if bbox is not None else ...
604 array = archive.get_array(array_model, strip_header=_strip_header, slices=slices)
605 sky_projection = self.sky_projection.deserialize(archive) if self.sky_projection is not None else None
606 return Image(
607 array,
608 yx0=self.yx0 if bbox is None else bbox.start,
609 unit=unit,
610 sky_projection=sky_projection,
611 )._finish_deserialize(self)
613 def deserialize_component(self, component: str, archive: InputArchive[Any], **kwargs: Any) -> Any:
614 if kwargs:
615 raise InvalidParameterError(f"Unsupported parameters for Image components: {set(kwargs.keys())}.")
616 return super().deserialize_component(component, archive)