Coverage for tests/test_difference_image_extras.py: 55%
199 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-14 19:06 +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
14import dataclasses
15import logging
16import math
17import os
18import uuid
19from pathlib import Path
20from typing import Any
22import astropy.units as u
23import numpy as np
24import pytest
25from astro_metadata_translator import ObservationInfo
27from lsst.daf.butler import DimensionUniverse
28from lsst.images import (
29 Box,
30 DetectorFrame,
31 DifferenceImage,
32 DifferenceImageTemplateInfo,
33 Image,
34 MaskPlane,
35 MaskSchema,
36 Polygon,
37)
38from lsst.images.cameras import Detector
39from lsst.images.convolution_kernels import ConvolutionKernel, ImageBasisConvolutionKernel
40from lsst.images.fields import ChebyshevField
41from lsst.images.psfs import GaussianPointSpreadFunction
42from lsst.images.serialization import read_archive
43from lsst.images.tests import (
44 DP2_TEMPLATE_COADD_DATASETS,
45 DP2_VISIT_DETECTOR_DATA_ID,
46 RoundtripFits,
47 assert_values_equal,
48 current_fixture_path,
49 get_dp2_exposure_record,
50 make_random_sky_projection,
51 reset_afw_mask_planes, # noqa: F401
52)
54try:
55 from lsst.afw.image import Exposure as LegacyExposure
56 from lsst.afw.image import ImageD as LegacyImageD
57 from lsst.afw.math import Kernel as LegacyKernel
58 from lsst.afw.table import ExposureCatalog as LegacyExposureCatalog
59 from lsst.daf.base import PropertyList as LegacyPropertyList
60 from lsst.geom import Extent2I as LegacyExtent2I
61 from lsst.meas.algorithms import CoaddPsf as LegacyCoaddPsf
62except ImportError:
63 type LegacyExposure = Any # type: ignore[no-redef]
64 type LegacyImageD = Any # type: ignore[no-redef]
65 type LegacyKernel = Any # type: ignore[no-redef]
66 type LegacyPropertyList = Any # type: ignore[no-redef]
67 type LegacyCoaddPsf = Any # type: ignore[no-redef]
68 type LegacyExposureCatalog = Any # type: ignore[no-redef]
69 type LegacyExtent2I = Any # type: ignore[no-redef]
72EXTERNAL_DATA_DIR = os.environ.get("TESTDATA_IMAGES_DIR", None)
75@dataclasses.dataclass
76class _LegacyTestData:
77 kernel: LegacyKernel
78 template_metadata: LegacyPropertyList
79 template_psf: LegacyCoaddPsf
80 exposure: LegacyExposure
81 detector_frame: DetectorFrame
82 exposure_record: Any = None
85@pytest.fixture
86def legacy_test_data(reset_afw_mask_planes: None) -> _LegacyTestData: # noqa: F811
87 """Return a struct of legacy test objects loaded from EXTERNAL_DATA_DIR.
89 Skips if TESTDATA_IMAGES_DIR is unset or afw is unavailable.
90 """
91 # reset_afw_mask_planes will have already skipped if afw is not available.
92 from lsst.afw.image import ExposureFitsReader
94 if EXTERNAL_DATA_DIR is None: 94 ↛ 96line 94 didn't jump to line 96 because the condition on line 94 was always true
95 pytest.skip("TESTDATA_IMAGES_DIR is not in the environment.")
96 kernel_filename = os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "difference_kernel.fits")
97 template_filename = os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "template_detector.fits")
98 exposure_filename = os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "difference_image.fits")
99 kernel = LegacyKernel.readFits(kernel_filename)
100 template_reader = ExposureFitsReader(template_filename)
101 template_metadata = template_reader.readMetadata()
102 template_psf = template_reader.readPsf()
103 exposure = ExposureFitsReader(exposure_filename).read()
104 exposure_record = get_dp2_exposure_record(DimensionUniverse())
105 detector_frame = DetectorFrame(
106 **DP2_VISIT_DETECTOR_DATA_ID, bbox=Box.from_legacy(exposure.getDetector().getBBox())
107 )
108 return _LegacyTestData(
109 kernel=kernel,
110 template_metadata=template_metadata,
111 template_psf=template_psf,
112 exposure=exposure,
113 detector_frame=detector_frame,
114 exposure_record=exposure_record,
115 )
118def compare_kernel_to_legacy(kernel: ConvolutionKernel, legacy_kernel: LegacyKernel) -> None:
119 """Assert that a ConvolutionKernel matches a legacy Kernel at sampled
120 points.
121 """
122 xy_array = kernel.bounds.bbox.meshgrid(3)
123 legacy_im = LegacyImageD(kernel.kernel_bbox.to_legacy())
124 for x, y in zip(xy_array.x.flat, xy_array.y.flat):
125 x = round(x)
126 y = round(y)
127 im = kernel.compute_kernel_image(x=x, y=y)
128 legacy_im.array[...] = 0.0
129 legacy_kernel.computeImage(legacy_im, doNormalize=False, x=x, y=y)
130 assert_values_equal(im.array, legacy_im.array, rtol=1e-15, atol=1e-15)
133def _sanity_check_template_info(
134 template_info: list[DifferenceImageTemplateInfo], detector_frame: DetectorFrame
135) -> None:
136 """Check that a list of DifferenceImageTemplateInfo looks plausible."""
137 assert len(template_info) == 9
138 assert {info.dataset_id for info in template_info} == set(DP2_TEMPLATE_COADD_DATASETS.keys())
139 assert {
140 frozenset({"skymap": info.skymap, "tract": info.tract, "patch": info.patch, "band": "r"}.items())
141 for info in template_info
142 } == {frozenset(v.items()) for v in DP2_TEMPLATE_COADD_DATASETS.values()}
143 assert not any(info.psf_shape_flag for info in template_info)
144 assert not any(math.isnan(info.psf_shape_xx) for info in template_info)
145 assert not any(math.isnan(info.psf_shape_yy) for info in template_info)
146 assert not any(math.isnan(info.psf_shape_xy) for info in template_info)
147 assert all(detector_frame.bbox.contains(info.bounds.bbox) for info in template_info)
148 # Patches overlap, so total area is a bit more than detector area.
149 assert sum(info.bounds.area for info in template_info) < 1.5 * detector_frame.bbox.area
152def _make_difference_image(legacy_test_data: _LegacyTestData) -> DifferenceImage:
153 """Return a DifferenceImage with kernel and template components
154 attached.
155 """
156 difference_image = DifferenceImage.from_legacy(
157 legacy_test_data.exposure, exposure_record=legacy_test_data.exposure_record
158 )
159 difference_image.kernel = ImageBasisConvolutionKernel.from_legacy(legacy_test_data.kernel)
160 difference_image.templates = DifferenceImageTemplateInfo.from_legacy(
161 legacy_test_data.detector_frame,
162 legacy_test_data.template_psf,
163 legacy_test_data.template_metadata,
164 DP2_TEMPLATE_COADD_DATASETS,
165 )
166 return difference_image
169def test_roundtrip(legacy_test_data: _LegacyTestData) -> None:
170 """Test round-tripping a DifferenceImage with extra components through
171 FITS.
172 """
173 difference_image = _make_difference_image(legacy_test_data)
174 with RoundtripFits(difference_image, storage_class="DifferenceImage") as roundtrip:
175 pass
176 compare_kernel_to_legacy(roundtrip.result.kernel, legacy_test_data.kernel)
177 _sanity_check_template_info(roundtrip.result.templates, legacy_test_data.detector_frame)
180def test_kernel_component_read(legacy_test_data: _LegacyTestData) -> None:
181 """Verify the kernel component of a DifferenceImage can be read on its
182 own.
184 Requires a butler; skips when `lsst.daf.butler` is absent. Butler-free
185 assertions live in `test_roundtrip`.
186 """
187 difference_image = _make_difference_image(legacy_test_data)
188 with RoundtripFits(difference_image, storage_class="DifferenceImage") as roundtrip:
189 compare_kernel_to_legacy(roundtrip.get("kernel"), legacy_test_data.kernel)
192def test_difference_kernel(legacy_test_data: _LegacyTestData) -> None:
193 """Test converting a legacy difference kernel to and from the new type."""
194 kernel = ImageBasisConvolutionKernel.from_legacy(legacy_test_data.kernel)
195 compare_kernel_to_legacy(kernel, legacy_test_data.kernel)
196 legacy_kernel_2 = kernel.to_legacy()
197 compare_kernel_to_legacy(kernel, legacy_kernel_2)
200def test_template_info(legacy_test_data: _LegacyTestData) -> None:
201 """Test extracting template information from legacy template_detector
202 components.
203 """
204 template_info = DifferenceImageTemplateInfo.from_legacy(
205 legacy_test_data.detector_frame,
206 legacy_test_data.template_psf,
207 legacy_test_data.template_metadata,
208 DP2_TEMPLATE_COADD_DATASETS,
209 )
210 _sanity_check_template_info(template_info, legacy_test_data.detector_frame)
213def test_template_info_no_overlap_is_skipped(
214 legacy_test_data: _LegacyTestData, caplog: pytest.LogCaptureFixture
215) -> None:
216 """Test that template coadd whose patch does not overlap the detector is
217 skipped in DifferenceImageTemplateInfo.from_legacy, instead of raising
218 NoOverlapError.
219 """
220 psf = legacy_test_data.template_psf
221 n_components = psf.getComponentCount()
222 assert n_components > 1
224 # Rebuild the CoaddPsf from its per-component getters, shifting component
225 # 0's patch bbox far from the detector so the bbox intersection in
226 # from_legacy is empty (NoOverlapError). Keep the tract/patch columns so
227 # from_legacy can still resolve its butler info.
228 schema = LegacyExposureCatalog.Table.makeMinimalSchema()
229 schema.addField("weight", type="D")
230 schema.addField("tract", type="I")
231 schema.addField("patch", type="I")
232 catalog = LegacyExposureCatalog(schema)
233 shift = LegacyExtent2I(50000, 50000)
234 for n in range(n_components):
235 record = catalog.addNew()
236 record.setId(psf.getId(n))
237 record.setWcs(psf.getWcs(n))
238 record.setPsf(psf.getPsf(n))
239 record.setValidPolygon(psf.getValidPolygon(n))
240 bbox = psf.getBBox(n)
241 if n == 0:
242 bbox = bbox.shiftedBy(shift)
243 record.setBBox(bbox)
244 record.set("weight", psf.getWeight(n))
245 record.set("tract", psf.getTract(n))
246 record.set("patch", psf.getPatch(n))
247 shifted_psf = LegacyCoaddPsf(catalog, psf.getCoaddWcs(), psf.getAveragePosition())
249 skipped_tract = psf.getTract(0)
250 skipped_patch = psf.getPatch(0)
251 with caplog.at_level(logging.ERROR):
252 template_info = DifferenceImageTemplateInfo.from_legacy(
253 legacy_test_data.detector_frame,
254 shifted_psf,
255 legacy_test_data.template_metadata,
256 DP2_TEMPLATE_COADD_DATASETS,
257 log=logging.getLogger("test_template_info_no_overlap"),
258 )
259 assert len(template_info) == n_components - 1
260 assert all((info.tract, info.patch) != (skipped_tract, skipped_patch) for info in template_info)
261 assert not any(info.psf_shape_flag for info in template_info)
262 assert all(legacy_test_data.detector_frame.bbox.contains(info.bounds.bbox) for info in template_info)
263 assert any(
264 f"No overlap with tract={skipped_tract}, patch={skipped_patch}" in record.message
265 for record in caplog.records
266 )
269LOCAL_DATA_DIR = os.path.join(os.path.dirname(__file__), "data")
272def test_difference_image_repr_str_pinned() -> None:
273 """Pin the exact str and repr output of a DifferenceImage."""
274 rng = np.random.default_rng(500)
275 det_frame = DetectorFrame(instrument="Inst", visit=1234, detector=1, bbox=Box.factory[1:4096, 1:4096])
276 mask_schema = MaskSchema([MaskPlane("M1", "D1")])
277 obs_info = ObservationInfo(instrument="LSSTCam", detector_num=4, physical_filter="r1")
278 detector = read_archive(os.path.join(LOCAL_DATA_DIR, "detector.json"), Detector)
279 image = Image(42, shape=(1024, 1024), unit=u.nJy)
280 sky_projection = make_random_sky_projection(rng, det_frame, det_frame.bbox)
281 di = DifferenceImage(
282 image,
283 psf=GaussianPointSpreadFunction(2.5, stamp_size=33, bounds=Box.factory[-10:10, -12:13]),
284 mask_schema=mask_schema,
285 sky_projection=sky_projection,
286 detector=detector,
287 obs_info=obs_info,
288 band="r",
289 )
290 assert str(di) == "DifferenceImage(Image([y=0:1024, x=0:1024], int64), ['M1'])"
291 assert repr(di) == (
292 "DifferenceImage(Image(..., bbox=Box(y=Interval(start=0, stop=1024), x=Interval(start=0, stop=1024)),"
293 " dtype=dtype('int64')), mask_schema=MaskSchema([MaskPlane(name='M1', description='D1')],"
294 " dtype=dtype('uint8')))"
295 )
298FIXTURE_DIR = Path(__file__).parent / "data" / "schemas"
301def _make_template(
302 tract: int = 1,
303 patch: int = 2,
304 *,
305 skymap: str = "sky",
306 dataset_run: str = "run",
307 psf_shape_xx: float = 4.0,
308 psf_shape_yy: float = 4.0,
309 psf_shape_xy: float = 0.0,
310 psf_shape_flag: bool = False,
311) -> DifferenceImageTemplateInfo:
312 """Return a template info struct with defaults for everything a test does
313 not care about.
314 """
315 return DifferenceImageTemplateInfo(
316 skymap=skymap,
317 tract=tract,
318 patch=patch,
319 dataset_id=uuid.uuid4(),
320 dataset_run=dataset_run,
321 bounds=Polygon(x_vertices=[-0.5, 3.5, -0.5], y_vertices=[-0.5, -0.5, 3.5]),
322 psf_shape_xx=psf_shape_xx,
323 psf_shape_yy=psf_shape_yy,
324 psf_shape_xy=psf_shape_xy,
325 psf_shape_flag=psf_shape_flag,
326 )
329def test_describe_templates_hoists_shared_values() -> None:
330 """A value every template shares becomes a field instead of a column."""
331 templates = [_make_template(1, 2), _make_template(1, 3)]
332 fields, table = DifferenceImageTemplateInfo._describe_templates(templates)
333 assert {field.label: field.value for field in fields} == {"skymap": "sky", "template run": "run"}
334 assert table.title == "Templates"
335 assert table.columns == ["Tract", "Patch", "PSF \N{GREEK SMALL LETTER SIGMA}", "Dataset ID"]
336 assert [row[:3] for row in table.rows] == [[1, 2, "2.000"], [1, 3, "2.000"]]
339def test_describe_templates_keeps_varying_values_as_columns() -> None:
340 """A value that differs between templates stays in the table."""
341 templates = [_make_template(skymap="a", dataset_run="r1"), _make_template(skymap="b", dataset_run="r1")]
342 fields, table = DifferenceImageTemplateInfo._describe_templates(templates)
343 # Only the run is shared, so only the run is hoisted.
344 assert [field.label for field in fields] == ["template run"]
345 assert table.columns[0] == "Skymap"
346 assert [row[0] for row in table.rows] == ["a", "b"]
349def test_describe_templates_marks_unusable_psf_shapes() -> None:
350 """A flagged or degenerate PSF shape reports no radius, and the flag
351 column appears only when some template carries it.
352 """
353 good = _make_template(psf_shape_xx=9.0, psf_shape_yy=9.0)
354 flagged = _make_template(psf_shape_flag=True)
355 degenerate = _make_template(psf_shape_xx=1.0, psf_shape_yy=1.0, psf_shape_xy=1.0)
356 _, table = DifferenceImageTemplateInfo._describe_templates([good, degenerate])
357 assert "PSF flag" not in table.columns
358 assert [row[2] for row in table.rows] == ["3.000", "n/a"]
359 _, table = DifferenceImageTemplateInfo._describe_templates([good, flagged])
360 flag_column = table.columns.index("PSF flag")
361 assert [row[2] for row in table.rows] == ["3.000", "n/a"]
362 assert [row[flag_column] for row in table.rows] == ["", "set"]
365def test_difference_image_describe_reports_templates_and_kernel() -> None:
366 """A deserialized difference image describes both of the parts it adds to
367 a visit image.
368 """
369 difference_image = read_archive(str(current_fixture_path(FIXTURE_DIR, "difference_image", variant="dp2")))
370 report = difference_image.describe()
371 tables = {table.title: table for table in report.tables}
372 assert len(tables["Templates"].rows) == len(difference_image.templates)
373 assert "skymap" in {field.label for field in report.fields}
374 assert report.children["kernel"].type_name == "ImageBasisConvolutionKernel"
375 # Both renderers run over the whole thing without error.
376 assert isinstance(report._repr_html_(), str)
377 report.__rich__()
380def test_difference_image_describe_brief_skips_templates_and_kernel() -> None:
381 """Brief reports feed repr and str, which name neither part."""
382 difference_image = read_archive(str(current_fixture_path(FIXTURE_DIR, "difference_image", variant="dp2")))
383 report = difference_image.describe(brief=True)
384 assert report.tables == []
385 assert "kernel" not in report.children
388def test_convolution_kernel_describes_itself() -> None:
389 """A kernel summarizes its basis rather than listing it."""
390 spatial = [ChebyshevField(Box.factory[0:10, 0:20], np.array([[1.0]])) for _ in range(3)]
391 kernel = ImageBasisConvolutionKernel(np.ones((3, 5, 5)), spatial)
392 report = kernel.describe()
393 values = {field.label: field.value for field in report.fields}
394 assert values["basis images"] == 3
395 assert values["spatial variation"] == "ChebyshevField"
396 assert report.children == {}
397 assert str(kernel) == "ImageBasisConvolutionKernel with 3 basis images over [y=0:10, x=0:20]"