Coverage for tests/test_difference_image_extras.py: 38%

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 logging 

16import math 

17import os 

18from typing import Any 

19 

20import astropy.units as u 

21import numpy as np 

22import pytest 

23from astro_metadata_translator import ObservationInfo 

24 

25from lsst.images import ( 

26 Box, 

27 DetectorFrame, 

28 DifferenceImage, 

29 DifferenceImageTemplateInfo, 

30 Image, 

31 MaskPlane, 

32 MaskSchema, 

33) 

34from lsst.images.cameras import Detector 

35from lsst.images.convolution_kernels import ConvolutionKernel, ImageBasisConvolutionKernel 

36from lsst.images.psfs import GaussianPointSpreadFunction 

37from lsst.images.serialization import read_archive 

38from lsst.images.tests import ( 

39 DP2_TEMPLATE_COADD_DATASETS, 

40 DP2_VISIT_DETECTOR_DATA_ID, 

41 RoundtripFits, 

42 assert_values_equal, 

43 make_random_sky_projection, 

44 reset_afw_mask_planes, # noqa: F401 

45) 

46 

47try: 

48 from lsst.afw.image import Exposure as LegacyExposure 

49 from lsst.afw.image import ImageD as LegacyImageD 

50 from lsst.afw.math import Kernel as LegacyKernel 

51 from lsst.afw.table import ExposureCatalog as LegacyExposureCatalog 

52 from lsst.daf.base import PropertyList as LegacyPropertyList 

53 from lsst.geom import Extent2I as LegacyExtent2I 

54 from lsst.meas.algorithms import CoaddPsf as LegacyCoaddPsf 

55except ImportError: 

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

57 type LegacyImageD = Any # type: ignore[no-redef] 

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

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

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

61 type LegacyExposureCatalog = Any # type: ignore[no-redef] 

62 type LegacyExtent2I = Any # type: ignore[no-redef] 

63 

64 

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

66 

67 

68@dataclasses.dataclass 

69class _LegacyTestData: 

70 kernel: LegacyKernel 

71 template_metadata: LegacyPropertyList 

72 template_psf: LegacyCoaddPsf 

73 exposure: LegacyExposure 

74 detector_frame: DetectorFrame 

75 

76 

77@pytest.fixture 

78def legacy_test_data(reset_afw_mask_planes: None) -> _LegacyTestData: # noqa: F811 

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

80 

81 Skips if TESTDATA_IMAGES_DIR is unset or afw is unavailable. 

82 """ 

83 # reset_afw_mask_planes will have already skipped if afw is not available. 

84 from lsst.afw.image import ExposureFitsReader 

85 

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

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

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

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

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

91 kernel = LegacyKernel.readFits(kernel_filename) 

92 template_reader = ExposureFitsReader(template_filename) 

93 template_metadata = template_reader.readMetadata() 

94 template_psf = template_reader.readPsf() 

95 exposure = ExposureFitsReader(exposure_filename).read() 

96 detector_frame = DetectorFrame( 

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

98 ) 

99 return _LegacyTestData( 

100 kernel=kernel, 

101 template_metadata=template_metadata, 

102 template_psf=template_psf, 

103 exposure=exposure, 

104 detector_frame=detector_frame, 

105 ) 

106 

107 

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

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

110 points. 

111 """ 

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

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

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

115 x = round(x) 

116 y = round(y) 

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

118 legacy_im.array[...] = 0.0 

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

120 assert_values_equal(im.array, legacy_im.array, rtol=1e-15, atol=1e-15) 

121 

122 

123def _sanity_check_template_info( 

124 template_info: list[DifferenceImageTemplateInfo], detector_frame: DetectorFrame 

125) -> None: 

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

127 assert len(template_info) == 9 

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

