Coverage for python / lsst / images / _visit_image.py: 27%
227 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-02 07:16 -0800
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
14import warnings
16__all__ = ("VisitImage", "VisitImageSerializationModel")
18import functools
19from collections.abc import Callable, Mapping, MutableMapping
20from typing import Any, Literal, cast, overload
22import astropy.io.fits
23import astropy.units
24import astropy.wcs
25import pydantic
26from astro_metadata_translator import ObservationInfo, VisitInfoTranslator
28from ._geom import Box
29from ._image import Image, ImageSerializationModel
30from ._mask import Mask, MaskPlane, MaskSchema, MaskSerializationModel, get_legacy_visit_image_mask_planes
31from ._masked_image import MaskedImage, MaskedImageSerializationModel
32from ._transforms import DetectorFrame, Projection, ProjectionAstropyView, ProjectionSerializationModel
33from .fits import FitsOpaqueMetadata
34from .psfs import (
35 PiffSerializationModel,
36 PiffWrapper,
37 PointSpreadFunction,
38 PSFExSerializationModel,
39 PSFExWrapper,
40)
41from .serialization import ArchiveReadError, InputArchive, OutputArchive
44def _obs_info_from_md(md: MutableMapping[str, Any], visit_info: Any = None) -> ObservationInfo:
45 # Try to get an ObservationInfo from the primary header as if
46 # it's a raw header. Else fallback.
47 try:
48 obs_info = ObservationInfo.from_header(md, quiet=True)
49 except ValueError:
50 # Not known translator. Must fall back to visit info. If we have
51 # an actual VisitInfo, serialize it since we know that it will be
52 # complete.
53 if visit_info is not None:
54 from lsst.afw.image import setVisitInfoMetadata
55 from lsst.daf.base import PropertyList
57 pl = PropertyList()
58 setVisitInfoMetadata(pl, visit_info)
59 # Merge so that we still have access to butler provenance.
60 md.update(pl)
62 # Try the given header looking for VisitInfo hints.
63 # We get lots of warnings if nothing can be found. Currently
64 # no way to disable those without capturing them.
65 obs_info = ObservationInfo.from_header(md, translator_class=VisitInfoTranslator, quiet=True)
66 return obs_info
69class VisitImage(MaskedImage):
70 """A calibrated single-visit image.
72 Parameters
73 ----------
74 image
75 The main image plane. If this has a `Projection`, it will be used
76 for all planes unless a ``projection`` is passed separately.
77 mask
78 A bitmask image that annotates the main image plane. Must have the
79 same bounding box as ``image`` if provided. Any attached projection
80 is replaced (possibly by `None`).
81 variance
82 The per-pixel uncertainty of the main image as an image of variance
83 values. Must have the same bounding box as ``image`` if provided, and
84 its units must be the square of ``image.unit`` or `None`.
85 Values default to ``1.0``. Any attached projection is replaced
86 (possibly by `None`).
87 mask_schema
88 Schema for the mask plane. Must be provided if and only if ``mask`` is
89 not provided.
90 opaque_metadata
91 Opaque metadata obtained from reading this object from storage. It may
92 be provided when writing to storage to propagate that metadata and/or
93 preserve file-format-specific options (e.g. compression parameters).
94 projection
95 Projection that maps the pixel grid to the sky. Can only be `None` if
96 a projection is already attached to ``image``.
97 psf
98 Point-spread function model for this image, or an exception explaining
99 why it could not be read (to be raised if the PSF is requested later).
100 obs_info
101 General information about this visit in standardized form.
102 """
104 def __init__(
105 self,
106 image: Image,
107 *,
108 mask: Mask | None = None,
109 variance: Image | None = None,
110 mask_schema: MaskSchema | None = None,
111 projection: Projection[DetectorFrame] | None = None,
112 obs_info: ObservationInfo | None = None,
113 psf: PointSpreadFunction | ArchiveReadError,
114 ):
115 super().__init__(
116 image,
117 mask=mask,
118 variance=variance,
119 mask_schema=mask_schema,
120 projection=projection,
121 obs_info=obs_info,
122 )
123 if self.image.unit is None:
124 raise TypeError("The image component of a VisitImage must have units.")
125 if self.image.projection is None:
126 raise TypeError("The projection component of a VisitImage cannot be None.")
127 if self.image.obs_info is None:
128 raise TypeError("The observation info component of a VisitImage cannot be None.")
129 if not isinstance(self.image.projection.pixel_frame, DetectorFrame):
130 raise TypeError("The projection's pixel frame must be a DetectorFrame for VisitImage.")
131 self._psf = psf
133 @property
134 def unit(self) -> astropy.units.UnitBase:
135 """The units of the image plane (`astropy.units.Unit`)."""
136 return cast(astropy.units.UnitBase, super().unit)
138 @property
139 def projection(self) -> Projection[DetectorFrame]:
140 """The projection that maps the pixel grid to the sky
141 (`Projection` [`DetectorFrame`]).
142 """
143 return cast(Projection[DetectorFrame], super().projection)
145 @property
146 def obs_info(self) -> ObservationInfo:
147 """General information about this observation in standard form.
148 (`~astro_metadata_translator.ObservationInfo`).
149 """
150 obs_info = self.image.obs_info
151 assert obs_info is not None
152 return obs_info
154 @property
155 def astropy_wcs(self) -> ProjectionAstropyView:
156 """An Astropy WCS for the pixel arrays (`ProjectionAstropyView`).
158 Notes
159 -----
160 As expected for Astropy WCS objects, this defines pixel coordinates
161 such that the first row and column in the arrays are ``(0, 0)``, not
162 ``bbox.start``, as is the case for `projection`.
164 This object satisfies the `astropy.wcs.wcsapi.BaseHighLevelWCS` and
165 `astropy.wcs.wcsapi.BaseLowLevelWCS` interfaces, but it is not an
166 `astropy.wcs.WCS` (use `fits_wcs` for that).
167 """
168 return cast(ProjectionAstropyView, super().astropy_wcs)
170 @property
171 def psf(self) -> PointSpreadFunction:
172 """The point-spread function model for this image
173 (`.psfs.PointSpreadFunction`).
174 """
175 if isinstance(self._psf, ArchiveReadError):
176 raise self._psf
177 return self._psf
179 def __getitem__(self, bbox: Box) -> VisitImage:
180 result = VisitImage(
181 self.image[bbox],
182 mask=self.mask[bbox],
183 variance=self.variance[bbox],
184 projection=self.projection,
185 psf=self.psf,
186 obs_info=self.obs_info,
187 )
188 if opaque_metadata := getattr(self, "_opaque_metadata", None):
189 result._opaque_metadata = opaque_metadata.subset(bbox)
190 return result
192 def __str__(self) -> str:
193 return f"VisitImage({self.image!s}, {list(self.mask.schema.names)})"
195 def __repr__(self) -> str:
196 return f"VisitImage({self.image!r}, mask_schema={self.mask.schema!r})"
198 def copy(self) -> VisitImage:
199 """Deep-copy the visit image."""
200 result = VisitImage(
201 image=self._image.copy(),
202 mask=self._mask.copy(),
203 variance=self._variance.copy(),
204 psf=self._psf,
205 obs_info=self.obs_info,
206 )
207 if self._opaque_metadata is not None:
208 result._opaque_metadata = self._opaque_metadata.copy()
209 return result
211 def serialize(self, archive: OutputArchive[Any]) -> VisitImageSerializationModel:
212 serialized_image = archive.serialize_direct(
213 "image", functools.partial(self.image.serialize, save_projection=False)
214 )
215 serialized_mask = archive.serialize_direct(
216 "mask", functools.partial(self.mask.serialize, save_projection=False)
217 )
218 serialized_variance = archive.serialize_direct(
219 "variance", functools.partial(self.variance.serialize, save_projection=False)
220 )
221 serialized_projection = (
222 archive.serialize_direct("projection", self.projection.serialize)
223 if self.projection is not None
224 else None
225 )
226 serialized_psf: PiffSerializationModel | PSFExSerializationModel
227 match self._psf:
228 # MyPy is able to figure things out here with this match statement,
229 # but not a single isinstance check on both types.
230 case PiffWrapper():
231 serialized_psf = archive.serialize_direct("psf", self._psf.serialize)
232 case PSFExWrapper():
233 serialized_psf = archive.serialize_direct("psf", self._psf.serialize)
234 case _:
235 raise TypeError(
236 f"Cannot serialize VisitImage with unrecognized PSF type {type(self._psf).__name__}."
237 )
238 return VisitImageSerializationModel(
239 image=serialized_image,
240 mask=serialized_mask,
241 variance=serialized_variance,
242 projection=serialized_projection,
243 psf=serialized_psf,
244 obs_info=self.obs_info,
245 )
247 # Type-checkers want the model argument to only require
248 # MaskedImageSerializationModel[Any], and they'd be absolutely right if
249 # this were a regular instance method. But whether Liskov substitution
250 # applies to classmethods and staticmethods is sort of context-dependent,
251 # and here we do not want it to.
252 @staticmethod
253 def deserialize(
254 model: VisitImageSerializationModel[Any], # type: ignore[override]
255 archive: InputArchive[Any],
256 *,
257 bbox: Box | None = None,
258 ) -> VisitImage:
259 masked_image = MaskedImage.deserialize(model, archive, bbox=bbox)
260 psf = model.deserialize_psf(archive)
261 projection = Projection.deserialize(model.projection, archive)
262 result = VisitImage(
263 masked_image.image,
264 mask=masked_image.mask,
265 variance=masked_image.variance,
266 psf=psf,
267 projection=projection,
268 obs_info=model.obs_info,
269 )
270 return result
272 @staticmethod
273 def _get_archive_tree_type[P: pydantic.BaseModel](
274 pointer_type: type[P],
275 ) -> type[VisitImageSerializationModel[P]]:
276 """Return the serialization model type for this object for an archive
277 type that uses the given pointer type.
278 """
279 return VisitImageSerializationModel[pointer_type] # type: ignore
281 # write_fits and read_fits inherited from MaskedImage.
283 @staticmethod
284 def from_legacy(
285 legacy: Any,
286 *,
287 unit: astropy.units.Unit | None = None,
288 plane_map: Mapping[str, MaskPlane] | None = None,
289 instrument: str | None = None,
290 visit: int | None = None,
291 component: str | None = None,
292 ) -> VisitImage:
293 """Convert from an `lsst.afw.image.Exposure` instance.
295 Parameters
296 ----------
297 legacy
298 An `lsst.afw.image.Exposure` instance that will share image and
299 variance (but not mask) pixel data with the returned object.
300 unit
301 Units of the image. If not provided, the ``BUNIT`` metadata
302 key will be used, if available.
303 plane_map
304 A mapping from legacy mask plane name to the new plane name and
305 description. If `None` (default)
306 `get_legacy_visit_image_mask_planes` is used.
307 instrument
308 Name of the instrument. Extracted from the metadata if not
309 provided.
310 visit
311 ID of the visit. Extracted from the metadata if not provided.
312 """
313 if plane_map is None:
314 plane_map = get_legacy_visit_image_mask_planes()
315 md = legacy.getMetadata()
316 obs_info = _obs_info_from_md(md, visit_info=legacy.info.getVisitInfo())
317 instrument = _extract_or_check_header("LSST BUTLER DATAID INSTRUMENT", instrument, md, str)
318 visit = _extract_or_check_header("LSST BUTLER DATAID VISIT", visit, md, int)
319 unit = _extract_or_check_header("BUNIT", unit, md, lambda x: astropy.units.Unit(x, format="fits"))
320 legacy_wcs = legacy.getWcs()
321 if legacy_wcs is None:
322 raise ValueError("Exposure does not have a SkyWcs.")
323 legacy_detector = legacy.getDetector()
324 if legacy_detector is None:
325 raise ValueError("Exposure does not have a Detector.")
326 detector_bbox = Box.from_legacy(legacy_detector.getBBox())
327 opaque_fits_metadata = FitsOpaqueMetadata()
328 primary_header = astropy.io.fits.Header()
329 with warnings.catch_warnings():
330 # Silence warnings about long keys becoming HIERARCH.
331 warnings.simplefilter("ignore", category=astropy.io.fits.verify.VerifyWarning)
332 primary_header.update(md.toOrderedDict())
333 opaque_fits_metadata.extract_legacy_primary_header(primary_header)
334 projection = Projection.from_legacy(
335 legacy_wcs,
336 DetectorFrame(
337 instrument=instrument,
338 visit=visit,
339 detector=legacy_detector.getId(),
340 bbox=detector_bbox,
341 ),
342 )
343 legacy_psf = legacy.getPsf()
344 if legacy_psf is None:
345 raise ValueError("Exposure file does not have a Psf.")
346 psf = PointSpreadFunction.from_legacy(legacy_psf, bounds=detector_bbox)
347 masked_image = MaskedImage.from_legacy(legacy.getMaskedImage(), unit=unit, plane_map=plane_map)
348 result = VisitImage(
349 image=masked_image.image.view(unit=unit),
350 mask=masked_image.mask,
351 variance=masked_image.variance,
352 projection=projection,
353 psf=psf,
354 obs_info=obs_info,
355 )
356 result._opaque_metadata = masked_image._opaque_metadata
357 return result
359 @overload # type: ignore[override]
360 @staticmethod
361 def read_legacy( 361 ↛ exitline 361 didn't return from function 'read_legacy' because
362 filename: str,
363 *,
364 component: Literal["bbox"],
365 ) -> Box: ...
367 @overload
368 @staticmethod
369 def read_legacy( 369 ↛ exitline 369 didn't return from function 'read_legacy' because
370 filename: str,
371 *,
372 preserve_quantization: bool = False,
373 instrument: str | None = None,
374 visit: int | None = None,
375 component: Literal["image"],
376 ) -> Image: ...
378 @overload
379 @staticmethod
380 def read_legacy( 380 ↛ exitline 380 didn't return from function 'read_legacy' because
381 filename: str,
382 *,
383 plane_map: Mapping[str, MaskPlane] | None = None,
384 instrument: str | None = None,
385 visit: int | None = None,
386 component: Literal["mask"],
387 ) -> Mask: ...
389 @overload
390 @staticmethod
391 def read_legacy( 391 ↛ exitline 391 didn't return from function 'read_legacy' because
392 filename: str,
393 *,
394 preserve_quantization: bool = False,
395 instrument: str | None = None,
396 visit: int | None = None,
397 component: Literal["variance"],
398 ) -> Image: ...
400 @overload
401 @staticmethod
402 def read_legacy( 402 ↛ exitline 402 didn't return from function 'read_legacy' because
403 filename: str,
404 *,
405 instrument: str | None = None,
406 visit: int | None = None,
407 component: Literal["projection"],
408 ) -> Projection[DetectorFrame]: ...
410 @overload
411 @staticmethod
412 def read_legacy( 412 ↛ exitline 412 didn't return from function 'read_legacy' because
413 filename: str,
414 *,
415 component: Literal["psf"],
416 ) -> PointSpreadFunction: ...
418 @overload
419 @staticmethod
420 def read_legacy( 420 ↛ exitline 420 didn't return from function 'read_legacy' because
421 filename: str,
422 *,
423 preserve_quantization: bool = False,
424 plane_map: Mapping[str, MaskPlane] | None = None,
425 instrument: str | None = None,
426 visit: int | None = None,
427 component: None = None,
428 ) -> VisitImage: ...
430 @staticmethod
431 def read_legacy( # type: ignore[override]
432 filename: str,
433 *,
434 preserve_quantization: bool = False,
435 plane_map: Mapping[str, MaskPlane] | None = None,
436 instrument: str | None = None,
437 visit: int | None = None,
438 component: Literal["bbox", "image", "mask", "variance", "projection", "psf"] | None = None,
439 ) -> Any:
440 """Read a FITS file written by `lsst.afw.image.Exposure.writeFits`.
442 Parameters
443 ----------
444 filename
445 Full name of the file.
446 preserve_quantization
447 If `True`, ensure that writing the masked image back out again will
448 exactly preserve quantization-compressed pixel values. This causes
449 the image and variance plane arrays to be marked as read-only and
450 stores the original binary table data for those planes in memory.
451 If the `MaskedImage` is copied, the precompressed pixel values are
452 not transferred to the copy.
453 plane_map
454 A mapping from legacy mask plane name to the new plane name and
455 description. If `None` (default)
456 `get_legacy_visit_image_mask_planes` is used.
457 instrument
458 Name of the instrument. Read from the primary header if not
459 provided.
460 visit
461 ID of the visit. Read from the primary header if not
462 provided.
463 component
464 A component to read instead of the full image.
465 """
466 from lsst.afw.image import ExposureFitsReader
468 reader = ExposureFitsReader(filename)
469 if component == "bbox":
470 return Box.from_legacy(reader.readBBox())
471 legacy_detector = reader.readDetector()
472 if legacy_detector is None:
473 raise ValueError(f"Exposure file {filename!r} does not have a Detector.")
474 detector_bbox = Box.from_legacy(legacy_detector.getBBox())
475 if component in (None, "image", "mask", "variance", "projection"):
476 legacy_wcs = reader.readWcs()
477 if legacy_wcs is None:
478 raise ValueError(f"Exposure file {filename!r} does not have a SkyWcs.")
479 if component in (None, "psf"):
480 legacy_psf = reader.readPsf()
481 if legacy_psf is None:
482 raise ValueError(f"Exposure file {filename!r} does not have a Psf.")
483 psf = PointSpreadFunction.from_legacy(legacy_psf, bounds=detector_bbox)
484 if component == "psf":
485 return psf
486 assert component in (None, "image", "mask", "variance", "projection"), component # for MyPy
487 with astropy.io.fits.open(filename) as hdu_list:
488 primary_header = hdu_list[0].header
489 if instrument is None:
490 try:
491 instrument = primary_header["LSST BUTLER DATAID INSTRUMENT"]
492 except LookupError:
493 raise ValueError(
494 "Instrument could not be found in butler data ID FITS headers and must be provided."
495 ) from None
496 if visit is None:
497 try:
498 visit = primary_header["LSST BUTLER DATAID VISIT"]
499 except LookupError:
500 raise ValueError(
501 "Visit ID could not be found in butler data ID FITS headers and must be provided."
502 ) from None
503 projection = Projection.from_legacy(
504 legacy_wcs,
505 DetectorFrame(
506 instrument=instrument,
507 visit=visit,
508 detector=legacy_detector.getId(),
509 bbox=detector_bbox,
510 ),
511 )
512 if component == "projection":
513 return projection
514 if plane_map is None:
515 plane_map = get_legacy_visit_image_mask_planes()
516 assert component != "psf", component # for MyPy
517 from_masked_image = MaskedImage._read_legacy_hdus(
518 hdu_list,
519 filename,
520 preserve_quantization=preserve_quantization,
521 plane_map=plane_map,
522 component=component,
523 )
524 obs_info = _obs_info_from_md(primary_header)
525 if component is not None:
526 # This is the image, mask, or variance; attach the projection
527 # and return
528 return from_masked_image.view(projection=projection)
529 result = VisitImage(
530 from_masked_image.image,
531 mask=from_masked_image.mask,
532 variance=from_masked_image.variance,
533 projection=projection,
534 psf=psf,
535 obs_info=obs_info,
536 )
537 result._opaque_metadata = from_masked_image._opaque_metadata
538 return result
541class VisitImageSerializationModel[P: pydantic.BaseModel](MaskedImageSerializationModel[P]):
542 """A Pydantic model used to represent a serialized `VisitImage`."""
544 # Inherited attributes are duplicated because that improves the docs
545 # (some limitation in the sphinx/pydantic integration), and these are
546 # important docs.
548 image: ImageSerializationModel[P] = pydantic.Field(description="The main data image.")
549 mask: MaskSerializationModel[P] = pydantic.Field(
550 description="Bitmask that annotates the main image's pixels."
551 )
552 variance: ImageSerializationModel[P] = pydantic.Field(
553 description="Per-pixel variance estimates for the main image."
554 )
555 projection: ProjectionSerializationModel[P] = pydantic.Field(
556 description="Projection that maps the pixel grid to the sky.",
557 )
558 psf: PiffSerializationModel | PSFExSerializationModel | Any = pydantic.Field(
559 union_mode="left_to_right", description="PSF model for the image."
560 )
561 obs_info: ObservationInfo = pydantic.Field(
562 description="Standardized description of visit metadata",
563 )
565 def deserialize_psf(self, archive: InputArchive[Any]) -> PointSpreadFunction | ArchiveReadError:
566 """Finish deserializing the PSF model, or *return* any exception
567 raised in the attempt.
568 """
569 try:
570 match self.psf:
571 case PiffSerializationModel():
572 return PiffWrapper.deserialize(self.psf, archive)
573 case PSFExSerializationModel():
574 return PSFExWrapper.deserialize(self.psf, archive)
575 case _:
576 raise ArchiveReadError("PSF model type not recognized.")
577 except ArchiveReadError as err:
578 return err
581def _extract_or_check_header[T](
582 key: str, given_value: T | None, header: Any, coerce: Callable[[Any], T]
583) -> T:
584 hdr_value: T | None = None
585 if (hdr_raw_value := header.get(key)) is not None:
586 hdr_value = coerce(hdr_raw_value)
587 if given_value is not None:
588 if hdr_value is not None and hdr_value != given_value:
589 raise ValueError(
590 f"Given value {given_value!r} does not match {hdr_value!r} from header key {key}."
591 )
592 return given_value
593 else:
594 if hdr_value is None:
595 raise ValueError(f"No FITS header key {key} found.")
596 return hdr_value