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