Coverage for python/lsst/images/_visit_image.py: 46%
429 statements
« prev ^ index » next coverage.py v7.16.1, created at 2026-09-22 10:09 +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__ = ("VisitImage", "VisitImageSerializationModel")
16import functools
17import logging
18import warnings
19from collections.abc import Callable, Mapping, MutableMapping
20from types import EllipsisType
21from typing import TYPE_CHECKING, Any, ClassVar, Literal, cast
23import astropy.io.fits
24import astropy.units
25import numpy as np
26import pydantic
27from astro_metadata_translator import ObservationInfo, VisitInfoTranslator
29from ._backgrounds import BackgroundMap, BackgroundMapSerializationModel
30from ._concrete_bounds import BoundsSerializationModel
31from ._geom import Bounds, Box
32from ._image import Image, ImageSerializationModel
33from ._mask import Mask, MaskPlane, MaskSchema, MaskSerializationModel, get_legacy_visit_image_mask_planes
34from ._masked_image import MaskedImage, MaskedImageSerializationModel
35from ._observation_summary_stats import ObservationSummaryStats
36from ._polygon import Polygon
37from ._transforms import (
38 DetectorFrame,
39 SkyProjection,
40 SkyProjectionAstropyView,
41 SkyProjectionSerializationModel,
42)
43from .aperture_corrections import (
44 ApertureCorrectionMap,
45 ApertureCorrectionMapSerializationModel,
46 aperture_corrections_from_legacy,
47 aperture_corrections_to_legacy,
48)
49from .cameras import Detector, DetectorSerializationModel
50from .describe import DescribeOptions, FieldRole, Report, ReportField
51from .fields import BaseField, Field, FieldSerializationModel, field_from_legacy_photo_calib
52from .fits import FitsOpaqueMetadata
53from .psfs import (
54 GaussianPointSpreadFunction,
55 GaussianPSFSerializationModel,
56 LegacyPointSpreadFunction,
57 PiffSerializationModel,
58 PiffWrapper,
59 PointSpreadFunction,
60 PSFExSerializationModel,
61 PSFExWrapper,
62)
63from .serialization import ArchiveReadError, InputArchive, InvalidParameterError, MetadataValue, OutputArchive
64from .utils import is_none
66if TYPE_CHECKING:
67 try:
68 from lsst.afw.cameraGeom import Detector as LegacyDetector
69 from lsst.afw.image import Exposure as LegacyExposure
70 from lsst.afw.image import FilterLabel as LegacyFilterLabel
71 from lsst.afw.image import VisitInfo as LegacyVisitInfo
72 except ImportError:
73 type LegacyDetector = Any # type: ignore[no-redef]
74 type LegacyExposure = Any # type: ignore[no-redef]
75 type LegacyFilterLabel = Any # type: ignore[no-redef]
76 type LegacyVisitInfo = Any # type: ignore[no-redef]
78_LOG = logging.getLogger("lsst.images")
81class VisitImage(MaskedImage):
82 """A calibrated single-visit image.
84 Parameters
85 ----------
86 image
87 The main image plane. If this has a `SkyProjection`, it will be used
88 for all planes unless a ``sky_projection`` is passed separately.
89 mask
90 A bitmask image that annotates the main image plane. Must have the
91 same bounding box as ``image`` if provided. Any attached
92 ``sky_projection`` is replaced (possibly by `None`).
93 variance
94 The per-pixel uncertainty of the main image as an image of variance
95 values. Must have the same bounding box as ``image`` if provided, and
96 its units must be the square of ``image.unit`` or `None`.
97 Values default to ``1.0``. Any attached ``sky_projection`` is replaced
98 (possibly by `None`).
99 mask_schema
100 Schema for the mask plane. Must be provided if and only if ``mask`` is
101 not provided.
102 sky_projection
103 Projection that maps the pixel grid to the sky. Can only be `None` if
104 a ``sky_projection`` is already attached to ``image``.
105 bounds
106 The region where this image's pixels and other properties are valid.
107 If not provided, the bounding box of the image is used. Other
108 components (``psf``, ``sky_projection``, ``aperture_corrections``,
109 etc.) are assumed to have their own bounds which may or may not be the
110 same as the image bounds. If ``bounds`` extends beyond the image
111 bounding box, the intersection between ``bounds`` and the image
112 bounding box is used instead.
113 obs_info
114 General information about this visit in standardized form.
115 summary_stats
116 Summary statistics associated with this visit. Initialized to default
117 values if not provided.
118 photometric_scaling
119 Field that can be used to multiply a post-ISR image units to yield
120 calibrated image units. This may be a scaling that was already
121 applied (so dividing by it will recover the post-ISR units) or a
122 scaling that has not been applied, depending on ``image.unit``.
123 psf
124 Point-spread function model for this image, or an exception explaining
125 why it could not be read (to be raised if the PSF is requested later).
126 detector
127 Geometry and electronic information for the detector attached to this
128 image.
129 aperture_corrections : `dict` [`str`, `~fields.BaseField`]
130 Mapping from photometry algorithm name to the aperture correction for
131 that algorithm.
132 backgrounds
133 Background models associated with this image.
134 band
135 Name of the passband the image was observed with (this is a shorter,
136 less specific version of ``obs_info.physical_filter``).
137 metadata
138 Arbitrary flexible metadata to associate with the image.
139 """
141 def __init__(
142 self,
143 image: Image,
144 *,
145 mask: Mask | None = None,
146 variance: Image | None = None,
147 mask_schema: MaskSchema | None = None,
148 sky_projection: SkyProjection[DetectorFrame] | None = None,
149 bounds: Bounds | None = None,
150 obs_info: ObservationInfo | None = None,
151 summary_stats: ObservationSummaryStats | None = None,
152 photometric_scaling: Field | None = None,
153 psf: PointSpreadFunction | ArchiveReadError,
154 detector: Detector,
155 aperture_corrections: ApertureCorrectionMap | None = None,
156 backgrounds: BackgroundMap | None = None,
157 band: str,
158 metadata: dict[str, MetadataValue] | None = None,
159 ) -> None:
160 super().__init__(
161 image,
162 mask=mask,
163 variance=variance,
164 mask_schema=mask_schema,
165 sky_projection=sky_projection,
166 metadata=metadata,
167 )
168 if self.image.unit is None:
169 raise TypeError("The image component of a VisitImage must have units.")
170 if self.image.sky_projection is None:
171 raise TypeError("The sky_projection component of a VisitImage cannot be None.")
172 if obs_info is None:
173 raise TypeError("The observation info component of a VisitImage cannot be None.")
174 if obs_info.physical_filter is None: 174 ↛ 175line 174 didn't jump to line 175 because the condition on line 174 was never true
175 raise ValueError("The obs_info.physical_filter attribute of a VisitImage cannot be None.")
176 self._obs_info = obs_info
177 if not isinstance(self.image.sky_projection.pixel_frame, DetectorFrame):
178 raise TypeError("The sky_projection's pixel frame must be a DetectorFrame for VisitImage.")
179 if summary_stats is None:
180 summary_stats = ObservationSummaryStats()
181 self._summary_stats = summary_stats
182 if photometric_scaling is not None and photometric_scaling.unit is None: 182 ↛ 183line 182 didn't jump to line 183 because the condition on line 182 was never true
183 raise TypeError("If a photometric_scaling is provided, it must have units.")
184 self._photometric_scaling = photometric_scaling
185 self._psf = psf
186 self._detector = detector
187 self._aperture_corrections = aperture_corrections if aperture_corrections is not None else {}
188 self._bounds = bounds if bounds is not None else self.bbox
189 if not self.bbox.contains(self._bounds.bbox):
190 self._bounds = self._bounds.intersection(self.bbox)
191 self._backgrounds = backgrounds if backgrounds is not None else BackgroundMap()
192 self._band = band
194 @property
195 def unit(self) -> astropy.units.UnitBase:
196 """The units of the image plane (`astropy.units.Unit`)."""
197 return cast(astropy.units.UnitBase, super().unit)
199 @property
200 def sky_projection(self) -> SkyProjection[DetectorFrame]:
201 """The projection that maps the pixel grid to the sky
202 (`SkyProjection` [`DetectorFrame`]).
203 """
204 return cast(SkyProjection[DetectorFrame], super().sky_projection)
206 @property
207 def bounds(self) -> Bounds:
208 """The region where pixels are valid (`Bounds`)."""
209 return self._bounds
211 @property
212 def obs_info(self) -> ObservationInfo:
213 """General information about this observation in standard form.
214 (`~astro_metadata_translator.ObservationInfo`).
215 """
216 return self._obs_info
218 @property
219 def physical_filter(self) -> str:
220 """Full name of the physical bandpass filter (`str`)."""
221 assert self._obs_info.physical_filter is not None, "Guaranteed at construction."
222 return self._obs_info.physical_filter
224 @property
225 def band(self) -> str:
226 """Short name of the bandpass filter (`str`)."""
227 return self._band
229 @property
230 def astropy_wcs(self) -> SkyProjectionAstropyView:
231 """An Astropy WCS for the pixel arrays (`SkyProjectionAstropyView`).
233 Notes
234 -----
235 As expected for Astropy WCS objects, this defines pixel coordinates
236 such that the first row and column in the arrays are ``(0, 0)``, not
237 ``bbox.start``, as is the case for `sky_projection`.
239 This object satisfies the `astropy.wcs.wcsapi.BaseHighLevelWCS` and
240 `astropy.wcs.wcsapi.BaseLowLevelWCS` interfaces, but it is not an
241 `astropy.wcs.WCS` (use `fits_wcs` for that).
242 """
243 return cast(SkyProjectionAstropyView, super().astropy_wcs)
245 @property
246 def summary_stats(self) -> ObservationSummaryStats:
247 """Optional summary statistics for this observation
248 (`ObservationSummaryStats`).
249 """
250 return self._summary_stats
252 @property
253 def photometric_scaling(self) -> Field | None:
254 """Field that multiplies a post-ISR image to yield the calibrated
255 image (~`fields.BaseField`).
256 """
257 return self._photometric_scaling
259 @photometric_scaling.setter
260 def photometric_scaling(self, value: Field) -> None:
261 if value.unit is None: 261 ↛ 262line 261 didn't jump to line 262 because the condition on line 261 was never true
262 raise TypeError("The photometric_scaling for a VisitImage must have units.")
263 self._photometric_scaling = value
265 @property
266 def psf(self) -> PointSpreadFunction:
267 """The point-spread function model for this image
268 (`.psfs.PointSpreadFunction`).
269 """
270 if isinstance(self._psf, ArchiveReadError):
271 raise self._psf
272 return self._psf
274 @property
275 def detector(self) -> Detector:
276 """Geometry and electronic information about the detector
277 (`.cameras.Detector`).
278 """
279 return self._detector
281 @property
282 def aperture_corrections(self) -> ApertureCorrectionMap:
283 """A mapping from photometry algorithm name to the aperture correction
284 field for that algorithm (`dict` [`str`, `~.fields.BaseField`]).
285 """
286 return self._aperture_corrections
288 @property
289 def backgrounds(self) -> BackgroundMap:
290 """A mapping of backgrounds associated with this image
291 (`BackgroundMap`).
292 """
293 return self._backgrounds
295 def __getitem__(self, bbox: Box | EllipsisType) -> VisitImage:
296 bbox, _ = self._handle_getitem_args(bbox)
297 return self._transfer_metadata(
298 VisitImage(
299 self.image[bbox],
300 mask=self.mask[bbox],
301 variance=self.variance[bbox],
302 sky_projection=self.sky_projection,
303 psf=self._psf,
304 obs_info=self.obs_info,
305 bounds=self._bounds, # don't need to intersect here, because __init__ will do that.
306 summary_stats=self.summary_stats,
307 detector=self._detector,
308 photometric_scaling=self._photometric_scaling,
309 aperture_corrections=self.aperture_corrections,
310 backgrounds=self._backgrounds,
311 band=self._band,
312 ),
313 bbox=bbox,
314 )
316 def _describe(self, options: DescribeOptions = DescribeOptions(), /) -> Report:
317 """Return a `Report` describing this visit image.
319 Parameters
320 ----------
321 options : `DescribeOptions`, optional
322 Rendering options; forwarded to all children. Child construction
323 can be expensive or raise for an unreadable component, so
324 `DescribeOptions.brief` skips it.
325 """
326 # The image and mask schema are rendered as children below, and the
327 # image renders as ``Image(bbox, dtype)``, which would restate the
328 # shared bbox. Both are REPR_ONLY so only repr sees them.
329 fields = [
330 ReportField(
331 label="image",
332 value=self.image,
333 repr_value=repr(self.image),
334 positional=True,
335 role=FieldRole.REPR_ONLY,
336 ),
337 ReportField(
338 label="mask_schema",
339 value=self.mask.schema,
340 repr_value=repr(self.mask.schema),
341 role=FieldRole.REPR_ONLY,
342 ),
343 ReportField(label="band", value=self.band, role=FieldRole.DERIVED),
344 ReportField(label="physical_filter", value=self.physical_filter, role=FieldRole.DERIVED),
345 ReportField(label="bbox", value=self.bbox, repr_value=repr(self.bbox), role=FieldRole.DERIVED),
346 ]
347 summary = f"VisitImage({self.image!s}, {list(self.mask.schema.names)})"
348 if options.brief:
349 return Report(type_name="VisitImage", summary=summary, fields=fields)
350 child = options.for_child()
351 plane = options.for_child("sky_projection", "bbox")
352 children: dict[str, Report] = {
353 "image": self.image._describe(plane),
354 "mask": self.mask._describe(plane),
355 "variance": self.variance._describe(plane),
356 "sky_projection": self.sky_projection._describe(child, bbox=self.bbox),
357 "psf": self.psf._describe(child),
358 "detector": self.detector._describe(child),
359 "summary_stats": self.summary_stats._describe(child),
360 "backgrounds": self.backgrounds._describe(child),
361 }
362 if self.photometric_scaling is not None: 362 ↛ 363line 362 didn't jump to line 363 because the condition on line 362 was never true
363 children["photometric_scaling"] = self.photometric_scaling._describe(child)
364 return Report(
365 type_name="VisitImage",
366 summary=summary,
367 fields=fields,
368 children=children,
369 )
371 def copy(self, *, copy_detector: bool = False) -> VisitImage:
372 """Deep-copy the visit image.
374 Parameters
375 ----------
376 copy_detector
377 Whether to deep-copy the `detector` attribute.
378 """
379 return self._transfer_metadata(
380 VisitImage(
381 image=self._image.copy(),
382 mask=self._mask.copy(),
383 variance=self._variance.copy(),
384 psf=self._psf,
385 obs_info=self.obs_info,
386 bounds=self._bounds,
387 summary_stats=self.summary_stats.model_copy(),
388 detector=self._detector.copy() if copy_detector else self._detector,
389 photometric_scaling=self._photometric_scaling,
390 aperture_corrections=self.aperture_corrections.copy(),
391 backgrounds=self._backgrounds.copy(),
392 band=self.band,
393 ),
394 copy=True,
395 )
397 def convert_unit(
398 self,
399 unit: astropy.units.UnitBase = astropy.units.nJy,
400 copy: Literal["as-needed"] | bool = True,
401 copy_detector: bool = False,
402 ) -> VisitImage:
403 """Return an equivalent image with different pixel units.
405 Parameters
406 ----------
407 unit
408 The unit to transform to. This may be any of the following:
410 - any unit directly relatable to the current units via Astropy;
411 - any unit relatable to the product of the current units with the
412 `photometric_scaling` (i.e. if the current image is in
413 instrumental units but we know how to calibrate them)
414 - any unit relatable to the quotient of the current units with the
415 `photometric_scaling` (i.e. if the current image is in
416 calibrated units and we want to revert back to instrumental
417 units).
418 copy
419 Whether to copy the images and other components. If `True`, all
420 components that aren't controlled by some other argument will
421 always be deep-copied. If `False`, the operation will fail if the
422 image is not already in the right units. If ``as-needed``, only
423 the image and variance will be copied, and only if they are not
424 already in the right units.
425 copy_detector
426 Whether to deep-copy the `detector` attribute.
428 Returns
429 -------
430 `VisitImage`
431 An image with the given units.
432 """
433 if copy not in (True, False, "as-needed"): 433 ↛ 434line 433 didn't jump to line 434 because the condition on line 433 was never true
434 raise TypeError(f"Invalid value for 'copy' parameter: {copy!r}.")
435 if (factor := _get_unit_conversion_factor(self.unit, unit)) is not None: 435 ↛ 436line 435 didn't jump to line 436 because the condition on line 435 was never true
436 if factor == 1.0:
437 if copy is True: # not "as-needed"
438 return self.copy()
439 else:
440 return self[...]
441 elif copy is False:
442 raise astropy.units.UnitConversionError(
443 f"Units must be converted ({self.unit} -> {unit}), but copy=False."
444 )
445 image = Image(
446 self._image.array * factor, bbox=self.bbox, sky_projection=self.sky_projection, unit=unit
447 )
448 variance = Image(
449 self._variance.array * factor**2,
450 bbox=self.bbox,
451 unit=unit**2,
452 )
453 elif self._photometric_scaling is None: 453 ↛ 454line 453 didn't jump to line 454 because the condition on line 453 was never true
454 raise astropy.units.UnitConversionError(
455 "VisitImage.photometric_scaling is None, and there "
456 f"is no constant conversion from {self.unit} to {unit}."
457 )
458 else:
459 if copy is False: 459 ↛ 460line 459 didn't jump to line 460 because the condition on line 459 was never true
460 raise astropy.units.UnitConversionError(
461 f"Photometric scaling must be applied to go from ={self.unit} to {unit}, but copy=False."
462 )
463 scaling = self._photometric_scaling
464 assert scaling.unit is not None, "Checked at construction."
465 if (constant_factor := _get_unit_conversion_factor(self.unit * scaling.unit, unit)) is not None:
466 if constant_factor != 1.0: 466 ↛ 467line 466 didn't jump to line 467 because the condition on line 466 was never true
467 scaling = scaling * constant_factor
468 scaling_array = scaling.render(self.bbox, dtype=self.image.array.dtype).array
469 elif (constant_factor := _get_unit_conversion_factor(self.unit / scaling.unit, unit)) is not None: 469 ↛ 475line 469 didn't jump to line 475 because the condition on line 469 was always true
470 if constant_factor != 1.0: 470 ↛ 471line 470 didn't jump to line 471 because the condition on line 470 was never true
471 scaling = scaling / constant_factor
472 scaling_array = scaling.render(self.bbox, dtype=self.image.array.dtype).array
473 np.true_divide(1.0, scaling_array, out=scaling_array)
474 else:
475 raise astropy.units.UnitConversionError(
476 f"photometric_scaling with units {scaling.unit} does not "
477 f"provide a path from {self.unit} to {unit}."
478 )
479 # We needed to allocate a new array to evaluate the scaling field,
480 # and then we need to allocate another to hold its square for the
481 # variance scaling. But then we can multiply those arrays in-place
482 # to get the output image and variance to avoid yet more
483 # allocations (note we can't instead multiply the visit image's
484 # image and variance arrays in place because they might have other
485 # references that are still associated with the old units).
486 image = Image(scaling_array, bbox=self.bbox, unit=unit)
487 variance = Image(np.square(scaling_array), bbox=self.bbox, unit=unit**2)
488 image.array *= self._image.array
489 variance.array *= self._variance.array
490 copy_components = copy is True
491 return self._transfer_metadata(
492 VisitImage(
493 image=image,
494 mask=self._mask if not copy_components else self._mask.copy(),
495 variance=variance,
496 sky_projection=self.sky_projection, # never copied; immutable
497 obs_info=self.obs_info if not copy_components else self.obs_info.model_copy(),
498 psf=self._psf, # never copied; immutable
499 bounds=self._bounds, # never copied; immutable
500 summary_stats=self.summary_stats if not copy_components else self.summary_stats.model_copy(),
501 detector=self._detector if not copy_detector else self._detector.copy(),
502 photometric_scaling=self._photometric_scaling, # never copied; immutable
503 aperture_corrections=(
504 self.aperture_corrections if not copy_components else self.aperture_corrections.copy()
505 ),
506 backgrounds=self.backgrounds if not copy_components else self.backgrounds.copy(),
507 band=self.band,
508 )
509 )
511 def serialize(self, archive: OutputArchive[Any]) -> VisitImageSerializationModel[Any]:
512 return self._serialize_impl(VisitImageSerializationModel, archive)
514 # This is slightly bad Liskov substitution - we're demanding M be a
515 # VisitImageSerializationModel, not just a MaskedImageSerializationModel,
516 # but that's because we know only `serialize` will call it.
517 def _serialize_impl[M: VisitImageSerializationModel[Any]]( # type: ignore[override]
518 self, model_type: type[M], archive: OutputArchive[Any]
519 ) -> M:
520 result = super()._serialize_impl(model_type, archive)
521 match self._psf:
522 # MyPy is able to figure things out here with this match statement,
523 # but not a single isinstance check on the three types.
524 case PiffWrapper(): 524 ↛ 525line 524 didn't jump to line 525 because the pattern on line 524 never matched
525 result.psf = archive.serialize_direct("psf", self._psf.serialize)
526 case PSFExWrapper(): 526 ↛ 527line 526 didn't jump to line 527 because the pattern on line 526 never matched
527 result.psf = archive.serialize_direct("psf", self._psf.serialize)
528 case GaussianPointSpreadFunction(): 528 ↛ 530line 528 didn't jump to line 530 because the pattern on line 528 always matched
529 result.psf = archive.serialize_direct("psf", self._psf.serialize)
530 case _:
531 raise TypeError(
532 f"Cannot serialize VisitImage with unrecognized PSF type {type(self._psf).__name__}."
533 )
534 assert result.sky_projection is not None, "VisitImage always has a sky_projection."
535 result.obs_info = self.obs_info
536 result.summary_stats = self.summary_stats
537 result.bounds = self._bounds.serialize() if self._bounds != self.bbox else None
538 result.detector = archive.serialize_direct("detector", self._detector.serialize)
539 result.band = self.band
540 result.photometric_scaling = (
541 # MyPy can't quite follow the type union through the serialize
542 # method return types.
543 archive.serialize_direct(
544 "photometric_scaling",
545 self._photometric_scaling.serialize,
546 ) # type: ignore[assignment]
547 if self._photometric_scaling is not None
548 else None
549 )
550 result.aperture_corrections = archive.serialize_direct(
551 "aperture_corrections",
552 functools.partial(ApertureCorrectionMapSerializationModel.serialize, self.aperture_corrections),
553 )
554 result.backgrounds = archive.serialize_direct("backgrounds", self._backgrounds.serialize)
555 return result
557 @staticmethod
558 def _get_archive_tree_type[P: pydantic.BaseModel](
559 pointer_type: type[P],
560 ) -> type[VisitImageSerializationModel[P]]:
561 """Return the serialization model type for this object for an archive
562 type that uses the given pointer type.
563 """
564 return VisitImageSerializationModel[pointer_type] # type: ignore
566 @staticmethod
567 def from_legacy( # type: ignore[override]
568 legacy: LegacyExposure,
569 *,
570 unit: astropy.units.UnitBase | None = None,
571 plane_map: Mapping[str, MaskPlane] | None = None,
572 instrument: str | None = None,
573 visit: int | None = None,
574 ) -> VisitImage:
575 """Convert from an `lsst.afw.image.Exposure` instance.
577 Parameters
578 ----------
579 legacy
580 An `lsst.afw.image.Exposure` instance that will share image and
581 variance (but not mask) pixel data with the returned object.
582 unit
583 Units of the image. If not provided, the ``BUNIT`` metadata
584 key will be used, if available.
585 plane_map
586 A mapping from legacy mask plane name to the new plane name and
587 description. If `None` (default)
588 `get_legacy_visit_image_mask_planes` is used.
589 instrument
590 Name of the instrument. Extracted from the metadata if not
591 provided.
592 visit
593 ID of the visit. Extracted from the metadata if not provided.
594 """
595 if plane_map is None: 595 ↛ 597line 595 didn't jump to line 597 because the condition on line 595 was always true
596 plane_map = get_legacy_visit_image_mask_planes()
597 md = legacy.getMetadata()
598 obs_info = _obs_info_from_md(md, visit_info=legacy.info.getVisitInfo())
599 instrument = _extract_or_check_header(
600 "LSST BUTLER DATAID INSTRUMENT", instrument, md, obs_info.instrument, str
601 )
602 visit = _extract_or_check_header("LSST BUTLER DATAID VISIT", visit, md, obs_info.exposure_id, int)
603 legacy_wcs = legacy.getWcs()
604 if legacy_wcs is None: 604 ↛ 605line 604 didn't jump to line 605 because the condition on line 604 was never true
605 raise ValueError("Exposure does not have a SkyWcs.")
606 legacy_detector = legacy.getDetector()
607 if legacy_detector is None: 607 ↛ 608line 607 didn't jump to line 608 because the condition on line 607 was never true
608 raise ValueError("Exposure does not have a Detector.")
609 detector_bbox = Box.from_legacy(legacy_detector.getBBox())
611 # Update the ObservationInfo from other components.
612 obs_info = _update_obs_info_from_legacy(obs_info, legacy_detector, legacy.info.getFilter())
614 opaque_fits_metadata = FitsOpaqueMetadata()
615 primary_header = astropy.io.fits.Header()
616 with warnings.catch_warnings():
617 # Silence warnings about long keys becoming HIERARCH.
618 warnings.simplefilter("ignore", category=astropy.io.fits.verify.VerifyWarning)
619 for name in md.getOrderedNames():
620 # Some keys may be set more than once.
621 # Write one card per value in those cases.
622 for value in md.getArray(name):
623 primary_header.append((name, value), end=True)
624 metadata = opaque_fits_metadata.extract_legacy_primary_header(primary_header)
625 instrumental_unit = opaque_fits_metadata.get_instrumental_unit() or astropy.units.electron
626 hdr_unit: astropy.units.UnitBase | None = None
627 if hdr_unit_str := md.get("BUNIT"): 627 ↛ 633line 627 didn't jump to line 633 because the condition on line 627 was always true
628 hdr_unit = astropy.units.Unit(hdr_unit_str, format="FITS")
629 if hdr_unit == astropy.units.adu and instrumental_unit == astropy.units.electron: 629 ↛ 632line 629 didn't jump to line 632 because the condition on line 629 was never true
630 # Fix incorrect BUNIT='adu' in LSST
631 # preliminary_visit_image.
632 hdr_unit = astropy.units.electron
633 if unit is None: 633 ↛ 635line 633 didn't jump to line 635 because the condition on line 633 was always true
634 unit = hdr_unit
635 elif hdr_unit is not None and hdr_unit != unit:
636 raise ValueError(f"BUNIT value {hdr_unit} disagrees with given unit {unit}.")
637 sky_projection = SkyProjection.from_legacy(
638 legacy_wcs,
639 DetectorFrame(
640 instrument=instrument,
641 visit=visit,
642 detector=legacy_detector.getId(),
643 bbox=detector_bbox,
644 ),
645 )
646 legacy_psf = legacy.getPsf()
647 if legacy_psf is None: 647 ↛ 648line 647 didn't jump to line 648 because the condition on line 647 was never true
648 raise ValueError("Exposure file does not have a Psf.")
649 psf = PointSpreadFunction.from_legacy(legacy_psf, bounds=detector_bbox)
650 masked_image = MaskedImage.from_legacy(legacy.getMaskedImage(), unit=unit, plane_map=plane_map)
651 legacy_summary_stats = legacy.info.getSummaryStats()
652 legacy_ap_corr_map = legacy.info.getApCorrMap()
653 legacy_polygon = legacy.info.getValidPolygon()
654 legacy_photo_calib = legacy.info.getPhotoCalib()
655 detector = Detector.from_legacy(
656 legacy_detector, instrument=instrument, visit=visit, is_raw_assembled=True
657 )
658 _reconcile_detector_serial(obs_info, detector)
659 result = VisitImage(
660 image=masked_image.image.view(unit=unit),
661 mask=masked_image.mask,
662 variance=masked_image.variance,
663 sky_projection=sky_projection,
664 psf=psf,
665 obs_info=obs_info,
666 summary_stats=(
667 ObservationSummaryStats.from_legacy(legacy_summary_stats)
668 if legacy_summary_stats is not None
669 else None
670 ),
671 detector=detector,
672 aperture_corrections=(
673 aperture_corrections_from_legacy(legacy_ap_corr_map)
674 if legacy_ap_corr_map is not None
675 else None
676 ),
677 bounds=Polygon.from_legacy(legacy_polygon) if legacy_polygon is not None else None,
678 photometric_scaling=(
679 field_from_legacy_photo_calib(
680 legacy_photo_calib, bounds=detector_bbox, instrumental_unit=instrumental_unit
681 )
682 if legacy_photo_calib is not None
683 else None
684 ),
685 band=legacy.info.getFilter().bandLabel,
686 metadata=metadata,
687 )
688 result.metadata["id"] = legacy.info.getId()
689 result._opaque_metadata = opaque_fits_metadata
690 return result
692 def to_legacy(
693 self, *, copy: bool | None = None, plane_map: Mapping[str, MaskPlane] | None = None
694 ) -> LegacyExposure:
695 """Convert to an `lsst.afw.image.Exposure` instance.
697 Parameters
698 ----------
699 copy
700 If `True`, always copy the image and variance pixel data.
701 If `False`, return a view, and raise `TypeError` if the pixel data
702 is read-only (this is not supported by afw). If `None`, only copy
703 if the pixel data is read-only. Mask pixel data is always copied.
704 plane_map
705 A mapping from legacy mask plane name to the new plane name and
706 description. If `None` (default),
707 `get_legacy_visit_image_mask_planes` is used.
708 """
709 from lsst.afw.image import Exposure as LegacyExposure
710 from lsst.afw.image import FilterLabel as LegacyFilterLabel
711 from lsst.obs.base.makeRawVisitInfoViaObsInfo import MakeRawVisitInfoViaObsInfo
713 if plane_map is None: 713 ↛ 715line 713 didn't jump to line 715 because the condition on line 713 was always true
714 plane_map = get_legacy_visit_image_mask_planes()
715 legacy_masked_image = super().to_legacy(copy=copy, plane_map=plane_map)
716 result = LegacyExposure(legacy_masked_image, dtype=self.image.array.dtype)
717 result_info = result.info
718 result_info.setId(self.metadata.get("id"))
719 result_info.setWcs(self.sky_projection.to_legacy())
720 result_info.setDetector(self.detector.to_legacy())
721 result_info.setFilter(LegacyFilterLabel.fromBandPhysical(self.band, self.obs_info.physical_filter))
722 if self._photometric_scaling is not None: 722 ↛ 723line 722 didn't jump to line 723 because the condition on line 722 was never true
723 result_info.setPhotoCalib(self._photometric_scaling.to_legacy_photo_calib(self.unit))
724 else:
725 result_info.setPhotoCalib(BaseField.make_legacy_photo_calib(self.unit))
726 self._fill_legacy_metadata(result_info.getMetadata())
727 if isinstance(self._psf, LegacyPointSpreadFunction): 727 ↛ 729line 727 didn't jump to line 729 because the condition on line 727 was always true
728 result_info.setPsf(self._psf.legacy_psf)
729 elif isinstance(self._psf, PiffWrapper):
730 result_info.setPsf(self._psf.to_legacy())
731 if isinstance(self.bounds, Polygon): 731 ↛ 732line 731 didn't jump to line 732 because the condition on line 731 was never true
732 result_info.setValidPolygon(self.bounds.to_legacy())
733 if self.aperture_corrections: 733 ↛ 734line 733 didn't jump to line 734 because the condition on line 733 was never true
734 result_info.setApCorrMap(aperture_corrections_to_legacy(self.aperture_corrections))
735 result_info.setVisitInfo(MakeRawVisitInfoViaObsInfo.observationInfo2visitInfo(self.obs_info))
736 result_info.setSummaryStats(self.summary_stats.to_legacy())
737 return result
739 @staticmethod
740 def read_legacy( # type: ignore[override]
741 filename: str,
742 *,
743 preserve_quantization: bool = False,
744 plane_map: Mapping[str, MaskPlane] | None = None,
745 instrument: str | None = None,
746 visit: int | None = None,
747 component: Literal[
748 "bbox",
749 "image",
750 "mask",
751 "variance",
752 "sky_projection",
753 "psf",
754 "detector",
755 "photometric_scaling",
756 "obs_info",
757 "summary_stats",
758 "aperture_corrections",
759 ]
760 | None = None,
761 ) -> Any:
762 """Read a FITS file written by `lsst.afw.image.Exposure.writeFits`.
764 Parameters
765 ----------
766 filename
767 Full name of the file.
768 preserve_quantization
769 If `True`, ensure that writing the masked image back out again will
770 exactly preserve quantization-compressed pixel values. This causes
771 the image and variance plane arrays to be marked as read-only and
772 stores the original binary table data for those planes in memory.
773 If the `MaskedImage` is copied, the precompressed pixel values are
774 not transferred to the copy.
775 plane_map
776 A mapping from legacy mask plane name to the new plane name and
777 description. If `None` (default)
778 `get_legacy_visit_image_mask_planes` is used.
779 instrument
780 Name of the instrument. Read from the primary header if not
781 provided.
782 visit
783 ID of the visit. Read from the primary header if not
784 provided.
785 component
786 A component to read instead of the full image.
787 """
788 from lsst.afw.image import ExposureFitsReader
790 reader = ExposureFitsReader(filename)
791 if component == "bbox":
792 return Box.from_legacy(reader.readBBox())
793 legacy_detector = reader.readDetector()
794 if legacy_detector is None:
795 raise ValueError(f"Exposure file {filename!r} does not have a Detector.")
796 detector_bbox = Box.from_legacy(legacy_detector.getBBox())
797 legacy_wcs = None
798 if component in (None, "image", "mask", "variance", "sky_projection"):
799 legacy_wcs = reader.readWcs()
800 if legacy_wcs is None:
801 raise ValueError(f"Exposure file {filename!r} does not have a SkyWcs.")
802 legacy_exposure_info = reader.readExposureInfo()
803 summary_stats = None
804 if component in (None, "summary_stats"):
805 legacy_stats = legacy_exposure_info.getSummaryStats()
806 if legacy_stats is not None:
807 summary_stats = ObservationSummaryStats.from_legacy(legacy_stats)
808 if component == "summary_stats":
809 return summary_stats
810 if component in (None, "psf"):
811 legacy_psf = reader.readPsf()
812 if legacy_psf is None:
813 raise ValueError(f"Exposure file {filename!r} does not have a Psf.")
814 psf = PointSpreadFunction.from_legacy(legacy_psf, bounds=detector_bbox)
815 if component == "psf":
816 return psf
817 aperture_corrections: ApertureCorrectionMap = {}
818 if component in (None, "aperture_corrections"):
819 legacy_ap_corr_map = reader.readApCorrMap()
820 if legacy_ap_corr_map is not None:
821 aperture_corrections = aperture_corrections_from_legacy(legacy_ap_corr_map)
822 if component == "aperture_corrections":
823 return aperture_corrections
824 assert component in (
825 None,
826 "image",
827 "mask",
828 "variance",
829 "sky_projection",
830 "obs_info",
831 "detector",
832 "photometric_scaling",
833 ), component # for MyPy
834 filter_label = reader.readFilter()
835 with astropy.io.fits.open(filename) as hdu_list:
836 primary_header = hdu_list[0].header
837 obs_info = _obs_info_from_md(primary_header)
838 obs_info = _update_obs_info_from_legacy(obs_info, legacy_detector, filter_label)
839 if component == "obs_info":
840 return obs_info
841 instrument = _extract_or_check_header(
842 "LSST BUTLER DATAID INSTRUMENT", instrument, primary_header, obs_info.instrument, str
843 )
844 visit = _extract_or_check_header(
845 "LSST BUTLER DATAID VISIT", visit, primary_header, obs_info.exposure_id, int
846 )
847 opaque_metadata = FitsOpaqueMetadata()
848 # This extraction is destructive, so we need to be sure to pass
849 # this opaque_metadata down to MaskedImage._read_legacy_hdus
850 # so it doesn't try to extract it again.
851 metadata = opaque_metadata.extract_legacy_primary_header(primary_header)
852 if (instrumental_unit := opaque_metadata.get_instrumental_unit()) is None:
853 instrumental_unit = astropy.units.electron
854 photometric_scaling: Field | None = None
855 if component in (None, "photometric_scaling"):
856 legacy_photo_calib = reader.readPhotoCalib()
857 if legacy_photo_calib is not None:
858 photometric_scaling = field_from_legacy_photo_calib(
859 legacy_photo_calib, bounds=detector_bbox, instrumental_unit=instrumental_unit
860 )
861 if component == "photometric_scaling":
862 return photometric_scaling
863 if component in ("detector", None):
864 detector = Detector.from_legacy(
865 legacy_detector, instrument=instrument, visit=visit, is_raw_assembled=True
866 )
867 _reconcile_detector_serial(obs_info, detector)
868 if component == "detector":
869 return detector
870 assert component != "detector", "MyPy can't work this out from the above."
871 sky_projection = SkyProjection.from_legacy(
872 legacy_wcs,
873 DetectorFrame(
874 instrument=instrument,
875 visit=visit,
876 detector=legacy_detector.getId(),
877 bbox=detector_bbox,
878 ),
879 )
880 if component == "sky_projection":
881 return sky_projection
882 if plane_map is None:
883 plane_map = get_legacy_visit_image_mask_planes()
884 from_masked_image = MaskedImage._read_legacy_hdus(
885 hdu_list,
886 filename,
887 opaque_metadata=opaque_metadata,
888 preserve_quantization=preserve_quantization,
889 plane_map=plane_map,
890 component=component,
891 )
892 if component is not None:
893 # This is the image, mask, or variance; attach the sky_projection
894 # and obs_info and return
895 return from_masked_image.view(sky_projection=sky_projection)
896 legacy_polygon = reader.readValidPolygon()
897 result = VisitImage(
898 from_masked_image.image,
899 mask=from_masked_image.mask,
900 variance=from_masked_image.variance,
901 sky_projection=sky_projection,
902 psf=psf,
903 detector=detector,
904 obs_info=obs_info,
905 summary_stats=summary_stats,
906 aperture_corrections=aperture_corrections,
907 bounds=Polygon.from_legacy(legacy_polygon) if legacy_polygon is not None else None,
908 photometric_scaling=photometric_scaling,
909 band=filter_label.bandLabel,
910 metadata=metadata,
911 )
912 result._opaque_metadata = from_masked_image._opaque_metadata
913 result.metadata["id"] = reader.readExposureId()
914 return result
917class VisitImageSerializationModel[P: pydantic.BaseModel](MaskedImageSerializationModel[P]):
918 """A Pydantic model used to represent a serialized `VisitImage`."""
920 SCHEMA_NAME: ClassVar[str] = "visit_image"
921 SCHEMA_VERSION: ClassVar[str] = "1.0.0"
922 MIN_READ_VERSION: ClassVar[int] = 1
923 PUBLIC_TYPE: ClassVar[type] = VisitImage
925 # Inherited attributes are duplicated because that improves the docs
926 # (some limitation in the sphinx/pydantic integration), and these are
927 # important docs.
929 image: ImageSerializationModel[P] = pydantic.Field(description="The main data image.")
930 mask: MaskSerializationModel[P] = pydantic.Field(
931 description="Bitmask that annotates the main image's pixels."
932 )
933 variance: ImageSerializationModel[P] = pydantic.Field(
934 description="Per-pixel variance estimates for the main image."
935 )
936 sky_projection: SkyProjectionSerializationModel[P] = pydantic.Field(
937 description="Projection that maps the pixel grid to the sky.",
938 )
939 psf: PiffSerializationModel | PSFExSerializationModel | GaussianPSFSerializationModel | Any = (
940 pydantic.Field(union_mode="left_to_right", description="PSF model for the image.")
941 )
942 obs_info: ObservationInfo = pydantic.Field(
943 description="Standardized description of visit metadata",
944 )
945 photometric_scaling: FieldSerializationModel | None = pydantic.Field(
946 default=None,
947 description="Scaling that can be used to multiply a post-ISR image to yield calibrated pixel values.",
948 )
949 summary_stats: ObservationSummaryStats = pydantic.Field(
950 description="Summary statistics for the observation."
951 )
952 detector: DetectorSerializationModel = pydantic.Field(
953 description="Geometry and electronic information for the detector."
954 )
955 aperture_corrections: ApertureCorrectionMapSerializationModel = pydantic.Field(
956 default_factory=ApertureCorrectionMapSerializationModel,
957 description="Aperture corrections, keyed by flux algorithm.",
958 )
959 bounds: BoundsSerializationModel | None = pydantic.Field(
960 default=None,
961 description="Pixel validity region, if different from the image bounding box.",
962 exclude_if=is_none,
963 )
964 backgrounds: BackgroundMapSerializationModel = pydantic.Field(
965 default_factory=BackgroundMapSerializationModel,
966 description="Background models associated with this image.",
967 )
968 band: str = pydantic.Field(description="Short name of the bandpass filter.")
970 def deserialize(
971 self, archive: InputArchive[Any], *, bbox: Box | None = None, **kwargs: Any
972 ) -> VisitImage:
973 if kwargs: 973 ↛ 974line 973 didn't jump to line 974 because the condition on line 973 was never true
974 raise InvalidParameterError(f"Unrecognized parameters for VisitImage: {set(kwargs.keys())}.")
975 masked_image = super().deserialize(archive, bbox=bbox)
976 try:
977 psf = self.psf.deserialize(archive)
978 except ArchiveReadError as err:
979 # Defer this until/unless somebody actually asks for the PSF.
980 psf = err
981 detector = self.detector.deserialize(archive)
982 aperture_corrections = self.aperture_corrections.deserialize(archive)
983 photometric_scaling = (
984 self.photometric_scaling.deserialize(archive) if self.photometric_scaling is not None else None
985 )
986 return VisitImage(
987 masked_image.image,
988 mask=masked_image.mask,
989 variance=masked_image.variance,
990 psf=psf,
991 sky_projection=masked_image.sky_projection,
992 obs_info=self.obs_info,
993 summary_stats=self.summary_stats,
994 detector=detector,
995 aperture_corrections=aperture_corrections,
996 photometric_scaling=photometric_scaling,
997 bounds=self.bounds.deserialize() if self.bounds is not None else None,
998 backgrounds=self.backgrounds.deserialize(archive),
999 band=self.band,
1000 )._finish_deserialize(self)
1002 def deserialize_component(self, component: str, archive: InputArchive[Any], **kwargs: Any) -> Any:
1003 if kwargs and component not in ("image", "mask", "variance", "masked_image"): 1003 ↛ 1004line 1003 didn't jump to line 1004 because the condition on line 1003 was never true
1004 raise InvalidParameterError(
1005 f"Unsupported parameters for VisitImage component {component}: {set(kwargs.keys())}."
1006 )
1007 if component == "masked_image":
1008 return super().deserialize(archive, **kwargs)
1009 return super().deserialize_component(component, archive, **kwargs)
1012def _obs_info_from_md(
1013 md: MutableMapping[str, Any], visit_info: LegacyVisitInfo | None = None
1014) -> ObservationInfo:
1015 # Try to get an ObservationInfo from the primary header as if
1016 # it's a raw header. Else fallback.
1017 try:
1018 obs_info = ObservationInfo.from_header(md, quiet=True)
1019 except ValueError:
1020 # Not known translator. Must fall back to visit info. If we have
1021 # an actual VisitInfo, serialize it since we know that it will be
1022 # complete.
1023 if visit_info is not None: 1023 ↛ 1035line 1023 didn't jump to line 1035 because the condition on line 1023 was always true
1024 from lsst.afw.image import setVisitInfoMetadata
1025 from lsst.daf.base import PropertyList
1027 pl = PropertyList()
1028 setVisitInfoMetadata(pl, visit_info)
1029 # Merge so that we still have access to butler provenance.
1030 md.update(pl)
1032 # Try the given header looking for VisitInfo hints.
1033 # We get lots of warnings if nothing can be found. Currently
1034 # no way to disable those without capturing them.
1035 obs_info = ObservationInfo.from_header(md, translator_class=VisitInfoTranslator, quiet=True)
1036 return obs_info
1039def _update_obs_info_from_legacy(
1040 obs_info: ObservationInfo,
1041 detector: LegacyDetector | None = None,
1042 filter_label: LegacyFilterLabel | None = None,
1043) -> ObservationInfo:
1044 extra_md: dict[str, str | int] = {}
1046 if filter_label is not None and filter_label.hasBandLabel(): 1046 ↛ 1053line 1046 didn't jump to line 1053 because the condition on line 1046 was always true
1047 extra_md["physical_filter"] = filter_label.physicalLabel
1049 # Fill in detector metadata, check for consistency.
1050 # ObsInfo detector name and group can not be derived from
1051 # the getName() information without knowing how the components
1052 # are separated.
1053 if detector is not None: 1053 ↛ 1060line 1053 didn't jump to line 1060 because the condition on line 1053 was always true
1054 detector_md = {
1055 "detector_num": detector.getId(),
1056 "detector_unique_name": detector.getName(),
1057 }
1058 extra_md.update(detector_md)
1060 obs_info_updates: dict[str, str | int] = {}
1061 for k, v in extra_md.items():
1062 current = getattr(obs_info, k)
1063 if current is None: 1063 ↛ 1066line 1063 didn't jump to line 1066 because the condition on line 1063 was always true
1064 obs_info_updates[k] = v
1065 continue
1066 if current != v:
1067 raise RuntimeError(
1068 f"ObservationInfo contains value for '{k}' that is inconsistent "
1069 f"with given legacy object: {v} != {current}"
1070 )
1072 if obs_info_updates: 1072 ↛ 1074line 1072 didn't jump to line 1074 because the condition on line 1072 was always true
1073 obs_info = obs_info.model_copy(update=obs_info_updates)
1074 return obs_info
1077def _reconcile_detector_serial(obs_info: ObservationInfo, detector: Detector) -> None:
1078 # Some LSSTCam detector serial numbers are/were incorrect in the camera
1079 # geometry (DM-55080), so if they conflict it's the ObservationInfo (from
1080 # the headers) that's correct.
1081 if obs_info.detector_serial is not None and detector.serial != obs_info.detector_serial: 1081 ↛ 1082line 1081 didn't jump to line 1082 because the condition on line 1081 was never true
1082 _LOG.warning(
1083 "Detector serial from ObservationInfo (%s) for detector %d does not agree "
1084 "with camera geometry %s; assuming the former is correct.",
1085 obs_info.detector_serial,
1086 detector.id,
1087 detector.serial,
1088 )
1089 detector._attributes.serial = obs_info.detector_serial
1092def _extract_or_check_value[T](
1093 key: str,
1094 given_value: T | None,
1095 *sources: tuple[str, T | None],
1096) -> T:
1097 # Compare given value against multiple sources. If given value is not
1098 # supplied return the first non-None value in the reference sources.
1099 if given_value is not None: 1099 ↛ 1111line 1099 didn't jump to line 1111 because the condition on line 1099 was always true
1100 for source_name, source_value in sources:
1101 if source_value is not None and source_value != given_value: 1101 ↛ 1102line 1101 didn't jump to line 1102 because the condition on line 1101 was never true
1102 raise ValueError(
1103 f"Given value {given_value!r} does not match {source_value!r} from {source_name}."
1104 )
1105 if source_value is not None:
1106 # Only check the first non-None source rather than checking
1107 # all supplied values.
1108 break
1109 return given_value
1111 for _, source_value in sources:
1112 if source_value is not None:
1113 return source_value
1115 raise ValueError(f"No value found for {key}.")
1118def _extract_or_check_header[T](
1119 key: str, given_value: T | None, header: Any, obs_info_value: T | None, coerce: Callable[[Any], T]
1120) -> T:
1121 hdr_value: T | None = None
1122 if (hdr_raw_value := header.get(key)) is not None: 1122 ↛ 1123line 1122 didn't jump to line 1123 because the condition on line 1122 was never true
1123 hdr_value = coerce(hdr_raw_value)
1124 return _extract_or_check_value(
1125 key, given_value, ("ObservationInfo", obs_info_value), (f"header key {key}", hdr_value)
1126 )
1129def _get_unit_conversion_factor(
1130 original: astropy.units.UnitBase, new: astropy.units.UnitBase
1131) -> float | None:
1132 try:
1133 return original.to(new)
1134 except astropy.units.UnitConversionError:
1135 return None