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