129 assert { 

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

131 for info in template_info 

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

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

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

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

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

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

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

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

140 

141 

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

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

144 attached. 

145 """ 

146 difference_image = DifferenceImage.from_legacy(legacy_test_data.exposure) 

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

148 difference_image.templates = DifferenceImageTemplateInfo.from_legacy( 

149 legacy_test_data.detector_frame, 

150 legacy_test_data.template_psf, 

151 legacy_test_data.template_metadata, 

152 DP2_TEMPLATE_COADD_DATASETS, 

153 ) 

154 return difference_image 

155 

156 

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

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

159 FITS. 

160 """ 

161 difference_image = _make_difference_image(legacy_test_data) 

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

163 pass 

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

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

166 

167 

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

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

170 own. 

171 

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

173 assertions live in `test_roundtrip`. 

174 """ 

175 difference_image = _make_difference_image(legacy_test_data) 

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

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

178 

179 

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

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

182 kernel = ImageBasisConvolutionKernel.from_legacy(legacy_test_data.kernel) 

183 compare_kernel_to_legacy(kernel, legacy_test_data.kernel) 

184 legacy_kernel_2 = kernel.to_legacy() 

185 compare_kernel_to_legacy(kernel, legacy_kernel_2) 

186 

187 

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

189 """Test extracting template information from legacy template_detector 

190 components. 

191 """ 

192 template_info = DifferenceImageTemplateInfo.from_legacy( 

193 legacy_test_data.detector_frame, 

194 legacy_test_data.template_psf, 

195 legacy_test_data.template_metadata, 

196 DP2_TEMPLATE_COADD_DATASETS, 

197 ) 

198 _sanity_check_template_info(template_info, legacy_test_data.detector_frame) 

199 

200 

201def test_template_info_no_overlap_is_skipped( 

202 legacy_test_data: _LegacyTestData, caplog: pytest.LogCaptureFixture 

203) -> None: 

204 """Test that template coadd whose patch does not overlap the detector is 

205 skipped in DifferenceImageTemplateInfo.from_legacy, instead of raising 

206 NoOverlapError. 

207 """ 

208 psf = legacy_test_data.template_psf 

209 n_components = psf.getComponentCount() 

210 assert n_components > 1 

211 

212 # Rebuild the CoaddPsf from its per-component getters, shifting component 

213 # 0's patch bbox far from the detector so the bbox intersection in 

214 # from_legacy is empty (NoOverlapError). Keep the tract/patch columns so 

215 # from_legacy can still resolve its butler info. 

216 schema = LegacyExposureCatalog.Table.makeMinimalSchema() 

217 schema.addField("weight", type="D") 

218 schema.addField("tract", type="I") 

219 schema.addField("patch", type="I") 

220 catalog = LegacyExposureCatalog(schema) 

221 shift = LegacyExtent2I(50000, 50000) 

222 for n in range(n_components): 

223 record = catalog.addNew() 

224 record.setId(psf.getId(n)) 

225 record.setWcs(psf.getWcs(n)) 

226 record.setPsf(psf.getPsf(n)) 

227 record.setValidPolygon(psf.getValidPolygon(n)) 

228 bbox = psf.getBBox(n) 

229 if n == 0: 

230 bbox = bbox.shiftedBy(shift) 

231 record.setBBox(bbox) 

232 record.set("weight", psf.getWeight(n)) 

233 record.set("tract", psf.getTract(n)) 

234 record.set("patch", psf.getPatch(n)) 

235 shifted_psf = LegacyCoaddPsf(catalog, psf.getCoaddWcs(), psf.getAveragePosition()) 

236 

237 skipped_tract = psf.getTract(0) 

238 skipped_patch = psf.getPatch(0) 

239 with caplog.at_level(logging.ERROR): 

240 template_info = DifferenceImageTemplateInfo.from_legacy( 

241 legacy_test_data.detector_frame, 

242 shifted_psf, 

243 legacy_test_data.template_metadata, 

244 DP2_TEMPLATE_COADD_DATASETS, 

245 log=logging.getLogger("test_template_info_no_overlap"), 

246 ) 

247 assert len(template_info) == n_components - 1 

248 assert all((info.tract, info.patch) != (skipped_tract, skipped_patch) for info in template_info) 

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

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

251 assert any( 

252 f"No overlap with tract={skipped_tract}, patch={skipped_patch}" in record.message 

253 for record in caplog.records 

254 ) 

255 

256 

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

258 

259 

260def test_difference_image_repr_str_pinned() -> None: 

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

262 rng = np.random.default_rng(500) 

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

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

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

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

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

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

269 di = DifferenceImage( 

270 image, 

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

272 mask_schema=mask_schema, 

273 sky_projection=sky_projection, 

274 detector=detector, 

275 obs_info=obs_info, 

276 band="r", 

277 ) 

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

279 assert repr(di) == ( 

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

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

282 " dtype=dtype('uint8')))" 

283 )