Coverage for tests/test_difference_image_extras.py: 48%

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. 

11 

12from __future__ import annotations 

13 

14import dataclasses 

15import math 

16import os 

17from typing import Any 

18 

19import astropy.units as u 

20import numpy as np 

21import pytest 

22from astro_metadata_translator import ObservationInfo 

23 

24from lsst.images import ( 

25 Box, 

26 DetectorFrame, 

27 DifferenceImage, 

28 DifferenceImageTemplateInfo, 

29 Image, 

30 MaskPlane, 

31 MaskSchema, 

32) 

33from lsst.images.cameras import Detector 

34from lsst.images.convolution_kernels import ConvolutionKernel, ImageBasisConvolutionKernel 

35from lsst.images.psfs import GaussianPointSpreadFunction 

36from lsst.images.serialization import read_archive 

37from lsst.images.tests import ( 

38 DP2_TEMPLATE_COADD_DATASETS, 

39 DP2_VISIT_DETECTOR_DATA_ID, 

40 RoundtripFits, 

41 assert_close, 

42 make_random_sky_projection, 

43) 

44 

45try: 

46 from lsst.afw.image import Exposure as LegacyExposure 

47 from lsst.afw.image import ImageD as LegacyImageD 

48 from lsst.afw.math import Kernel as LegacyKernel 

49 from lsst.daf.base import PropertyList as LegacyPropertyList 

50 from lsst.meas.algorithms import CoaddPsf as LegacyCoaddPsf 

51except ImportError: 

52 type LegacyExposure = Any # type: ignore[no-redef] 

53 type LegacyKernel = Any # type: ignore[no-redef] 

54 type LegacyPropertyList = Any # type: ignore[no-redef] 

55 type LegacyCoaddPsf = Any # type: ignore[no-redef] 

56 

57 

58EXTERNAL_DATA_DIR = os.environ.get("TESTDATA_IMAGES_DIR", None) 

59 

60 

61@dataclasses.dataclass 

62class _LegacyTestData: 

63 kernel: LegacyKernel 

64 template_metadata: LegacyPropertyList 

65 template_psf: LegacyCoaddPsf 

66 exposure: LegacyExposure 

67 detector_frame: DetectorFrame 

68 

69 

70@pytest.fixture(scope="session") 

71def legacy_test_data() -> _LegacyTestData: 

72 """Return a struct of legacy test objects loaded from EXTERNAL_DATA_DIR. 

73 

74 Skips if TESTDATA_IMAGES_DIR is unset or afw is unavailable. 

75 """ 

76 if EXTERNAL_DATA_DIR is None: 76 ↛ 78line 76 didn't jump to line 78 because the condition on line 76 was always true

77 pytest.skip("TESTDATA_IMAGES_DIR is not in the environment.") 

78 kernel_filename = os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "difference_kernel.fits") 

79 template_filename = os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "template_detector.fits") 

80 exposure_filename = os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "difference_image.fits") 

81 try: 

82 from lsst.afw.image import ExposureFitsReader 

83 except ImportError: 

84 pytest.skip("afw not available; cannot read legacy difference image or components") 

85 kernel = LegacyKernel.readFits(kernel_filename) 

86 template_reader = ExposureFitsReader(template_filename) 

87 template_metadata = template_reader.readMetadata() 

88 template_psf = template_reader.readPsf() 

89 exposure = ExposureFitsReader(exposure_filename).read() 

90 detector_frame = DetectorFrame( 

91 **DP2_VISIT_DETECTOR_DATA_ID, bbox=Box.from_legacy(exposure.getDetector().getBBox()) 

92 ) 

93 return _LegacyTestData( 

94 kernel=kernel, 

95 template_metadata=template_metadata, 

96 template_psf=template_psf, 

97 exposure=exposure, 

98 detector_frame=detector_frame, 

99 ) 

100 

101 

102def compare_kernel_to_legacy(kernel: ConvolutionKernel, legacy_kernel: LegacyKernel) -> None: 

103 """Assert that a ConvolutionKernel matches a legacy Kernel at sampled 

104 points. 

105 """ 

106 xy_array = kernel.bounds.bbox.meshgrid(3) 

107 legacy_im = LegacyImageD(kernel.kernel_bbox.to_legacy()) 

108 for x, y in zip(xy_array.x.flat, xy_array.y.flat): 

109 x = round(x) 

110 y = round(y) 

111 im = kernel.compute_kernel_image(x=x, y=y) 

112 legacy_im.array[...] = 0.0 

113 legacy_kernel.computeImage(legacy_im, doNormalize=False, x=x, y=y) 

114 assert_close(im.array, legacy_im.array, rtol=1e-15, atol=1e-15) 

115 

116 

117def _sanity_check_template_info( 

118 template_info: list[DifferenceImageTemplateInfo], detector_frame: DetectorFrame 

119) -> None: 

120 """Check that a list of DifferenceImageTemplateInfo looks plausible.""" 

121 assert len(template_info) == 9 

122 assert {info.dataset_id for info in template_info} == set(DP2_TEMPLATE_COADD_DATASETS.keys()) 

123 assert { 

124 frozenset({"skymap": info.skymap, "tract": info.tract, "patch": info.patch, "band": "r"}.items()) 

125 for info in template_info 

126 } == {frozenset(v.items()) for v in DP2_TEMPLATE_COADD_DATASETS.values()} 

