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