Coverage for python/lsst/images/_difference_image.py: 6%
161 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-30 08:38 +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 uuid
19from collections.abc import Iterable, Mapping
20from types import EllipsisType
21from typing import TYPE_CHECKING, Any, ClassVar, Literal, cast
23import astropy.units
24import pydantic
25from astro_metadata_translator import ObservationInfo
27from ._backgrounds import BackgroundMap
28from ._geom import Bounds, Box
29from ._image import Image
30from ._mask import Mask, MaskPlane, MaskSchema, get_legacy_difference_image_mask_planes
31from ._observation_summary_stats import ObservationSummaryStats
32from ._polygon import Polygon
33from ._transforms import DetectorFrame, SkyProjection, TractFrame, Transform
34from ._visit_image import VisitImage, VisitImageSerializationModel
35from .aperture_corrections import (
36 ApertureCorrectionMap,
37)
38from .cameras import Detector
39from .convolution_kernels import ConvolutionKernel, ConvolutionKernelSerializationModel
40from .describe import DescribeOptions, Report
41from .fields import Field
42from .psfs import (
43 PointSpreadFunction,
44)
45from .serialization import (
46 ArchiveReadError,
47 InputArchive,
48 InvalidParameterError,
49 MetadataValue,
50 OutputArchive,
51)
53if TYPE_CHECKING:
54 from lsst.daf.butler import DataId
56 try:
57 from lsst.afw.geom import SkyWcs as LegacySkyWcs
58 from lsst.afw.image import Exposure as LegacyExposure
59 from lsst.geom import Box2I as LegacyBox2I
60 from lsst.meas.algorithms import CoaddPsf as LegacyCoaddPsf
61 except ImportError:
62 type LegacyBox2I = Any # type: ignore[no-redef]
63 type LegacyExposure = Any # type: ignore[no-redef]
64 type LegacyCoaddPsf = Any # type: ignore[no-redef]
65 type LegacySkyWcs = Any # type: ignore[no-redef]
68class DifferenceImage(VisitImage):
69 """An image that is the PSF-matched difference of two other images.
71 Parameters
72 ----------
73 image
74 The main image plane. If this has a `SkyProjection`, it will be used
75 for all planes unless a ``sky_projection`` is passed separately.
76 mask
77 A bitmask image that annotates the main image plane. Must have the
78 same bounding box as ``image`` if provided. Any attached
79 ``sky_projection`` is replaced (possibly by `None`).
80 variance
81 The per-pixel uncertainty of the main image as an image of variance
82 values. Must have the same bounding box as ``image`` if provided, and
83 its units must be the square of ``image.unit`` or `None`.
84 Values default to ``1.0``. Any attached sky_projection is replaced
85 (possibly by `None`).
86 mask_schema
87 Schema for the mask plane. Must be provided if and only if ``mask`` is
88 not provided.
89 sky_projection
90 Projection that maps the pixel grid to the sky. Can only be `None` if
91 a ``sky_projection`` is already attached to ``image``.
92 bounds
93 The region where this image's pixels and other properties are valid.
94 If not provided, the bounding box of the image is used. Other
95 components (``psf``, ``sky_projection``, ``aperture_corrections``,
96 etc.) are assumed to have their own bounds which may or may not be the
97 same as the image bounds. If ``bounds`` extends beyond the image
98 bounding box, the intersection between ``bounds`` and the image
99 bounding box is used instead.
100 obs_info
101 General information about this visit in standardized form.
102 summary_stats
103 Summary statistics associated with this visit. Initialized to default
104 values if not provided.
105 photometric_scaling
106 Field that can be used to multiply a post-ISR image units to yield
107 calibrated image units. This may be a scaling that was already
108 applied (so dividing by it will recover the post-ISR units) or a
109 scaling that has not been applied, depending on ``image.unit``.
110 psf
111 Point-spread function model for this image, or an exception explaining
112 why it could not be read (to be raised if the PSF is requested later).
113 detector
114 Geometry and electronic information for the detector attached to this
115 image.
116 aperture_corrections : `dict` [`str`, `~fields.BaseField`]
117 Mapping from photometry algorithm name to the aperture correction for
118 that algorithm.
119 backgrounds
120 Background models associated with this image.
121 band
122 Name of the passband the image was observed with (this is a shorter,
123 less specific version of ``obs_info.physical_filter``).
124 kernel
125 The convolution kernel used to match the (warped) template to the
126 science image.
127 templates
128 Information about the template coadds that went into this difference
129 image.
130 metadata
131 Arbitrary flexible metadata to associate with the image.
133 Notes
134 -----
135 This class assumes that the difference has been performed on the pixel
136 grid of the 'science image' (i.e. a single observation, like `VisitImage`),
137 and most of the attributes of `DifferenceImage` correspond to the science
138 image. The 'template image' is assumed to be comprised of one or more
139 resampled coadd images stitched together.
141 The `DifferenceImage` class can also be used to represent the stitched
142 template itself; while this makes the naming a bit confusing, the type has
143 the right state to play this role.
144 """
146 def __init__(
147 self,
148 image: Image,
149 *,
150 mask: Mask | None = None,
151 variance: Image | None = None,
152 mask_schema: MaskSchema | None = None,
153 sky_projection: SkyProjection[DetectorFrame] | None = None,
154 bounds: Bounds | None = None,
155 obs_info: ObservationInfo | None = None,
156 summary_stats: ObservationSummaryStats | None = None,
157 photometric_scaling: Field | None = None,
158 psf: PointSpreadFunction | ArchiveReadError,
159 detector: Detector,
160 aperture_corrections: ApertureCorrectionMap | None = None,
161 backgrounds: BackgroundMap | None = None,
162 band: str,
163 kernel: ConvolutionKernel | None = None,
164 templates: Iterable[DifferenceImageTemplateInfo] | None = None,
165 metadata: dict[str, MetadataValue] | None = None,
166 ) -> None:
167 super().__init__(
168 image,
169 mask=mask,
170 variance=variance,
171 mask_schema=mask_schema,
172 sky_projection=sky_projection,
173 bounds=bounds,
174 obs_info=obs_info,
175 summary_stats=summary_stats,
176 photometric_scaling=photometric_scaling,
177 psf=psf,
178 detector=detector,
179 aperture_corrections=aperture_corrections,
180 backgrounds=backgrounds,
181 band=band,
182 metadata=metadata,
183 )
184 self._kernel = kernel
185 self._templates = list(templates) if templates is not None else None
187 @staticmethod
188 def _from_visit_image(
189 visit_image: VisitImage,
190 kernel: ConvolutionKernel | None,
191 templates: Iterable[DifferenceImageTemplateInfo] | None,
192 ) -> DifferenceImage:
193 return visit_image._transfer_metadata(
194 DifferenceImage(
195 visit_image.image,
196 mask=visit_image.mask,
197 variance=visit_image.variance,
198 sky_projection=visit_image.sky_projection,
199 bounds=visit_image.bounds,
200 obs_info=visit_image.obs_info,
201 summary_stats=visit_image.summary_stats,
202 photometric_scaling=visit_image.photometric_scaling,
203 psf=visit_image._psf, # get private attr to avoid triggering on ArchiveReadError early.
204 detector=visit_image.detector,
205 aperture_corrections=visit_image.aperture_corrections,
206 backgrounds=visit_image.backgrounds,
207 kernel=kernel,
208 templates=templates,
209 band=visit_image.band,
210 ),
211 )
213 @property
214 def kernel(self) -> ConvolutionKernel:
215 """The convolution kernel used to match the (warped) template
216 to the science image (`.convolution_kernels.ConvolutionKernel`).
217 """
218 if self._kernel is None:
219 raise AttributeError("This difference image does not have a kernel attached.")
220 return self._kernel
222 @kernel.setter
223 def kernel(self, kernel: ConvolutionKernel) -> None:
224 self._kernel = kernel
226 @kernel.deleter
227 def kernel(self) -> None:
228 self._kernel = None
230 @property
231 def templates(self) -> list[DifferenceImageTemplateInfo]:
232 """Information about the template coadds that went into this
233 difference image (`list` [`DifferenceImageTemplateInfo`]).
234 """
235 if self._templates is None:
236 raise AttributeError("This difference image does not have any template information attached.")
237 return self._templates
239 @templates.setter
240 def templates(self, templates: Iterable[DifferenceImageTemplateInfo]) -> None:
241 self._templates = list(templates)
243 @templates.deleter
244 def templates(self) -> None:
245 self._templates = None
247 def __getitem__(self, bbox: Box | EllipsisType) -> DifferenceImage:
248 if bbox is ...:
249 return self
250 return self._from_visit_image(
251 super().__getitem__(bbox), kernel=self._kernel, templates=self._templates
252 )
254 def _describe(self, options: DescribeOptions = DescribeOptions(), /) -> Report:
255 """Return a `Report` describing this difference image.
257 Parameters
258 ----------
259 options : `DescribeOptions`, optional
260 Rendering options; forwarded to the base-class report.
261 """
262 report = super()._describe(options)
263 report.type_name = "DifferenceImage"
264 report.summary = f"DifferenceImage({self.image!s}, {list(self.mask.schema.names)})"
265 # ConvolutionKernel does not implement _describe; omit it from the
266 # report until a describe method is added to that class.
267 return report
269 def copy(self, *, copy_detector: bool = False) -> DifferenceImage:
270 """Deep-copy the difference image.
272 Parameters
273 ----------
274 copy_detector
275 Whether to deep-copy the `detector` attribute.
276 """
277 return self._from_visit_image(
278 super().copy(copy_detector=copy_detector), kernel=self._kernel, templates=self._templates
279 )
281 def convert_unit(
282 self,
283 unit: astropy.units.UnitBase = astropy.units.nJy,
284 copy: Literal["as-needed"] | bool = True,
285 copy_detector: bool = False,
286 ) -> DifferenceImage:
287 """Return an equivalent image with different pixel units.
289 Parameters
290 ----------
291 unit
292 The unit to transform to. This may be any of the following:
294 - any unit directly relatable to the current units via Astropy;
295 - any unit relatable to the product of the current units with the
296 `photometric_scaling` (i.e. if the current image is in
297 instrumental units but we know how to calibrate them)
298 - any unit relatable to the quotient of the current units with the
299 `photometric_scaling` (i.e. if the current image is in
300 calibrated units and we want to revert back to instrumental
301 units).
302 copy
303 Whether to copy the images and other components. If `True`, all
304 components that aren't controlled by some other argument will
305 always be deep-copied. If `False`, the operation will fail if the
306 image is not already in the right units. If ``as-needed``, only
307 the image and variance will be copied, and only if they are not
308 already in the right units.
309 copy_detector
310 Whether to deep-copy the `detector` attribute.
312 Returns
313 -------
314 `DifferenceImage`
315 An image with the given units.
316 """
317 return self._from_visit_image(
318 super().convert_unit(unit, copy=copy, copy_detector=copy_detector),
319 kernel=self._kernel,
320 templates=self._templates,
321 )
323 def serialize(self, archive: OutputArchive[Any]) -> DifferenceImageSerializationModel[Any]:
324 result = self._serialize_impl(DifferenceImageSerializationModel, archive)
325 if self._kernel is not None:
326 result.kernel = archive.serialize_direct("kernel", self._kernel.serialize)
327 else:
328 result.kernel = None
329 result.templates = self._templates
330 return result
332 @staticmethod
333 def _get_archive_tree_type[P: pydantic.BaseModel](
334 pointer_type: type[P],
335 ) -> type[DifferenceImageSerializationModel[P]]:
336 """Return the serialization model type for this object for an archive
337 type that uses the given pointer type.
338 """
339 return DifferenceImageSerializationModel[pointer_type] # type: ignore
341 @staticmethod
342 def from_legacy( # type: ignore[override]
343 legacy: LegacyExposure,
344 *,
345 unit: astropy.units.UnitBase | None = None,
346 plane_map: Mapping[str, MaskPlane] | None = None,
347 instrument: str | None = None,
348 visit: int | None = None,
349 ) -> DifferenceImage:
350 """Convert from an `lsst.afw.image.Exposure` instance.
352 Parameters
353 ----------
354 legacy
355 An `lsst.afw.image.Exposure` instance that will share image and
356 variance (but not mask) pixel data with the returned object.
357 unit
358 Units of the image. If not provided, the ``BUNIT`` metadata
359 key will be used, if available.
360 plane_map
361 A mapping from legacy mask plane name to the new plane name and
362 description. If `None` (default)
363 `get_legacy_visit_image_mask_planes` is used.
364 instrument
365 Name of the instrument. Extracted from the metadata if not
366 provided.
367 visit
368 ID of the visit. Extracted from the metadata if not provided.
369 """
370 if plane_map is None:
371 plane_map = get_legacy_difference_image_mask_planes()
372 return DifferenceImage._from_visit_image(
373 VisitImage.from_legacy(
374 legacy, unit=unit, plane_map=plane_map, instrument=instrument, visit=visit
375 ),
376 kernel=None,
377 templates=None,
378 )
380 def to_legacy(
381 self, *, copy: bool | None = None, plane_map: Mapping[str, MaskPlane] | None = None
382 ) -> LegacyExposure:
383 """Convert to an `lsst.afw.image.Exposure` instance.
385 Parameters
386 ----------
387 copy
388 If `True`, always copy the image and variance pixel data.
389 If `False`, return a view, and raise `TypeError` if the pixel data
390 is read-only (this is not supported by afw). If `None`, only copy
391 if the pixel data is read-only. Mask pixel data is always copied.
392 plane_map
393 A mapping from legacy mask plane name to the new plane name and
394 description. If `None` (default),
395 `get_legacy_visit_image_mask_planes` is used.
396 """
397 if plane_map is None:
398 plane_map = get_legacy_difference_image_mask_planes()
399 return super().to_legacy(copy=copy, plane_map=plane_map)
401 @staticmethod
402 def read_legacy( # type: ignore[override]
403 filename: str,
404 *,
405 preserve_quantization: bool = False,
406 plane_map: Mapping[str, MaskPlane] | None = None,
407 instrument: str | None = None,
408 visit: int | None = None,
409 component: Literal[
410 "bbox",
411 "image",
412 "mask",
413 "variance",
414 "sky_projection",
415 "psf",
416 "detector",
417 "photometric_scaling",
418 "obs_info",
419 "summary_stats",
420 "aperture_corrections",
421 ]
422 | None = None,
423 ) -> Any:
424 """Read a FITS file written by `lsst.afw.image.Exposure.writeFits`.
426 Parameters
427 ----------
428 filename
429 Full name of the file.
430 preserve_quantization
431 If `True`, ensure that writing the masked image back out again will
432 exactly preserve quantization-compressed pixel values. This causes
433 the image and variance plane arrays to be marked as read-only and
434 stores the original binary table data for those planes in memory.
435 If the `MaskedImage` is copied, the precompressed pixel values are
436 not transferred to the copy.
437 plane_map
438 A mapping from legacy mask plane name to the new plane name and
439 description. If `None` (default)
440 `get_legacy_visit_image_mask_planes` is used.
441 instrument
442 Name of the instrument. Read from the primary header if not
443 provided.
444 visit
445 ID of the visit. Read from the primary header if not
446 provided.
447 component
448 A component to read instead of the full image.
449 """
450 if plane_map is None:
451 plane_map = get_legacy_difference_image_mask_planes()
452 result = VisitImage.read_legacy(
453 filename,
454 preserve_quantization=preserve_quantization,
455 plane_map=plane_map,
456 instrument=instrument,
457 visit=visit,
458 component=component,
459 )
460 if component is None:
461 return DifferenceImage._from_visit_image(result, kernel=None, templates=None)
462 return result
465class DifferenceImageTemplateInfo(pydantic.BaseModel, ser_json_inf_nan="constants"):
466 """Information about how a template image contributed to a difference
467 image.
468 """
470 skymap: str = pydantic.Field(description="Name of the skymap that defines the tract/patch tiling.")
471 tract: int = pydantic.Field(description="ID of the tract (each tract is a different projection).")
472 patch: int = pydantic.Field(
473 description="ID of the patch (all patches within a tract share a projection)."
474 )
475 dataset_id: uuid.UUID = pydantic.Field(
476 description="Universally unique butler identifier for this template.",
477 )
478 dataset_run: str = pydantic.Field(description="Name of the butler RUN collection for this template.")
479 bounds: Polygon = pydantic.Field(
480 description=(
481 "The approximate intersection of the template and the science image, "
482 "in the science image's pixel coordinate system."
483 )
484 )
485 psf_shape_xx: float = pydantic.Field(description="Second moment of the effective PSF of the template.")
486 psf_shape_yy: float = pydantic.Field(description="Second moment of the effective PSF of the template.")
487 psf_shape_xy: float = pydantic.Field(description="Second moment of the effective PSF of the template.")
488 psf_shape_flag: bool = pydantic.Field(
489 description="Flag set if the second moments of the effective template PSF could not be computed."
490 )
492 @staticmethod
493 def from_legacy(
494 detector_frame: DetectorFrame,
495 legacy_template_psf: LegacyCoaddPsf,
496 legacy_template_metadata: Mapping[str, Any],
497 coadd_data_ids_by_uuid: Mapping[uuid.UUID, DataId],
498 coadd_dataset_type: str = "template_coadd",
499 log: logging.Logger | None = None,
500 ) -> list[DifferenceImageTemplateInfo]:
501 """Construct a list of template information structs from information
502 stored in a legacy stitched template image.
504 Parameters
505 ----------
506 detector_frame
507 Coordinate system and bounding box of the science image.
508 legacy_template_psf
509 The lazy-evaluation PSF model for the stitched template; used to
510 extract the tract and patch IDs of the coadds actually used and
511 their PSF models.
512 legacy_template_metadata
513 The FITS-style metadata of the stitched template; used to extract
514 butler UUIDs and RUN collection names for all *potential* input
515 coadds.
516 coadd_data_ids_by_uuid
517 A mapping from butler dataset ID to ``{tract, patch, band}`` data
518 ID for all coadds that may have contributed to the template. May
519 be a much larger superset of the needed datasets.
520 coadd_dataset_type
521 The name of the coadd template dataset type.
522 log
523 Logger to use for diagnostic messages.
524 """
525 from lsst.afw.geom import makeWcsPairTransform
527 n_inputs = legacy_template_metadata["LSST BUTLER N_INPUTS"]
528 butler_info: dict[tuple[int, int], tuple[uuid.UUID, str]] = {}
529 skymap: str | None = None
530 for n in range(n_inputs):
531 if legacy_template_metadata[f"LSST BUTLER INPUT {n} DATASETTYPE"] == coadd_dataset_type:
532 input_id = uuid.UUID(legacy_template_metadata[f"LSST BUTLER INPUT {n} ID"])
533 input_run = legacy_template_metadata[f"LSST BUTLER INPUT {n} RUN"]
534 input_data_id = coadd_data_ids_by_uuid[input_id]
535 if skymap is None:
536 skymap = cast(str, input_data_id["skymap"])
537 elif skymap != input_data_id["skymap"]:
538 raise RuntimeError("Cannot handle multiple skymaps in the inputs to a single template.")
539 butler_info[cast(int, input_data_id["tract"]), cast(int, input_data_id["patch"])] = (
540 input_id,
541 input_run,
542 )
543 result: list[DifferenceImageTemplateInfo] = []
544 # A "component" of this PSF is an input {tract, patch} coadd.
545 for n in range(legacy_template_psf.getComponentCount()):
546 tract = legacy_template_psf.getTract(n)
547 patch = legacy_template_psf.getPatch(n)
548 dataset_id, dataset_run = butler_info[tract, patch]
549 patch_bbox = Box.from_legacy(legacy_template_psf.getBBox(n))
550 coadd_frame = TractFrame(
551 skymap=skymap,
552 tract=tract,
553 # This bbox is supposed to be the full tract bbox, but this
554 # frame is just a temporary and we don't have access to that.
555 # (If this ever becomes not-a-temporary, we could add a skymap
556 # argument).
557 bbox=patch_bbox,
558 )
559 detector_to_coadd = Transform.from_legacy(
560 makeWcsPairTransform(
561 # CoaddPsf method names did not anticipate being used for
562 # detector-level templates, so this is confusing:
563 legacy_template_psf.getCoaddWcs(), # this is the template_detector WCS!
564 legacy_template_psf.getWcs(n), # this is the template_coadd WCS!
565 ),
566 detector_frame,
567 coadd_frame,
568 )
569 coadd_to_detector = detector_to_coadd.inverted()
570 # We transform the detector bbox to each coadd frame, do the
571 # intersection there, and then transform the intersection back to
572 # the detector frame, because we do not trust detector WCSs beyond
573 # the detector bounding box; they can be polynomials that
574 # extrapolate badly. Coadd WCSs in contrast are simple projections.
575 tmp_bounds = (
576 Polygon.from_box(detector_frame.bbox).transform(detector_to_coadd).intersection(patch_bbox)
577 ).transform(coadd_to_detector)
578 # Unfortunately doing the intersection in the coadd coordinate
579 # system means the transformed intersection might not quite be
580 # contained by the detector bounding box, due to floating-point
581 # round-off error. Intersect one more time to tidy it up.
582 bounds = tmp_bounds.intersection(detector_frame.bbox)
583 assert isinstance(bounds, Polygon), (
584 "The operations above should not change the region's fundamental topology."
585 )
586 try:
587 psf_shape = legacy_template_psf.computeShape(bounds.centroid.to_legacy_float_point())
588 except Exception:
589 if log is not None:
590 log.exception(
591 "Could not compute PSF shape for template coadd with tract=%s, patch=%s", tract, patch
592 )
593 else:
594 raise
595 psf_shape = None
596 result.append(
597 DifferenceImageTemplateInfo(
598 skymap=skymap,
599 tract=tract,
600 patch=patch,
601 dataset_id=dataset_id,
602 dataset_run=dataset_run,
603 bounds=bounds,
604 psf_shape_xx=psf_shape.getIxx() if psf_shape is not None else math.nan,
605 psf_shape_yy=psf_shape.getIyy() if psf_shape is not None else math.nan,
606 psf_shape_xy=psf_shape.getIxy() if psf_shape is not None else math.nan,
607 psf_shape_flag=psf_shape is None,
608 )
609 )
610 result.sort(key=lambda item: (item.tract, item.patch))
611 return result
614class DifferenceImageSerializationModel[P: pydantic.BaseModel](VisitImageSerializationModel[P]):
615 """A Pydantic model used to represent a serialized `DifferenceImage`."""
617 SCHEMA_NAME: ClassVar[str] = "difference_image"
618 SCHEMA_VERSION: ClassVar[str] = "1.0.0.dev0"
619 MIN_READ_VERSION: ClassVar[int] = 1
620 PUBLIC_TYPE: ClassVar[type] = DifferenceImage
622 kernel: ConvolutionKernelSerializationModel | None = pydantic.Field(
623 description="The convolution kernel used to match the (warped) template to the science image."
624 )
625 templates: list[DifferenceImageTemplateInfo] | None = pydantic.Field(
626 description="Information about the template coadds that went into this difference image"
627 )
629 def deserialize(
630 self, archive: InputArchive[Any], *, bbox: Box | None = None, **kwargs: Any
631 ) -> DifferenceImage:
632 if kwargs: 632 ↛ 633line 632 didn't jump to line 633 because the condition on line 632 was never true
633 raise InvalidParameterError(f"Unrecognized parameters for DifferenceImage: {set(kwargs.keys())}.")
634 kernel = self.kernel.deserialize(archive) if self.kernel is not None else None
635 return DifferenceImage._from_visit_image(
636 super().deserialize(archive, bbox=bbox), kernel=kernel, templates=self.templates
637 )
639 def deserialize_component(self, component: str, archive: InputArchive[Any], **kwargs: Any) -> Any:
640 if kwargs and component not in ("image", "mask", "variance", "masked_image"):
641 raise InvalidParameterError(
642 f"Unsupported parameters for DifferenceImage component {component}: {set(kwargs.keys())}."
643 )
644 return super().deserialize_component(component, archive, **kwargs)