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