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