Coverage for python/lsst/images/_difference_image.py: 24%
207 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-17 09:26 +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__ = ("DifferenceImage", "DifferenceImageSerializationModel", "DifferenceImageTemplateInfo")
16import logging
17import math
18import operator
19import uuid
20from collections.abc import Iterable, Mapping
21from types import EllipsisType
22from typing import TYPE_CHECKING, Any, ClassVar, Literal, cast
24import astropy.units
25import pydantic
26from astro_metadata_translator import ObservationInfo
28from ._backgrounds import BackgroundMap
29from ._geom import Bounds, Box, NoOverlapError
30from ._image import Image
31from ._mask import Mask, MaskPlane, MaskSchema, get_legacy_difference_image_mask_planes
32from ._observation_summary_stats import ObservationSummaryStats
33from ._polygon import Polygon
34from ._transforms import DetectorFrame, SkyProjection, TractFrame, Transform
35from ._visit_image import VisitImage, VisitImageSerializationModel
36from .aperture_corrections import (
37 ApertureCorrectionMap,
38)
39from .cameras import Detector
40from .convolution_kernels import ConvolutionKernel, ConvolutionKernelSerializationModel
41from .describe import DescribeOptions, FieldRole, Report, ReportField, ReportTable
42from .fields import Field
43from .psfs import (
44 PointSpreadFunction,
45)
46from .serialization import (
47 ArchiveReadError,
48 InputArchive,
49 InvalidParameterError,
50 MetadataValue,
51 OutputArchive,
52)
54if TYPE_CHECKING:
55 from lsst.daf.butler import DataId, DimensionRecord
57 try:
58 from lsst.afw.geom import SkyWcs as LegacySkyWcs
59 from lsst.afw.image import Exposure as LegacyExposure
60 from lsst.geom import Box2I as LegacyBox2I
61 from lsst.meas.algorithms import CoaddPsf as LegacyCoaddPsf
62 except ImportError:
63 type LegacyBox2I = Any # type: ignore[no-redef]
64 type LegacyExposure = Any # type: ignore[no-redef]
65 type LegacyCoaddPsf = Any # type: ignore[no-redef]
66 type LegacySkyWcs = Any # type: ignore[no-redef]
69class DifferenceImage(VisitImage):
70 """An image that is the PSF-matched difference of two other images.
72 Parameters
73 ----------
74 image
75 The main image plane. If this has a `SkyProjection`, it will be used
76 for all planes unless a ``sky_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
80 ``sky_projection`` 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 sky_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 sky_projection
91 Projection that maps the pixel grid to the sky. Can only be `None` if
92 a ``sky_projection`` is already attached to ``image``.
93 bounds
94 The region where this image's pixels and other properties are valid.
95 If not provided, the bounding box of the image is used. Other
96 components (``psf``, ``sky_projection``, ``aperture_corrections``,
97 etc.) are assumed to have their own bounds which may or may not be the
98 same as the image bounds. If ``bounds`` extends beyond the image
99 bounding box, the intersection between ``bounds`` and the image
100 bounding box is used instead.
101 obs_info
102 General information about this visit in standardized form.
103 summary_stats
104 Summary statistics associated with this visit. Initialized to default
105 values if not provided.
106 photometric_scaling
107 Field that can be used to multiply a post-ISR image units to yield
108 calibrated image units. This may be a scaling that was already
109 applied (so dividing by it will recover the post-ISR units) or a
110 scaling that has not been applied, depending on ``image.unit``.
111 psf
112 Point-spread function model for this image, or an exception explaining
113 why it could not be read (to be raised if the PSF is requested later).
114 detector
115 Geometry and electronic information for the detector attached to this
116 image.
117 aperture_corrections : `dict` [`str`, `~fields.BaseField`]
118 Mapping from photometry algorithm name to the aperture correction for
119 that algorithm.
120 backgrounds
121 Background models associated with this image.
122 band
123 Name of the passband the image was observed with (this is a shorter,
124 less specific version of ``obs_info.physical_filter``).
125 kernel
126 The convolution kernel used to match the (warped) template to the
127 science image.
128 templates
129 Information about the template coadds that went into this difference
130 image.
131 metadata
132 Arbitrary flexible metadata to associate with the image.
134 Notes
135 -----
136 This class assumes that the difference has been performed on the pixel
137 grid of the 'science image' (i.e. a single observation, like `VisitImage`),
138 and most of the attributes of `DifferenceImage` correspond to the science
139 image. The 'template image' is assumed to be comprised of one or more
140 resampled coadd images stitched together.
142 The `DifferenceImage` class can also be used to represent the stitched
143 template itself; while this makes the naming a bit confusing, the type has
144 the right state to play this role.
145 """
147 def __init__(
148 self,
149 image: Image,
150 *,
151 mask: Mask | None = None,
152 variance: Image | None = None,
153 mask_schema: MaskSchema | None = None,
154 sky_projection: SkyProjection[DetectorFrame] | None = None,
155 bounds: Bounds | None = None,
156 obs_info: ObservationInfo | None = None,
157 summary_stats: ObservationSummaryStats | None = None,
158 photometric_scaling: Field | None = None,
159 psf: PointSpreadFunction | ArchiveReadError,
160 detector: Detector,
161 aperture_corrections: ApertureCorrectionMap | None = None,
162 backgrounds: BackgroundMap | None = None,
163 band: str,
164 kernel: ConvolutionKernel | None = None,
165 templates: Iterable[DifferenceImageTemplateInfo] | None = None,
166 metadata: dict[str, MetadataValue] | None = None,
167 ) -> None:
168 super().__init__(
169 image,
170 mask=mask,
171 variance=variance,
172 mask_schema=mask_schema,
173 sky_projection=sky_projection,
174 bounds=bounds,
175 obs_info=obs_info,
176 summary_stats=summary_stats,
177 photometric_scaling=photometric_scaling,
178 psf=psf,
179 detector=detector,
180 aperture_corrections=aperture_corrections,
181 backgrounds=backgrounds,
182 band=band,
183 metadata=metadata,
184 )
185 self._kernel = kernel
186 self._templates = list(templates) if templates is not None else None
188 @staticmethod
189 def _from_visit_image(
190 visit_image: VisitImage,
191 kernel: ConvolutionKernel | None,
192 templates: Iterable[DifferenceImageTemplateInfo] | None,
193 ) -> DifferenceImage:
194 return visit_image._transfer_metadata(
195 DifferenceImage(
196 visit_image.image,
197 mask=visit_image.mask,
198 variance=visit_image.variance,
199 sky_projection=visit_image.sky_projection,
200 bounds=visit_image.bounds,
201 obs_info=visit_image.obs_info,
202 summary_stats=visit_image.summary_stats,
203 photometric_scaling=visit_image.photometric_scaling,
204 psf=visit_image._psf, # get private attr to avoid triggering on ArchiveReadError early.
205 detector=visit_image.detector,
206 aperture_corrections=visit_image.aperture_corrections,
207 backgrounds=visit_image.backgrounds,
208 kernel=kernel,
209 templates=templates,
210 band=visit_image.band,
211 ),
212 )
214 @property
215 def kernel(self) -> ConvolutionKernel:
216 """The convolution kernel used to match the (warped) template
217 to the science image (`.convolution_kernels.ConvolutionKernel`).
218 """
219 if self._kernel is None:
220 raise AttributeError("This difference image does not have a kernel attached.")
221 return self._kernel
223 @kernel.setter
224 def kernel(self, kernel: ConvolutionKernel) -> None:
225 self._kernel = kernel
227 @kernel.deleter
228 def kernel(self) -> None:
229 self._kernel = None
231 @property
232 def templates(self) -> list[DifferenceImageTemplateInfo]:
233 """Information about the template coadds that went into this
234 difference image (`list` [`DifferenceImageTemplateInfo`]).
235 """
236 if self._templates is None: 236 ↛ 237line 236 didn't jump to line 237 because the condition on line 236 was never true
237 raise AttributeError("This difference image does not have any template information attached.")
238 return self._templates
240 @templates.setter
241 def templates(self, templates: Iterable[DifferenceImageTemplateInfo]) -> None:
242 self._templates = list(templates)
244 @templates.deleter
245 def templates(self) -> None:
246 self._templates = None
248 def __getitem__(self, bbox: Box | EllipsisType) -> DifferenceImage:
249 if bbox is ...:
250 return self
251 return self._from_visit_image(
252 super().__getitem__(bbox), kernel=self._kernel, templates=self._templates
253 )
255 def _describe(self, options: DescribeOptions = DescribeOptions(), /) -> Report:
256 """Return a `Report` describing this difference image.
258 Parameters
259 ----------
260 options : `DescribeOptions`, optional
261 Rendering options; forwarded to the base-class report.
262 """
263 report = super()._describe(options)
264 report.type_name = "DifferenceImage"
265 report.summary = f"DifferenceImage({self.image!s}, {list(self.mask.schema.names)})"
266 if options.brief:
267 return report
268 if self._templates: 268 ↛ 272line 268 didn't jump to line 272 because the condition on line 268 was always true
269 fields, table = DifferenceImageTemplateInfo._describe_templates(self._templates)
270 report.fields.extend(fields)
271 report.tables.append(table)
272 if self._kernel is not None: 272 ↛ 274line 272 didn't jump to line 274 because the condition on line 272 was always true
273 report.children["kernel"] = self._kernel._describe(options.for_child())
274 return report
276 def copy(self, *, copy_detector: bool = False) -> DifferenceImage:
277 """Deep-copy the difference image.
279 Parameters
280 ----------
281 copy_detector
282 Whether to deep-copy the `detector` attribute.
283 """
284 return self._from_visit_image(
285 super().copy(copy_detector=copy_detector), kernel=self._kernel, templates=self._templates
286 )
288 def convert_unit(
289 self,
290 unit: astropy.units.UnitBase = astropy.units.nJy,
291 copy: Literal["as-needed"] | bool = True,
292 copy_detector: bool = False,
293 ) -> DifferenceImage:
294 """Return an equivalent image with different pixel units.
296 Parameters
297 ----------
298 unit
299 The unit to transform to. This may be any of the following:
301 - any unit directly relatable to the current units via Astropy;
302 - any unit relatable to the product of the current units with the
303 `photometric_scaling` (i.e. if the current image is in
304 instrumental units but we know how to calibrate them)
305 - any unit relatable to the quotient of the current units with the
306 `photometric_scaling` (i.e. if the current image is in
307 calibrated units and we want to revert back to instrumental
308 units).
309 copy
310 Whether to copy the images and other components. If `True`, all
311 components that aren't controlled by some other argument will
312 always be deep-copied. If `False`, the operation will fail if the
313 image is not already in the right units. If ``as-needed``, only
314 the image and variance will be copied, and only if they are not
315 already in the right units.
316 copy_detector
317 Whether to deep-copy the `detector` attribute.
319 Returns
320 -------
321 `DifferenceImage`
322 An image with the given units.
323 """
324 return self._from_visit_image(
325 super().convert_unit(unit, copy=copy, copy_detector=copy_detector),
326 kernel=self._kernel,
327 templates=self._templates,
328 )
330 def serialize(self, archive: OutputArchive[Any]) -> DifferenceImageSerializationModel[Any]:
331 result = self._serialize_impl(DifferenceImageSerializationModel, archive)
332 if self._kernel is not None:
333 result.kernel = archive.serialize_direct("kernel", self._kernel.serialize)
334 else:
335 result.kernel = None
336 result.templates = self._templates
337 return result
339 @staticmethod
340 def _get_archive_tree_type[P: pydantic.BaseModel](
341 pointer_type: type[P],
342 ) -> type[DifferenceImageSerializationModel[P]]:
343 """Return the serialization model type for this object for an archive
344 type that uses the given pointer type.
345 """
346 return DifferenceImageSerializationModel[pointer_type] # type: ignore
348 @staticmethod
349 def from_legacy( # type: ignore[override]
350 legacy: LegacyExposure,
351 *,
352 exposure_record: DimensionRecord | None,
353 unit: astropy.units.UnitBase | None = None,
354 plane_map: Mapping[str, MaskPlane] | None = None,
355 ) -> DifferenceImage:
356 """Convert from an `lsst.afw.image.Exposure` instance.
358 Parameters
359 ----------
360 legacy
361 An `lsst.afw.image.Exposure` instance that will share image and
362 variance (but not mask) pixel data with the returned object.
363 exposure_record
364 The ``exposure`` dimension record for this observation, as a
365 source of additional required metadata. It is only consulted when
366 the legacy metadata does not carry an embedded
367 `~astro_metadata_translator.ObservationInfo` (flat ``OBS_INFO``
368 cards); in that case pass `None` explicitly to convert without
369 some ``ObservationInfo`` fields (with a warning).
370 unit
371 Units of the image. If not provided, the ``BUNIT`` metadata
372 key will be used, if available.
373 plane_map
374 A mapping from legacy mask plane name to the new plane name and
375 description. If `None` (default),
376 `get_legacy_difference_image_mask_planes` is used.
377 """
378 if plane_map is None:
379 plane_map = get_legacy_difference_image_mask_planes()
380 return DifferenceImage._from_visit_image(
381 VisitImage.from_legacy(legacy, unit=unit, plane_map=plane_map, exposure_record=exposure_record),
382 kernel=None,
383 templates=None,
384 )
386 def to_legacy(
387 self, *, copy: bool | None = None, plane_map: Mapping[str, MaskPlane] | None = None
388 ) -> LegacyExposure:
389 """Convert to an `lsst.afw.image.Exposure` instance.
391 Parameters
392 ----------
393 copy
394 If `True`, always copy the image and variance pixel data.
395 If `False`, return a view, and raise `TypeError` if the pixel data
396 is read-only (this is not supported by afw). If `None`, only copy
397 if the pixel data is read-only. Mask pixel data is always copied.
398 plane_map
399 A mapping from legacy mask plane name to the new plane name and
400 description. If `None` (default),
401 `get_legacy_difference_image_mask_planes` is used.
402 """
403 if plane_map is None:
404 plane_map = get_legacy_difference_image_mask_planes()
405 return super().to_legacy(copy=copy, plane_map=plane_map)
407 @staticmethod
408 def read_legacy( # type: ignore[override]
409 filename: str,
410 *,
411 exposure_record: DimensionRecord | None,
412 preserve_quantization: bool = False,
413 plane_map: Mapping[str, MaskPlane] | None = None,
414 component: Literal[
415 "bbox",
416 "image",
417 "mask",
418 "variance",
419 "sky_projection",
420 "psf",
421 "detector",
422 "photometric_scaling",
423 "obs_info",
424 "summary_stats",
425 "aperture_corrections",
426 "bounds",
427 "unit",
428 ]
429 | None = None,
430 ) -> Any:
431 """Read a FITS file written by `lsst.afw.image.Exposure.writeFits`.
433 Parameters
434 ----------
435 filename
436 Full name of the file.
437 exposure_record
438 The ``exposure`` dimension record for this observation, as a
439 source of additional required metadata. It is only consulted when
440 the legacy metadata does not carry an embedded
441 `~astro_metadata_translator.ObservationInfo` (flat ``OBS_INFO``
442 cards); in that case pass `None` explicitly to convert without
443 some ``ObservationInfo`` fields (with a warning).
444 preserve_quantization
445 If `True`, ensure that writing the masked image back out again will
446 exactly preserve quantization-compressed pixel values. This causes
447 the image and variance plane arrays to be marked as read-only and
448 stores the original binary table data for those planes in memory.
449 If the `MaskedImage` is copied, the precompressed pixel values are
450 not transferred to the copy.
451 plane_map
452 A mapping from legacy mask plane name to the new plane name and
453 description. If `None` (default),
454 `get_legacy_difference_image_mask_planes` is used.
455 component
456 A component to read instead of the full image.
457 """
458 if plane_map is None:
459 plane_map = get_legacy_difference_image_mask_planes()
460 result = VisitImage.read_legacy(
461 filename,
462 exposure_record=exposure_record,
463 preserve_quantization=preserve_quantization,
464 plane_map=plane_map,
465 component=component,
466 )
467 if component is None:
468 return DifferenceImage._from_visit_image(result, kernel=None, templates=None)
469 return result
472class DifferenceImageTemplateInfo(pydantic.BaseModel, ser_json_inf_nan="constants"):
473 """Information about how a template image contributed to a difference
474 image.
475 """
477 skymap: str = pydantic.Field(description="Name of the skymap that defines the tract/patch tiling.")
478 tract: int = pydantic.Field(description="ID of the tract (each tract is a different projection).")
479 patch: int = pydantic.Field(
480 description="ID of the patch (all patches within a tract share a projection)."
481 )
482 dataset_id: uuid.UUID = pydantic.Field(
483 description="Universally unique butler identifier for this template.",
484 )
485 dataset_run: str = pydantic.Field(description="Name of the butler RUN collection for this template.")
486 bounds: Polygon = pydantic.Field(
487 description=(
488 "The approximate intersection of the template and the science image, "
489 "in the science image's pixel coordinate system."
490 )
491 )
492 psf_shape_xx: float = pydantic.Field(description="Second moment of the effective PSF of the template.")
493 psf_shape_yy: float = pydantic.Field(description="Second moment of the effective PSF of the template.")
494 psf_shape_xy: float = pydantic.Field(description="Second moment of the effective PSF of the template.")
495 psf_shape_flag: bool = pydantic.Field(
496 description="Flag set if the second moments of the effective template PSF could not be computed."
497 )
499 @staticmethod
500 def from_legacy(
501 detector_frame: DetectorFrame,
502 legacy_template_psf: LegacyCoaddPsf,
503 legacy_template_metadata: Mapping[str, Any],
504 coadd_data_ids_by_uuid: Mapping[uuid.UUID, DataId],
505 coadd_dataset_type: str = "template_coadd",
506 log: logging.Logger | None = None,
507 ) -> list[DifferenceImageTemplateInfo]:
508 """Construct a list of template information structs from information
509 stored in a legacy stitched template image.
511 Parameters
512 ----------
513 detector_frame
514 Coordinate system and bounding box of the science image.
515 legacy_template_psf
516 The lazy-evaluation PSF model for the stitched template; used to
517 extract the tract and patch IDs of the coadds actually used and
518 their PSF models.
519 legacy_template_metadata
520 The FITS-style metadata of the stitched template; used to extract
521 butler UUIDs and RUN collection names for all *potential* input
522 coadds.
523 coadd_data_ids_by_uuid
524 A mapping from butler dataset ID to ``{tract, patch, band}`` data
525 ID for all coadds that may have contributed to the template. May
526 be a much larger superset of the needed datasets.
527 coadd_dataset_type
528 The name of the coadd template dataset type.
529 log
530 Logger to use for diagnostic messages.
532 Returns
533 -------
534 `list` [`DifferenceImageTemplateInfo`]
535 One struct for each input coadd that overlaps the science image,
536 sorted by tract and patch.
538 Raises
539 ------
540 KeyError
541 Raised if a PSF component's ``(tract, patch)`` has no matching
542 input coadd in the metadata, or if an input coadd's dataset ID is
543 missing from ``coadd_data_ids_by_uuid``.
544 RuntimeError
545 Raised if the input coadds come from more than one skymap.
547 See Also
548 --------
549 from_legacy_psf
550 """
551 n_inputs = legacy_template_metadata["LSST BUTLER N_INPUTS"]
552 butler_info: dict[tuple[int, int], tuple[uuid.UUID, str]] = {}
553 skymap: str | None = None
554 for n in range(n_inputs):
555 if legacy_template_metadata[f"LSST BUTLER INPUT {n} DATASETTYPE"] == coadd_dataset_type:
556 input_id = uuid.UUID(legacy_template_metadata[f"LSST BUTLER INPUT {n} ID"])
557 input_run = legacy_template_metadata[f"LSST BUTLER INPUT {n} RUN"]
558 input_data_id = coadd_data_ids_by_uuid[input_id]
559 if skymap is None:
560 skymap = cast(str, input_data_id["skymap"])
561 elif skymap != input_data_id["skymap"]:
562 raise RuntimeError("Cannot handle multiple skymaps in the inputs to a single template.")
563 butler_info[cast(int, input_data_id["tract"]), cast(int, input_data_id["patch"])] = (
564 input_id,
565 input_run,
566 )
567 return DifferenceImageTemplateInfo.from_legacy_psf(
568 detector_frame, legacy_template_psf, cast(str, skymap), butler_info, log=log
569 )
571 @staticmethod
572 def from_legacy_psf(
573 detector_frame: DetectorFrame,
574 legacy_template_psf: LegacyCoaddPsf,
575 skymap: str,
576 butler_info: Mapping[tuple[int, int], tuple[uuid.UUID, str]],
577 log: logging.Logger | None = None,
578 ) -> list[DifferenceImageTemplateInfo]:
579 """Construct a list of template information structs from a legacy
580 stitched template PSF and butler provenance held by the caller.
582 Parameters
583 ----------
584 detector_frame
585 Coordinate system and bounding box of the science image.
586 legacy_template_psf
587 The lazy-evaluation PSF model for the stitched template; used to
588 extract the tract and patch IDs of the coadds actually used and
589 their PSF models.
590 skymap
591 Name of the skymap that defines the tract and patch tiling.
592 butler_info
593 A mapping from ``(tract, patch)`` to the butler dataset ID and
594 RUN collection name of that coadd. May be a superset of the
595 coadds that contributed to the template.
596 log
597 Logger to use for diagnostic messages.
599 Returns
600 -------
601 `list` [`DifferenceImageTemplateInfo`]
602 One struct for each input coadd that overlaps the science image,
603 sorted by tract and patch.
605 Raises
606 ------
607 KeyError
608 Raised if a PSF component's ``(tract, patch)`` is missing from
609 ``butler_info``.
611 See Also
612 --------
613 from_legacy
615 Notes
616 -----
617 A task that builds the template has its input coadd references in
618 hand and can call this directly. `from_legacy` is for a caller that
619 has only the stored template, and reads the same provenance out of
620 its FITS metadata.
621 """
622 from lsst.afw.geom import makeWcsPairTransform
624 result: list[DifferenceImageTemplateInfo] = []
625 # A "component" of this PSF is an input {tract, patch} coadd.
626 for n in range(legacy_template_psf.getComponentCount()):
627 tract = legacy_template_psf.getTract(n)
628 patch = legacy_template_psf.getPatch(n)
629 dataset_id, dataset_run = butler_info[tract, patch]
630 patch_bbox = Box.from_legacy(legacy_template_psf.getBBox(n))
631 coadd_frame = TractFrame(
632 skymap=skymap,
633 tract=tract,
634 # This bbox is supposed to be the full tract bbox, but this
635 # frame is just a temporary and we don't have access to that.
636 # (If this ever becomes not-a-temporary, we could add a skymap
637 # argument).
638 bbox=patch_bbox,
639 )
640 detector_to_coadd = Transform.from_legacy(
641 makeWcsPairTransform(
642 # CoaddPsf method names did not anticipate being used for
643 # detector-level templates, so this is confusing:
644 legacy_template_psf.getCoaddWcs(), # this is the template_detector WCS!
645 legacy_template_psf.getWcs(n), # this is the template_coadd WCS!
646 ),
647 detector_frame,
648 coadd_frame,
649 )
650 coadd_to_detector = detector_to_coadd.inverted()
651 try:
652 # We transform the detector bbox to each coadd frame, do the
653 # intersection there, and then transform the intersection back
654 # to the detector frame, because we do not trust detector WCSs
655 # beyond the detector bounding box; they can be polynomials
656 # that extrapolate badly. Coadd WCSs in contrast are simple
657 # projections.
658 tmp_bounds = (
659 Polygon.from_box(detector_frame.bbox)
660 .transform(detector_to_coadd)
661 .intersection(patch_bbox)
662 ).transform(coadd_to_detector)
663 # Unfortunately doing the intersection in the coadd coordinate
664 # system means the transformed intersection might not quite be
665 # contained by the detector bounding box, due to floating-point
666 # round-off error. Intersect one more time to tidy it up.
667 bounds = tmp_bounds.intersection(detector_frame.bbox)
668 assert isinstance(bounds, Polygon), (
669 "The operations above should not change the region's fundamental topology."
670 )
671 except NoOverlapError:
672 if log is not None:
673 log.exception(
674 "No overlap with tract=%s, patch=%s; skipping provenance for that template.",
675 tract,
676 patch,
677 )
678 continue
679 try:
680 psf_shape = legacy_template_psf.computeShape(bounds.centroid.to_legacy_float_point())
681 except Exception:
682 if log is not None:
683 log.exception(
684 "Could not compute PSF shape for template coadd with tract=%s, patch=%s", tract, patch
685 )
686 else:
687 raise
688 psf_shape = None
689 result.append(
690 DifferenceImageTemplateInfo(
691 skymap=skymap,
692 tract=tract,
693 patch=patch,
694 dataset_id=dataset_id,
695 dataset_run=dataset_run,
696 bounds=bounds,
697 psf_shape_xx=psf_shape.getIxx() if psf_shape is not None else math.nan,
698 psf_shape_yy=psf_shape.getIyy() if psf_shape is not None else math.nan,
699 psf_shape_xy=psf_shape.getIxy() if psf_shape is not None else math.nan,
700 psf_shape_flag=psf_shape is None,
701 )
702 )
703 result.sort(key=lambda item: (item.tract, item.patch))
704 return result
706 _HOISTABLE_REPORT_VALUES: ClassVar[tuple[tuple[str, str, str], ...]] = (
707 ("skymap", "skymap", "Skymap"),
708 ("dataset_run", "template run", "Run"),
709 )
710 """Template attributes a report shows once above the templates table when
711 every template carries the same one, as ``(attribute, field label, column
712 header)``.
714 A value belongs here when it is wide enough to crowd the table and usually
715 uniform across a difference image's templates, and when it does not
716 identify the row: see the notes on `_describe_templates`. Adding an
717 attribute here is all it takes to cover it; a value that turns out to vary
718 falls back to its column on its own.
719 """
721 def _report_psf_sigma(self) -> str:
722 """Return the determinant radius of a template's effective PSF,
723 formatted for a report cell.
725 Returns
726 -------
727 sigma : `str`
728 Determinant radius in pixels, or ``"n/a"`` where the second moments
729 do not describe an ellipse.
731 Notes
732 -----
733 This is the same quantity as `ObservationSummaryStats.psfSigma`, so the
734 template PSFs can be compared directly against the science image's.
735 """
736 if self.psf_shape_flag:
737 return "n/a"
738 determinant = self.psf_shape_xx * self.psf_shape_yy - self.psf_shape_xy**2
739 if not math.isfinite(determinant) or determinant <= 0.0:
740 return "n/a"
741 return f"{determinant**0.25:.3f}"
743 @staticmethod
744 def _describe_templates(
745 templates: list[DifferenceImageTemplateInfo],
746 ) -> tuple[list[ReportField], ReportTable]:
747 """Return the report elements describing a difference image's
748 templates.
750 Parameters
751 ----------
752 templates
753 Templates to describe; must not be empty.
755 Returns
756 -------
757 fields : `list` [ `~lsst.images.ReportField` ]
758 The `_HOISTABLE_REPORT_VALUES` that every template shares.
759 table : `~lsst.images.ReportTable`
760 One row per template.
762 Notes
763 -----
764 Some of what a template carries is usually, but not dependably, the
765 same for every template of one difference image, and is long enough to
766 crowd out the columns that do differ. Each such value is shown once
767 above the table where every template shares it, and stays a column
768 where they do not, so nothing is lost when the usual case does not
769 hold. `_HOISTABLE_REPORT_VALUES` lists them.
771 What identifies a row is never hoisted, however uniform it happens to
772 be. A detector overlapping a single tract gives every template the
773 same tract, but tract and patch together name the coadd a row
774 describes, and splitting that pair between a field and a column would
775 leave each row unable to say what it is.
776 """
777 fields: list[ReportField] = []
778 columns: list[str] = []
779 getters: list[Any] = []
780 for attribute, label, column in DifferenceImageTemplateInfo._HOISTABLE_REPORT_VALUES:
781 getter = operator.attrgetter(attribute)
782 values = {getter(template) for template in templates}
783 if len(values) == 1:
784 fields.append(ReportField(label=label, value=values.pop(), role=FieldRole.DERIVED))
785 else:
786 columns.append(column)
787 getters.append(getter)
788 columns.extend(["Tract", "Patch", "PSF \N{GREEK SMALL LETTER SIGMA}"])
789 getters.extend([lambda t: t.tract, lambda t: t.patch, DifferenceImageTemplateInfo._report_psf_sigma])
790 if any(template.psf_shape_flag for template in templates):
791 columns.append("PSF flag")
792 getters.append(lambda t: "set" if t.psf_shape_flag else "")
793 columns.append("Dataset ID")
794 getters.append(lambda t: t.dataset_id)
795 return fields, ReportTable(
796 title="Templates",
797 columns=columns,
798 rows=[[getter(template) for getter in getters] for template in templates],
799 )
802class DifferenceImageSerializationModel[P: pydantic.BaseModel](VisitImageSerializationModel[P]):
803 """A Pydantic model used to represent a serialized `DifferenceImage`."""
805 SCHEMA_NAME: ClassVar[str] = "difference_image"
806 SCHEMA_VERSION: ClassVar[str] = "1.1.0"
807 MIN_READ_VERSION: ClassVar[int] = 1
808 PUBLIC_TYPE: ClassVar[type] = DifferenceImage
810 kernel: ConvolutionKernelSerializationModel | None = pydantic.Field(
811 description="The convolution kernel used to match the (warped) template to the science image."
812 )
813 templates: list[DifferenceImageTemplateInfo] | None = pydantic.Field(
814 description="Information about the template coadds that went into this difference image"
815 )
817 def deserialize(
818 self, archive: InputArchive[Any], *, bbox: Box | None = None, **kwargs: Any
819 ) -> DifferenceImage:
820 if kwargs: 820 ↛ 821line 820 didn't jump to line 821 because the condition on line 820 was never true
821 raise InvalidParameterError(f"Unrecognized parameters for DifferenceImage: {set(kwargs.keys())}.")
822 kernel = self.kernel.deserialize(archive) if self.kernel is not None else None
823 return DifferenceImage._from_visit_image(
824 super().deserialize(archive, bbox=bbox), kernel=kernel, templates=self.templates
825 )
827 def deserialize_component(self, component: str, archive: InputArchive[Any], **kwargs: Any) -> Any:
828 if kwargs and component not in ("image", "mask", "variance", "masked_image"):
829 raise InvalidParameterError(
830 f"Unsupported parameters for DifferenceImage component {component}: {set(kwargs.keys())}."
831 )
832 return super().deserialize_component(component, archive, **kwargs)