Coverage for python / lsst / images / _image.py: 25%
217 statements
« prev ^ index » next coverage.py v7.14.0, created at 2026-05-16 07:54 +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 Any, ClassVar, final
21import astropy.io.fits
22import astropy.units
23import astropy.wcs
24import numpy as np
25import numpy.typing as npt
26import pydantic
27from astro_metadata_translator import ObservationInfo
29from lsst.resources import ResourcePath, ResourcePathExpression
31from . import fits
32from ._generalized_image import GeneralizedImage
33from ._geom import YX, Box
34from ._transforms import Frame, Projection, ProjectionSerializationModel
35from .serialization import (
36 ArchiveTree,
37 ArrayReferenceModel,
38 ArrayReferenceQuantityModel,
39 InlineArrayModel,
40 InlineArrayQuantityModel,
41 InputArchive,
42 MetadataValue,
43 OutputArchive,
44 no_header_updates,
45)
46from .utils import is_none
49@final
50class Image(GeneralizedImage):
51 """A 2-d array that may be augmented with units and a nonzero origin.
53 Parameters
54 ----------
55 array_or_fill
56 Array or fill value for the image. If a fill value, ``bbox`` or
57 ``shape`` must be provided.
58 bbox
59 Bounding box for the image.
60 start
61 Logical coordinates of the first pixel in the array, ordered ``y``,
62 ``x`` (unless an `XY` instance is passed). Ignored if
63 ``bbox`` is provided. Defaults to zeros.
64 shape
65 Leading dimensions of the array, ordered ``y``, ``x`` (unless an `XY`
66 instance is passed). Only needed if ``array_or_fill`` is not an
67 array and ``bbox`` is not provided. Like the bbox, this does not
68 include the last dimension of the array.
69 dtype
70 Pixel data type override.
71 unit
72 Units for the image's pixel values.
73 projection
74 Projection that maps the pixel grid to the sky.
75 obs_info
76 General information about the associated observation in standardized
77 form.
78 metadata
79 Arbitrary flexible metadata to associate with the image.
81 Notes
82 -----
83 Indexing the `array` attribute of an `Image` does not take into account its
84 ``start`` offset, but accessing a subimage by indexing an `Image` with a
85 `Box` does, and the `bbox` of the subimage is set to match its location
86 within the original image.
88 Indexed assignment to a subimage requires consistency between the
89 coordinate systems and units of both operands, but it will automatically
90 select a subimage of the right-hand side and convert compatible units when
91 possible. In other words::
93 a[box] = b
95 is a shortcut for
97 a[box].quantity = b[box].quantity
99 An ellipsis (``...``) can be used instead of a `Box` to assign to the full
100 image.
101 """
103 def __init__(
104 self,
105 array_or_fill: np.ndarray | int | float = 0,
106 /,
107 *,
108 bbox: Box | None = None,
109 start: Sequence[int] | None = None,
110 shape: Sequence[int] | None = None,
111 dtype: npt.DTypeLike | None = None,
112 unit: astropy.units.UnitBase | None = None,
113 projection: Projection[Any] | None = None,
114 obs_info: ObservationInfo | None = None,
115 metadata: dict[str, MetadataValue] | None = None,
116 ):
117 super().__init__(metadata)
118 if isinstance(array_or_fill, np.ndarray):
119 if dtype is not None:
120 array = np.array(array_or_fill, dtype=dtype)
121 else:
122 array = array_or_fill
123 if bbox is None:
124 bbox = Box.from_shape(array.shape, start=start)
125 elif bbox.shape != array.shape:
126 raise ValueError(
127 f"Explicit bbox shape {bbox.shape} does not match array with shape {array.shape}."
128 )
129 if shape is not None and shape != array.shape:
130 raise ValueError(f"Explicit shape {shape} does not match array with shape {array.shape}.")
131 else:
132 if bbox is None:
133 if shape is None:
134 raise TypeError("No bbox, shape, or array provided.")
135 bbox = Box.from_shape(shape, start=start)
136 elif shape is not None and shape != bbox.shape:
137 raise ValueError(f"Explicit shape {shape} does not match bbox shape {bbox.shape}.")
138 array = np.full(bbox.shape, array_or_fill, dtype=dtype)
139 self._array: np.ndarray = array
140 self._bbox: Box = bbox
141 self._unit = unit
142 self._projection = projection
143 self._obs_info = obs_info
145 @property
146 def array(self) -> np.ndarray:
147 """The low-level array (`numpy.ndarray`).
149 Assigning to this attribute modifies the existing array in place; the
150 bounding box and underlying data pointer are never changed.
151 """
152 return self._array
154 @array.setter
155 def array(self, value: np.ndarray | int | float) -> None:
156 self._array[...] = value
158 @property
159 def quantity(self) -> astropy.units.Quantity:
160 """The low-level array with units (`astropy.units.Quantity`).
162 Assigning to this attribute modifies the existing array in place; the
163 bounding box and underlying data pointer are never changed.
164 """
165 return astropy.units.Quantity(self._array, self._unit, copy=False)
167 @quantity.setter
168 def quantity(self, value: astropy.units.Quantity) -> None:
169 self.quantity[...] = value
171 @property
172 def bbox(self) -> Box:
173 """Bounding box for the image (`Box`)."""
174 return self._bbox
176 @property
177 def unit(self) -> astropy.units.UnitBase | None:
178 """Units for the image's pixel values (`astropy.units.Unit` or
179 `None`).
180 """
181 return self._unit
183 @property
184 def projection(self) -> Projection[Any] | None:
185 """The projection that maps this image's pixel grid to the sky
186 (`Projection` | `None`).
188 Notes
189 -----
190 The pixel coordinates used by this projection account for the bounding
191 box ``start``; they are not just array indices.
192 """
193 return self._projection
195 @property
196 def obs_info(self) -> ObservationInfo | None:
197 """General information about the associated observation in standard
198 form. (`~astro_metadata_translator.ObservationInfo` | `None`).
199 """
200 return self._obs_info
202 def __getitem__(self, bbox: Box | EllipsisType) -> Image:
203 if bbox is ...:
204 return self
205 super().__getitem__(bbox)
206 indices = bbox.slice_within(self._bbox)
207 return self._transfer_metadata(
208 Image(
209 self._array[indices],
210 bbox=bbox,
211 unit=self._unit,
212 projection=self._projection,
213 obs_info=self._obs_info,
214 ),
215 bbox=bbox,
216 )
218 def __setitem__(self, bbox: Box | EllipsisType, value: Image) -> None:
219 self[bbox].quantity[...] = value.quantity
221 def __str__(self) -> str:
222 return f"Image({self.bbox!s}, {self.array.dtype.type.__name__})"
224 def __repr__(self) -> str:
225 return f"Image(..., bbox={self.bbox!r}, dtype={self.array.dtype!r})"
227 def __eq__(self, other: object) -> bool:
228 if not isinstance(other, Image):
229 return NotImplemented
230 return (
231 self._bbox == other._bbox
232 and self._unit == other._unit
233 and np.array_equal(self._array, other._array, equal_nan=True)
234 )
236 def copy(self) -> Image:
237 return self._transfer_metadata(
238 Image(
239 self._array.copy(),
240 bbox=self._bbox,
241 unit=self._unit,
242 projection=self._projection,
243 obs_info=self._obs_info,
244 ),
245 copy=True,
246 )
248 def view(
249 self,
250 *,
251 unit: astropy.units.UnitBase | None | EllipsisType = ...,
252 projection: Projection | None | EllipsisType = ...,
253 start: Sequence[int] | EllipsisType = ...,
254 obs_info: ObservationInfo | None | EllipsisType = ...,
255 ) -> Image:
256 """Make a view of the image, with optional updates."""
257 if unit is ...:
258 unit = self._unit
259 if projection is ...:
260 projection = self._projection
261 if start is ...:
262 start = self._bbox.start
263 if obs_info is ...:
264 obs_info = self._obs_info
265 return self._transfer_metadata(
266 Image(self._array, start=start, unit=unit, projection=projection, obs_info=obs_info)
267 )
269 def serialize[P: pydantic.BaseModel](
270 self,
271 archive: OutputArchive[P],
272 *,
273 update_header: Callable[[astropy.io.fits.Header], None] = no_header_updates,
274 save_projection: bool = True,
275 save_obs_info: bool = True,
276 add_offset_wcs: str | None = "A",
277 ) -> ImageSerializationModel[P]:
278 """Serialize the image to an output archive.
280 Parameters
281 ----------
282 archive
283 Archive to write to.
284 update_header
285 A callback that will be given the FITS header for the HDU
286 containing this image in order to add keys to it. This callback
287 may be provided but will not be called if the output format is not
288 FITS.
289 save_projection
290 If `True`, save the `Projection` attached to the image, if there
291 is one. This does not affect whether a FITS WCS corresponding to
292 the projection is written (it always is, if available, and if
293 ``add_offset_wcs`` is not ``" "``).
294 save_obs_info
295 If `True`, save the
296 `~astro_metadata_translator.ObservationInfo` attached to the
297 image, if there is one.
298 add_offset_wcs
299 A FITS WCS single-character suffix to use when adding a linear
300 WCS that maps the FITS array to the logical pixel coordinates
301 defined by ``bbox.start``. Set to `None` to not write this WCS.
302 If this is set to ``" "``, it will prevent the `Projection` from
303 being saved as a FITS WCS.
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.
313 pass
314 if self.projection is not None and add_offset_wcs != " ":
315 if self.fits_wcs:
316 header.update(self.fits_wcs.to_header(relax=True))
317 if add_offset_wcs is not None:
318 fits.add_offset_wcs(header, x=self.bbox.x.start, y=self.bbox.y.start, key=add_offset_wcs)
320 array_model = archive.add_array(self.array, update_header=_update_header)
321 serialized_projection: ProjectionSerializationModel[P] | None = None
322 if save_projection and self.projection is not None:
323 serialized_projection = archive.serialize_direct("projection", self.projection.serialize)
324 data = array_model if self.unit is None else array_model.with_units(self.unit)
325 return ImageSerializationModel.model_construct(
326 data=data,
327 start=list(self.bbox.start),
328 projection=serialized_projection,
329 obs_info=self._obs_info if save_obs_info else None,
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 def write_fits(
348 self,
349 filename: str,
350 *,
351 compression: fits.FitsCompressionOptions | None = fits.FitsCompressionOptions.DEFAULT,
352 compression_seed: int | None = None,
353 ) -> None:
354 """Write the image to a FITS file.
356 Parameters
357 ----------
358 filename
359 Name of the file to write to. Must be a local file.
360 compression
361 Compression options.
362 compression_seed
363 A FITS tile compression seed to use whenever the configured
364 compression seed is `None` or (for backwards compatibility) ``0``.
365 """
366 compression_options = {}
367 if compression is not fits.FitsCompressionOptions.DEFAULT:
368 compression_options[self._archive_default_name] = compression
369 fits.write(self, filename, compression_options, compression_seed=compression_seed)
371 @staticmethod
372 def read_fits(url: ResourcePathExpression, *, bbox: Box | None = None) -> Image:
373 """Read an image from a FITS file.
375 Parameters
376 ----------
377 url
378 URL of the file to read; may be any type supported by
379 `lsst.resources.ResourcePath`.
380 bbox
381 Bounding box of a subimage to read instead.
382 """
383 return fits.read(Image, url, bbox=bbox).deserialized
385 @staticmethod
386 def from_legacy(legacy: Any, unit: astropy.units.UnitBase | None = None) -> Image:
387 """Convert from an `lsst.afw.image.Image` instance.
389 Parameters
390 ----------
391 legacy
392 An `lsst.afw.image.Image` instance that will share pixel data with
393 the returned object.
394 unit
395 Units of the image.
396 """
397 return Image(legacy.array, start=(legacy.getY0(), legacy.getX0()), unit=unit)
399 def to_legacy(self, *, copy: bool | None = None) -> Any:
400 """Convert to an `lsst.afw.image.Image` instance.
402 Parameters
403 ----------
404 copy
405 If `True`, always copy the pixel data. If `False`, return a view,
406 and raise `TypeError` if the pixel data is read-only (this is not
407 supported by afw). If `None`, onyl if the pixel data is
408 read-only.
409 """
410 import lsst.afw.image
411 import lsst.geom
413 array = self._array
414 if copy:
415 array = array.copy()
416 elif not self._array.flags.writeable:
417 if copy is None:
418 array = array.copy()
419 else:
420 raise TypeError("Cannot create a legacy lsst.afw.image.Image view into a read-only array.")
422 return lsst.afw.image.Image(
423 array,
424 deep=False,
425 dtype=array.dtype.type,
426 xy0=lsst.geom.Point2I(self._bbox.x.min, self._bbox.y.min),
427 )
429 @staticmethod
430 def read_legacy(
431 uri: ResourcePathExpression,
432 *,
433 preserve_quantization: bool = False,
434 ext: str | int = 1,
435 fits_wcs_frame: Frame | None = None,
436 ) -> Image:
437 """Read a FITS file written by `lsst.afw.image.Image.writeFits`.
439 Parameters
440 ----------
441 uri
442 URI or file name.
443 preserve_quantization
444 If `True`, ensure that writing the image back out again will
445 exactly preserve quantization-compressed pixel values. This causes
446 the arrays to be marked as read-only and stores the original binary
447 table data for those planes in memory. If the `Image` is copied,
448 the precompressed pixel values are not transferred to the copy.
449 ext
450 Name or index of the FITS HDU to read.
451 fits_wcs_frame
452 If not `None` and the HDU containing the image has a FITS WCS,
453 attach a `Projection` to the returned image by converting that
454 WCS.
455 """
456 opaque_metadata = fits.FitsOpaqueMetadata()
457 with ExitStack() as exit_stack:
458 fs, fspath = ResourcePath(uri).to_fsspec()
459 stream = exit_stack.enter_context(fs.open(fspath))
460 hdu_list = exit_stack.enter_context(astropy.io.fits.open(stream))
461 opaque_metadata.extract_legacy_primary_header(hdu_list[0].header)
462 bintable_hdu: astropy.io.fits.BinTableHDU | None = None
463 if preserve_quantization:
464 bintable_stream = exit_stack.enter_context(fs.open(fspath))
465 bintable_hdu_list = exit_stack.enter_context(
466 astropy.io.fits.open(bintable_stream, disable_image_compression=True)
467 )
468 bintable_hdu = bintable_hdu_list[ext]
469 result = Image._read_legacy_hdu(
470 hdu_list[ext], opaque_metadata, preserve_bintable=bintable_hdu, fits_wcs_frame=fits_wcs_frame
471 )
472 result._opaque_metadata = opaque_metadata
473 return result
475 @staticmethod
476 def _read_legacy_hdu(
477 hdu: astropy.io.fits.ImageHDU | astropy.io.fits.CompImageHDU,
478 opaque_metadata: fits.FitsOpaqueMetadata,
479 *,
480 preserve_bintable: astropy.io.fits.BinTableHDU | None,
481 fits_wcs_frame: Frame | None = None,
482 ) -> Image:
483 unit: astropy.units.UnitBase | None = None
484 if (fits_unit := hdu.header.pop("BUNIT", None)) is not None:
485 unit = astropy.units.Unit(fits_unit, format="fits")
486 dx: int = hdu.header.pop("LTV1")
487 dy: int = hdu.header.pop("LTV2")
488 start = YX(y=-dy, x=-dx)
489 read_only: bool = False
490 if preserve_bintable is not None:
491 opaque_metadata.precompressed[hdu.name] = fits.PrecompressedImage.from_bintable(preserve_bintable)
492 read_only = True
493 projection: Projection | None = None
494 if fits_wcs_frame is not None:
495 try:
496 fits_wcs = astropy.wcs.WCS(hdu.header)
497 except KeyError:
498 pass
499 else:
500 projection = Projection.from_fits_wcs(
501 fits_wcs, pixel_frame=fits_wcs_frame, x0=start.x, y0=start.y
502 )
503 image = Image(hdu.data, start=start, unit=unit, projection=projection)
504 if read_only:
505 image._array.flags["WRITEABLE"] = False
506 fits.strip_wcs_cards(hdu.header)
507 hdu.header.strip()
508 hdu.header.remove("EXTTYPE", ignore_missing=True)
509 hdu.header.remove("INHERIT", ignore_missing=True)
510 hdu.header.remove("UZSCALE", ignore_missing=True)
511 opaque_metadata.add_header(hdu.header)
512 return image
515class ImageSerializationModel[P: pydantic.BaseModel](ArchiveTree):
516 """Pydantic model used to represent the serialized form of an `.Image`."""
518 data: ArrayReferenceQuantityModel | ArrayReferenceModel | InlineArrayModel | InlineArrayQuantityModel = (
519 pydantic.Field(description="Reference to pixel data.")
520 )
521 start: list[int] = pydantic.Field(
522 description="Coordinate of the first pixels in the array, ordered (y, x)."
523 )
524 projection: ProjectionSerializationModel[P] | None = pydantic.Field(
525 default=None,
526 exclude_if=is_none,
527 description="Projection that maps the logical pixel grid onto the sky.",
528 )
529 obs_info: ObservationInfo | None = pydantic.Field(
530 default=None,
531 exclude_if=is_none,
532 description="Standardized description of image metadata",
533 )
535 @property
536 def bbox(self) -> Box:
537 """The bounding box of the image."""
538 match self.data:
539 case ArrayReferenceQuantityModel() | InlineArrayQuantityModel():
540 shape = self.data.value.shape
541 case ArrayReferenceModel() | InlineArrayModel():
542 shape = self.data.shape
543 return Box.from_shape(shape, self.start)
545 def deserialize(
546 self,
547 archive: InputArchive[Any],
548 *,
549 bbox: Box | None = None,
550 strip_header: Callable[[astropy.io.fits.Header], None] = no_header_updates,
551 ) -> Image:
552 """Deserialize an image from an input archive.
554 Parameters
555 ----------
556 archive
557 Archive to read from.
558 bbox
559 Bounding box of a subimage to read instead.
560 strip_header
561 A callable that strips out any FITS header cards added by the
562 ``update_header`` argument in the corresponding call to
563 `Image.serialize`.
564 """
565 array_model: ArrayReferenceModel | InlineArrayModel
566 unit: astropy.units.UnitBase | None = None
567 if isinstance(self.data, ArrayReferenceQuantityModel | InlineArrayQuantityModel):
568 array_model = self.data.value
569 unit = self.data.unit
570 else:
571 array_model = self.data
573 def _strip_header(header: astropy.io.fits.Header) -> None:
574 if unit is not None:
575 header.pop("BUNIT", None)
576 fits.strip_wcs_cards(header)
577 strip_header(header)
579 slices = bbox.slice_within(self.bbox) if bbox is not None else ...
580 array = archive.get_array(array_model, strip_header=_strip_header, slices=slices)
581 projection = self.projection.deserialize(archive) if self.projection is not None else None
582 return Image(
583 array,
584 start=self.start if bbox is None else bbox.start,
585 unit=unit,
586 projection=projection,
587 obs_info=self.obs_info,
588 )._finish_deserialize(self)