Coverage for tests/test_difference_image_extras.py: 54%
207 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-20 09:29 +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 instrument="LSSTCam",
166 )
167 return difference_image
170def test_roundtrip(legacy_test_data: _LegacyTestData) -> None:
171 """Test round-tripping a DifferenceImage with extra components through
172 FITS.
173 """
174 difference_image = _make_difference_image(legacy_test_data)
175 with RoundtripFits(difference_image, storage_class="DifferenceImage") as roundtrip:
176 pass
177 compare_kernel_to_legacy(roundtrip.result.kernel, legacy_test_data.kernel)
178 _sanity_check_template_info(roundtrip.result.templates, legacy_test_data.detector_frame)
179 assert all(info.instrument == "LSSTCam" for info in roundtrip.result.templates)
182def test_kernel_component_read(legacy_test_data: _LegacyTestData) -> None:
183 """Verify the kernel component of a DifferenceImage can be read on its
184 own.
186 Requires a butler; skips when `lsst.daf.butler` is absent. Butler-free
187 assertions live in `test_roundtrip`.
188 """
189 difference_image = _make_difference_image(legacy_test_data)
190 with RoundtripFits(difference_image, storage_class="DifferenceImage") as roundtrip:
191 compare_kernel_to_legacy(roundtrip.get("kernel"), legacy_test_data.kernel)
194def test_difference_kernel(legacy_test_data: _LegacyTestData) -> None:
195 """Test converting a legacy difference kernel to and from the new type."""
196 kernel = ImageBasisConvolutionKernel.from_legacy(legacy_test_data.kernel)
197 compare_kernel_to_legacy(kernel, legacy_test_data.kernel)
198 legacy_kernel_2 = kernel.to_legacy()
199 compare_kernel_to_legacy(kernel, legacy_kernel_2)
202def test_template_info(legacy_test_data: _LegacyTestData) -> None:
203 """Test extracting template information from legacy template_detector
204 components.
205 """
206 template_info = DifferenceImageTemplateInfo.from_legacy(
207 legacy_test_data.detector_frame,
208 legacy_test_data.template_psf,
209 legacy_test_data.template_metadata,
210 DP2_TEMPLATE_COADD_DATASETS,
211 )
212 _sanity_check_template_info(template_info, legacy_test_data.detector_frame)
213 assert all(info.instrument == "UNKNOWN" for info in template_info)
216def test_template_info_instrument(legacy_test_data: _LegacyTestData) -> None:
217 """Test that the instrument passed to from_legacy is recorded on every
218 template.
219 """
220 template_info = DifferenceImageTemplateInfo.from_legacy(
221 legacy_test_data.detector_frame,
222 legacy_test_data.template_psf,
223 legacy_test_data.template_metadata,
224 DP2_TEMPLATE_COADD_DATASETS,
225 instrument="LSSTCam",
226 )
227 assert template_info
228 assert all(info.instrument == "LSSTCam" for info in template_info)
231def test_template_info_no_overlap_is_skipped(
232 legacy_test_data: _LegacyTestData, caplog: pytest.LogCaptureFixture
233) -> None:
234 """Test that template coadd whose patch does not overlap the detector is
235 skipped in DifferenceImageTemplateInfo.from_legacy, instead of raising
236 NoOverlapError.
237 """
238 psf = legacy_test_data.template_psf
239 n_components = psf.getComponentCount()
240 assert n_components > 1
242 # Rebuild the CoaddPsf from its per-component getters, shifting component
243 # 0's patch bbox far from the detector so the bbox intersection in
244 # from_legacy is empty (NoOverlapError). Keep the tract/patch columns so
245 # from_legacy can still resolve its butler info.
246 schema = LegacyExposureCatalog.Table.makeMinimalSchema()
247 schema.addField("weight", type="D")
248 schema.addField("tract", type="I")
249 schema.addField("patch", type="I")
250 catalog = LegacyExposureCatalog(schema)
251 shift = LegacyExtent2I(50000, 50000)
252 for n in range(n_components):
253 record = catalog.addNew()
254 record.setId(psf.getId(n))
255 record.setWcs(psf.getWcs(n))
256 record.setPsf(psf.getPsf(n))
257 record.setValidPolygon(psf.getValidPolygon(n))
258 bbox = psf.getBBox(n)
259 if n == 0:
260 bbox = bbox.shiftedBy(shift)
261 record.setBBox(bbox)
262 record.set("weight", psf.getWeight(n))
263 record.set("tract", psf.getTract(n))
264 record.set("patch", psf.getPatch(n))
265 shifted_psf = LegacyCoaddPsf(catalog, psf.getCoaddWcs(), psf.getAveragePosition())
267 skipped_tract = psf.getTract(0)
268 skipped_patch = psf.getPatch(0)
269 with caplog.at_level(logging.ERROR):
270 template_info = DifferenceImageTemplateInfo.from_legacy(
271 legacy_test_data.detector_frame,
272 shifted_psf,
273 legacy_test_data.template_metadata,
274 DP2_TEMPLATE_COADD_DATASETS,
275 log=logging.getLogger("test_template_info_no_overlap"),
276 )
277 assert len(template_info) == n_components - 1
278 assert all((info.tract, info.patch) != (skipped_tract, skipped_patch) for info in template_info)
279 assert not any(info.psf_shape_flag for info in template_info)
280 assert all(legacy_test_data.detector_frame.bbox.contains(info.bounds.bbox) for info in template_info)
281 assert any(
282 f"No overlap with tract={skipped_tract}, patch={skipped_patch}" in record.message
283 for record in caplog.records
284 )
287LOCAL_DATA_DIR = os.path.join(os.path.dirname(__file__), "data")
290def test_difference_image_repr_str_pinned() -> None:
291 """Pin the exact str and repr output of a DifferenceImage."""
292 rng = np.random.default_rng(500)
293 det_frame = DetectorFrame(instrument="Inst", visit=1234, detector=1, bbox=Box.factory[1:4096, 1:4096])
294 mask_schema = MaskSchema([MaskPlane("M1", "D1")])
295 obs_info = ObservationInfo(instrument="LSSTCam", detector_num=4, physical_filter="r1")
296 detector = read_archive(os.path.join(LOCAL_DATA_DIR, "detector.json"), Detector)
297 image = Image(42, shape=(1024, 1024), unit=u.nJy)
298 sky_projection = make_random_sky_projection(rng, det_frame, det_frame.bbox)
299 di = DifferenceImage(
300 image,
301 psf=GaussianPointSpreadFunction(2.5, stamp_size=33, bounds=Box.factory[-10:10, -12:13]),
302 mask_schema=mask_schema,
303 sky_projection=sky_projection,
304 detector=detector,
305 obs_info=obs_info,
306 band="r",
307 )
308 assert str(di) == "DifferenceImage(Image([y=0:1024, x=0:1024], int64), ['M1'])"
309 assert repr(di) == (
310 "DifferenceImage(Image(..., bbox=Box(y=Interval(start=0, stop=1024), x=Interval(start=0, stop=1024)),"
311 " dtype=dtype('int64')), mask_schema=MaskSchema([MaskPlane(name='M1', description='D1')],"
312 " dtype=dtype('uint8')))"
313 )
316FIXTURE_DIR = Path(__file__).parent / "data" / "schemas"
319def _make_template(
320 tract: int = 1,
321 patch: int = 2,
322 *,
323 skymap: str = "sky",
324 dataset_run: str = "run",
325 psf_shape_xx: float = 4.0,
326 psf_shape_yy: float = 4.0,
327 psf_shape_xy: float = 0.0,
328 psf_shape_flag: bool = False,
329) -> DifferenceImageTemplateInfo:
330 """Return a template info struct with defaults for everything a test does
331 not care about.
332 """
333 return DifferenceImageTemplateInfo(
334 skymap=skymap,
335 tract=tract,
336 patch=patch,
337 dataset_id=uuid.uuid4(),
338 dataset_run=dataset_run,
339 bounds=Polygon(x_vertices=[-0.5, 3.5, -0.5], y_vertices=[-0.5, -0.5, 3.5]),
340 psf_shape_xx=psf_shape_xx,
341 psf_shape_yy=psf_shape_yy,
342 psf_shape_xy=psf_shape_xy,
343 psf_shape_flag=psf_shape_flag,
344 )
347def test_template_info_instrument_defaults_to_unknown() -> None:
348 """A template without a recorded instrument, such as one read from a file
349 written before the field existed, has the instrument ``"UNKNOWN"``.
350 """
351 assert _make_template().instrument == "UNKNOWN"
354def test_describe_templates_hoists_shared_values() -> None:
355 """A value every template shares becomes a field instead of a column."""
356 templates = [_make_template(1, 2), _make_template(1, 3)]
357 fields, table = DifferenceImageTemplateInfo._describe_templates(templates)
358 assert {field.label: field.value for field in fields} == {"skymap": "sky", "template run": "run"}
359 assert table.title == "Templates"
360 assert table.columns == ["Tract", "Patch", "PSF \N{GREEK SMALL LETTER SIGMA}", "Dataset ID"]
361 assert [row[:3] for row in table.rows] == [[1, 2, "2.000"], [1, 3, "2.000"]]
364def test_describe_templates_keeps_varying_values_as_columns() -> None:
365 """A value that differs between templates stays in the table."""
366 templates = [_make_template(skymap="a", dataset_run="r1"), _make_template(skymap="b", dataset_run="r1")]
367 fields, table = DifferenceImageTemplateInfo._describe_templates(templates)
368 # Only the run is shared, so only the run is hoisted.
369 assert [field.label for field in fields] == ["template run"]
370 assert table.columns[0] == "Skymap"
371 assert [row[0] for row in table.rows] == ["a", "b"]
374def test_describe_templates_marks_unusable_psf_shapes() -> None:
375 """A flagged or degenerate PSF shape reports no radius, and the flag
376 column appears only when some template carries it.
377 """
378 good = _make_template(psf_shape_xx=9.0, psf_shape_yy=9.0)
379 flagged = _make_template(psf_shape_flag=True)
380 degenerate = _make_template(psf_shape_xx=1.0, psf_shape_yy=1.0, psf_shape_xy=1.0)
381 _, table = DifferenceImageTemplateInfo._describe_templates([good, degenerate])
382 assert "PSF flag" not in table.columns
383 assert [row[2] for row in table.rows] == ["3.000", "n/a"]
384 _, table = DifferenceImageTemplateInfo._describe_templates([good, flagged])
385 flag_column = table.columns.index("PSF flag")
386 assert [row[2] for row in table.rows] == ["3.000", "n/a"]
387 assert [row[flag_column] for row in table.rows] == ["", "set"]
390def test_difference_image_describe_reports_templates_and_kernel() -> None:
391 """A deserialized difference image describes both of the parts it adds to
392 a visit image.
393 """
394 difference_image = read_archive(str(current_fixture_path(FIXTURE_DIR, "difference_image", variant="dp2")))
395 report = difference_image.describe()
396 tables = {table.title: table for table in report.tables}
397 assert len(tables["Templates"].rows) == len(difference_image.templates)
398 assert "skymap" in {field.label for field in report.fields}
399 assert report.children["kernel"].type_name == "ImageBasisConvolutionKernel"
400 # Both renderers run over the whole thing without error.
401 assert isinstance(report._repr_html_(), str)
402 report.__rich__()
405def test_difference_image_describe_brief_skips_templates_and_kernel() -> None:
406 """Brief reports feed repr and str, which name neither part."""
407 difference_image = read_archive(str(current_fixture_path(FIXTURE_DIR, "difference_image", variant="dp2")))
408 report = difference_image.describe(brief=True)
409 assert report.tables == []
410 assert "kernel" not in report.children
413def test_convolution_kernel_describes_itself() -> None:
414 """A kernel summarizes its basis rather than listing it."""
415 spatial = [ChebyshevField(Box.factory[0:10, 0:20], np.array([[1.0]])) for _ in range(3)]
416 kernel = ImageBasisConvolutionKernel(np.ones((3, 5, 5)), spatial)
417 report = kernel.describe()
418 values = {field.label: field.value for field in report.fields}
419 assert values["basis images"] == 3
420 assert values["spatial variation"] == "ChebyshevField"
421 assert report.children == {}
422 assert str(kernel) == "ImageBasisConvolutionKernel with 3 basis images over [y=0:10, x=0:20]"