127 assert not any(info.psf_shape_flag for info in template_info) 

128 assert not any(math.isnan(info.psf_shape_xx) for info in template_info) 

129 assert not any(math.isnan(info.psf_shape_yy) for info in template_info) 

130 assert not any(math.isnan(info.psf_shape_xy) for info in template_info) 

131 assert all(detector_frame.bbox.contains(info.bounds.bbox) for info in template_info) 

132 # Patches overlap, so total area is a bit more than detector area. 

133 assert sum(info.bounds.area for info in template_info) < 1.5 * detector_frame.bbox.area 

134 

135 

136def _make_difference_image(legacy_test_data: _LegacyTestData) -> DifferenceImage: 

137 """Return a DifferenceImage with kernel and template components 

138 attached. 

139 """ 

140 difference_image = DifferenceImage.from_legacy(legacy_test_data.exposure) 

141 difference_image.kernel = ImageBasisConvolutionKernel.from_legacy(legacy_test_data.kernel) 

142 difference_image.templates = DifferenceImageTemplateInfo.from_legacy( 

143 legacy_test_data.detector_frame, 

144 legacy_test_data.template_psf, 

145 legacy_test_data.template_metadata, 

146 DP2_TEMPLATE_COADD_DATASETS, 

147 ) 

148 return difference_image 

149 

150 

151def test_roundtrip(legacy_test_data: _LegacyTestData) -> None: 

152 """Test round-tripping a DifferenceImage with extra components through 

153 FITS. 

154 """ 

155 difference_image = _make_difference_image(legacy_test_data) 

156 with RoundtripFits(difference_image, storage_class="DifferenceImage") as roundtrip: 

157 pass 

158 compare_kernel_to_legacy(roundtrip.result.kernel, legacy_test_data.kernel) 

159 _sanity_check_template_info(roundtrip.result.templates, legacy_test_data.detector_frame) 

160 

161 

162def test_kernel_component_read(legacy_test_data: _LegacyTestData) -> None: 

163 """Verify the kernel component of a DifferenceImage can be read on its 

164 own. 

165 

166 Requires a butler; skips when `lsst.daf.butler` is absent. Butler-free 

167 assertions live in `test_roundtrip`. 

168 """ 

169 difference_image = _make_difference_image(legacy_test_data) 

170 with RoundtripFits(difference_image, storage_class="DifferenceImage") as roundtrip: 

171 compare_kernel_to_legacy(roundtrip.get("kernel"), legacy_test_data.kernel) 

172 

173 

174def test_difference_kernel(legacy_test_data: _LegacyTestData) -> None: 

175 """Test converting a legacy difference kernel to and from the new type.""" 

176 kernel = ImageBasisConvolutionKernel.from_legacy(legacy_test_data.kernel) 

177 compare_kernel_to_legacy(kernel, legacy_test_data.kernel) 

178 legacy_kernel_2 = kernel.to_legacy() 

179 compare_kernel_to_legacy(kernel, legacy_kernel_2) 

180 

181 

182def test_template_info(legacy_test_data: _LegacyTestData) -> None: 

183 """Test extracting template information from legacy template_detector 

184 components. 

185 """ 

186 template_info = DifferenceImageTemplateInfo.from_legacy( 

187 legacy_test_data.detector_frame, 

188 legacy_test_data.template_psf, 

189 legacy_test_data.template_metadata, 

190 DP2_TEMPLATE_COADD_DATASETS, 

191 ) 

192 _sanity_check_template_info(template_info, legacy_test_data.detector_frame) 

193 

194 

195LOCAL_DATA_DIR = os.path.join(os.path.dirname(__file__), "data") 

196 

197 

198def test_difference_image_repr_str_pinned() -> None: 

199 """Pin the exact str and repr output of a DifferenceImage.""" 

200 rng = np.random.default_rng(500) 

201 det_frame = DetectorFrame(instrument="Inst", visit=1234, detector=1, bbox=Box.factory[1:4096, 1:4096]) 

202 mask_schema = MaskSchema([MaskPlane("M1", "D1")]) 

203 obs_info = ObservationInfo(instrument="LSSTCam", detector_num=4, physical_filter="r1") 

204 detector = read_archive(os.path.join(LOCAL_DATA_DIR, "detector.json"), Detector) 

205 image = Image(42, shape=(1024, 1024), unit=u.nJy) 

206 sky_projection = make_random_sky_projection(rng, det_frame, det_frame.bbox) 

207 di = DifferenceImage( 

208 image, 

209 psf=GaussianPointSpreadFunction(2.5, stamp_size=33, bounds=Box.factory[-10:10, -12:13]), 

210 mask_schema=mask_schema, 

211 sky_projection=sky_projection, 

212 detector=detector, 

213 obs_info=obs_info, 

214 band="r", 

215 ) 

216 assert str(di) == "DifferenceImage(Image([y=0:1024, x=0:1024], int64), ['M1'])" 

217 assert repr(di) == ( 

218 "DifferenceImage(Image(..., bbox=Box(y=Interval(start=0, stop=1024), x=Interval(start=0, stop=1024))," 

219 " dtype=dtype('int64')), mask_schema=MaskSchema([MaskPlane(name='M1', description='D1')]," 

220 " dtype=dtype('uint8')))" 

221 )