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