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