Coverage for python/lsst/images/_image.py: 57%
225 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-23 09:23 +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 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(
349 legacy: LegacyImage,
350 unit: astropy.units.UnitBase | None = None,
351 *,
352 sky_projection: SkyProjection[Any] | None = None,
353 ) -> Image:
354 """Convert from an `lsst.afw.image.Image` instance.
356 Parameters
357 ----------
358 legacy
359 An `lsst.afw.image.Image` instance that will share pixel data with
360 the returned object.
361 unit
362 Units of the image.
363 sky_projection
364 Projection from pixels to xky.
365 """
366 return Image(
367 legacy.array, yx0=YX(y=legacy.getY0(), x=legacy.getX0()), unit=unit, sky_projection=sky_projection
368 )
370 def to_legacy(self, *, copy: bool | None = None) -> LegacyImage:
371 """Convert to an `lsst.afw.image.Image` instance.
373 Parameters
374 ----------
375 copy
376 If `True`, always copy the pixel data. If `False`, return a view,
377 and raise `TypeError` if the pixel data is read-only (this is not
378 supported by afw). If `None`, only copy if the pixel data is
379 read-only.
380 """
381 import lsst.afw.image
382 import lsst.geom
384 array = self._array
385 if copy: 385 ↛ 386line 385 didn't jump to line 386 because the condition on line 385 was never true
386 array = array.copy()
387 elif not self._array.flags.writeable: 387 ↛ 388line 387 didn't jump to line 388 because the condition on line 387 was never true
388 if copy is None:
389 array = array.copy()
390 else:
391 raise TypeError("Cannot create a legacy lsst.afw.image.Image view into a read-only array.")
393 return lsst.afw.image.Image(
394 array,
395 deep=False,
396 dtype=array.dtype.type,
397 xy0=lsst.geom.Point2I(self._bbox.x.min, self._bbox.y.min),
398 )
400 @classmethod
401 def from_hdu_list(
402 cls,
403 hdu_list: astropy.io.fits.HDUList,
404 *,
405 fits_wcs_frame: Frame | None = DEFAULT_PIXEL_FRAME,
406 ) -> Image:
407 """Reconstruct an `~lsst.images.Image` from a cut-down ``lsst.images``
408 HDU list.
410 This reads only the first two HDUs (the primary HDU and the image
411 HDU), as written for the image-only cut-outs produced by
412 ``dax_images_cutout``: a real ``lsst.images`` file with its JSON-tree,
413 index, and any nested-archive HDUs dropped.
415 Parameters
416 ----------
417 hdu_list
418 HDU list whose first HDU is the primary header and whose second
419 HDU holds the image pixels.
420 fits_wcs_frame
421 Pixel-grid `~lsst.images.Frame` for the
422 `~lsst.images.SkyProjection` reconstructed from the image HDU's
423 FITS WCS. Defaults to a plain pixel frame; pass `None` to skip
424 attaching a projection.
426 Returns
427 -------
428 `~lsst.images.Image`
429 The reconstructed image, ready to be re-serialized as a normal
430 ``lsst.images`` file.
432 Notes
433 -----
434 The headers of the consumed HDUs are modified in place (WCS and other
435 interpreted cards are stripped), as in `read_legacy`.
436 """
437 opaque_metadata = fits.FitsOpaqueMetadata()
438 opaque_metadata.add_cutdown_primary_header(hdu_list[0].header)
439 result = cls._read_legacy_hdu(
440 hdu_list[1], opaque_metadata, preserve_bintable=None, fits_wcs_frame=fits_wcs_frame
441 )
442 result._opaque_metadata = opaque_metadata
443 return result
445 @staticmethod
446 def read_legacy(
447 uri: ResourcePathExpression,
448 *,
449 preserve_quantization: bool = False,
450 ext: str | int = 1,
451 fits_wcs_frame: Frame | None = None,
452 ) -> Image:
453 """Read a FITS file written by `lsst.afw.image.Image.writeFits`.
455 Parameters
456 ----------
457 uri
458 URI or file name.
459 preserve_quantization
460 If `True`, ensure that writing the image back out again will
461 exactly preserve quantization-compressed pixel values. This causes
462 the arrays to be marked as read-only and stores the original binary
463 table data for those planes in memory. If the `Image` is copied,
464 the precompressed pixel values are not transferred to the copy.
465 ext
466 Name or index of the FITS HDU to read.
467 fits_wcs_frame
468 If not `None` and the HDU containing the image has a FITS WCS,
469 attach a `SkyProjection` to the returned image by converting that
470 WCS.
471 """
472 opaque_metadata = fits.FitsOpaqueMetadata()
473 with ExitStack() as exit_stack:
474 fs, fspath = ResourcePath(uri).to_fsspec()
475 stream = exit_stack.enter_context(fs.open(fspath))
476 hdu_list = exit_stack.enter_context(astropy.io.fits.open(stream))
477 opaque_metadata.extract_legacy_primary_header(hdu_list[0].header)
478 bintable_hdu: astropy.io.fits.BinTableHDU | None = None
479 if preserve_quantization:
480 bintable_stream = exit_stack.enter_context(fs.open(fspath))
481 bintable_hdu_list = exit_stack.enter_context(
482 astropy.io.fits.open(bintable_stream, disable_image_compression=True)
483 )
484 bintable_hdu = bintable_hdu_list[ext]
485 result = Image._read_legacy_hdu(
486 hdu_list[ext], opaque_metadata, preserve_bintable=bintable_hdu, fits_wcs_frame=fits_wcs_frame
487 )
488 result._opaque_metadata = opaque_metadata
489 return result
491 @staticmethod
492 def _read_legacy_hdu(
493 hdu: astropy.io.fits.ImageHDU | astropy.io.fits.CompImageHDU,
494 opaque_metadata: fits.FitsOpaqueMetadata,
495 *,
496 preserve_bintable: astropy.io.fits.BinTableHDU | None,
497 fits_wcs_frame: Frame | None = None,
498 ) -> Image:
499 unit: astropy.units.UnitBase | None = None
500 if (fits_unit := hdu.header.pop("BUNIT", None)) is not None:
501 try:
502 unit = astropy.units.Unit(fits_unit, format="fits")
503 except ValueError:
504 # Accept non-FITS units by assuming Astropy can still figure
505 # them out if we don't specify the format.
506 unit = astropy.units.Unit(fits_unit)
507 if opaque_metadata.get_instrumental_unit() == astropy.units.electron: 507 ↛ 509line 507 didn't jump to line 509 because the condition on line 507 was never true
508 # Fix incorrect BUNIT='adu' in LSST preliminary_visit_image.
509 if unit == astropy.units.adu:
510 unit = astropy.units.electron
511 if unit == astropy.units.adu**2:
512 unit = astropy.units.electron**2
513 yx0 = fits.read_yx0(hdu.header)
514 hdu.header.remove("LTV1", ignore_missing=True)
515 hdu.header.remove("LTV2", ignore_missing=True)
516 read_only: bool = False
517 if preserve_bintable is not None: 517 ↛ 518line 517 didn't jump to line 518 because the condition on line 517 was never true
518 opaque_metadata.precompressed[hdu.name] = fits.PrecompressedImage.from_bintable(preserve_bintable)
519 read_only = True
520 sky_projection: SkyProjection | None = None
521 if fits_wcs_frame is not None:
522 try:
523 fits_wcs = astropy.wcs.WCS(hdu.header)
524 except KeyError:
525 pass
526 else:
527 sky_projection = SkyProjection.from_fits_wcs(
528 fits_wcs, pixel_frame=fits_wcs_frame, x0=yx0.x, y0=yx0.y
529 )
530 image = Image(hdu.data, yx0=yx0, unit=unit, sky_projection=sky_projection)
531 if read_only: 531 ↛ 532line 531 didn't jump to line 532 because the condition on line 531 was never true
532 image._array.flags["WRITEABLE"] = False
533 fits.strip_wcs_cards(hdu.header)
534 hdu.header.strip()
535 hdu.header.remove("EXTTYPE", ignore_missing=True)
536 hdu.header.remove("INHERIT", ignore_missing=True)
537 hdu.header.remove("UZSCALE", ignore_missing=True)
538 opaque_metadata.add_header(hdu.header)
539 return image
542class ImageSerializationModel[P: pydantic.BaseModel](ArchiveTree):
543 """Pydantic model used to represent the serialized form of an `.Image`."""
545 SCHEMA_NAME: ClassVar[str] = "image"
546 SCHEMA_VERSION: ClassVar[str] = "1.0.0"
547 MIN_READ_VERSION: ClassVar[int] = 1
548 PUBLIC_TYPE: ClassVar[type] = Image
550 data: ArrayReferenceQuantityModel | ArrayReferenceModel | InlineArrayModel | InlineArrayQuantityModel = (
551 pydantic.Field(description="Reference to pixel data.")
552 )
553 yx0: list[int] = pydantic.Field(
554 description="Coordinate of the first pixels in the array, ordered (y, x)."
555 )
556 sky_projection: SkyProjectionSerializationModel[P] | None = pydantic.Field(
557 default=None,
558 exclude_if=is_none,
559 description="Projection that maps the logical pixel grid onto the sky.",
560 )
562 @property
563 def bbox(self) -> Box:
564 """The bounding box of the image."""
565 match self.data:
566 case ArrayReferenceQuantityModel() | InlineArrayQuantityModel():
567 shape = self.data.value.shape
568 case ArrayReferenceModel() | InlineArrayModel(): 568 ↛ 570line 568 didn't jump to line 570 because the pattern on line 568 always matched
569 shape = self.data.shape
570 return Box.from_shape(shape, self.yx0)
572 def deserialize(
573 self,
574 archive: InputArchive[Any],
575 *,
576 bbox: Box | None = None,
577 strip_header: Callable[[astropy.io.fits.Header], None] = no_header_updates,
578 **kwargs: Any,
579 ) -> Image:
580 """Deserialize an image from an input archive.
582 Parameters
583 ----------
584 archive
585 Archive to read from.
586 bbox
587 Bounding box of a subimage to read instead.
588 strip_header
589 A callable that strips out any FITS header cards added by the
590 ``update_header`` argument in the corresponding call to
591 `Image.serialize`.
592 **kwargs
593 Unsupported keyword arguments are accepted only to provide better
594 error messages (raising `serialization.InvalidParameterError`).
595 """
596 if kwargs: 596 ↛ 597line 596 didn't jump to line 597 because the condition on line 596 was never true
597 raise InvalidParameterError(f"Unrecognized parameters for Image: {set(kwargs.keys())}.")
598 array_model: ArrayReferenceModel | InlineArrayModel
599 unit: astropy.units.UnitBase | None = None
600 if isinstance(self.data, ArrayReferenceQuantityModel | InlineArrayQuantityModel):
601 array_model = self.data.value
602 unit = self.data.unit
603 else:
604 array_model = self.data
606 def _strip_header(header: astropy.io.fits.Header) -> None:
607 if unit is not None:
608 header.pop("BUNIT", None)
609 fits.strip_wcs_cards(header)
610 strip_header(header)
612 slices = bbox.slice_within(self.bbox) if bbox is not None else ...
613 array = archive.get_array(array_model, strip_header=_strip_header, slices=slices)
614 sky_projection = self.sky_projection.deserialize(archive) if self.sky_projection is not None else None
615 return Image(
616 array,
617 yx0=self.yx0 if bbox is None else bbox.start,
618 unit=unit,
619 sky_projection=sky_projection,
620 )._finish_deserialize(self)
622 def deserialize_component(self, component: str, archive: InputArchive[Any], **kwargs: Any) -> Any:
623 if kwargs:
624 raise InvalidParameterError(f"Unsupported parameters for Image components: {set(kwargs.keys())}.")
625 return super().deserialize_component(component, archive)