Coverage for tests/test_cell_coadd.py: 25%

87 statements  

« prev     ^ index     » next       coverage.py v7.14.1, created at 2026-06-10 16:39 +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. 

11 

12from __future__ import annotations 

13 

14import os 

15import pickle 

16import unittest 

17from typing import Any 

18 

19import numpy as np 

20 

21from lsst.images import YX, Box, Interval, get_legacy_deep_coadd_mask_planes 

22from lsst.images.cells import CellCoadd, CellIJ 

23from lsst.images.fits import FitsCompressionOptions 

24from lsst.images.tests import ( 

25 DP2_COADD_DATA_ID, 

26 DP2_COADD_MISSING_CELL, 

27 RoundtripFits, 

28 RoundtripJson, 

29 assert_cell_coadds_equal, 

30 assert_masked_images_equal, 

31 assert_psfs_equal, 

32 compare_cell_coadd_to_legacy, 

33 compare_masked_image_to_legacy, 

34 compare_psf_to_legacy, 

35 compare_sky_projection_to_legacy_wcs, 

36) 

37 

38DATA_DIR = os.environ.get("TESTDATA_IMAGES_DIR", None) 

39 

40 

41@unittest.skipUnless(DATA_DIR is not None, "TESTDATA_IMAGES_DIR is not in the environment.") 

42class CellCoaddTestCase(unittest.TestCase): 

43 """Tests for the CellCoadd class and its many component classes.""" 

44 

45 @classmethod 

46 def setUpClass(cls) -> None: 

47 assert DATA_DIR is not None, "Guaranteed by decorator." 

48 cls.filename = os.path.join(DATA_DIR, "dp2", "legacy", "deep_coadd_cell_predetection.fits") 

49 cls.plane_map = get_legacy_deep_coadd_mask_planes() 

50 cls.missing_cell = CellIJ(**DP2_COADD_MISSING_CELL) 

51 try: 

52 from lsst.cell_coadds import MultipleCellCoadd 

53 

54 cls.legacy_cell_coadd = MultipleCellCoadd.read_fits(cls.filename) 

55 except ImportError: 

56 raise unittest.SkipTest("lsst.cell_coadds could not be imported.") from None 

57 with open(os.path.join(DATA_DIR, "dp2", "legacy", "skyMap.pickle"), "rb") as stream: 

58 cls.skymap = pickle.load(stream) 

59 cls.cell_coadd = CellCoadd.from_legacy( 

60 cls.legacy_cell_coadd, 

61 plane_map=cls.plane_map, 

62 tract_info=cls.skymap[DP2_COADD_DATA_ID["tract"]], 

63 ) 

64 

65 def make_psf_points(self, bbox: Box) -> YX[np.ndarray]: 

66 """Make arrays of points to test PSFs at, given a bbox that is assumed 

67 to be snapped to the cell_coadd grid. 

68 """ 

69 xc, yc = np.meshgrid( 

70 np.arange( 

71 bbox.x.start + self.cell_coadd.grid.cell_shape.x * 0.5, 

72 bbox.x.stop, 

73 self.cell_coadd.grid.cell_shape.x, 

74 ), 

75 np.arange( 

76 bbox.y.start + self.cell_coadd.grid.cell_shape.y * 0.5, 

77 bbox.y.stop, 

78 self.cell_coadd.grid.cell_shape.y, 

79 ), 

80 ) 

81 return YX( 

82 y=yc.ravel() + self.rng.uniform(-0.4, 0.4, size=yc.size), 

83 x=xc.ravel() + self.rng.uniform(-0.4, 0.4, size=xc.size), 

84 ) 

85 

86 def setUp(self) -> None: 

87 self.rng = np.random.default_rng(44) 

88 self.psf_points = self.make_psf_points(self.cell_coadd.bbox) 

89 

90 def test_from_legacy(self) -> None: 

91 """Test constructing a CellCoadd by converting a legacy 

92 lsst.cell_coadds.MultipleCellCoadd. 

93 """ 

94 self.assertEqual(self.cell_coadd.bounds.missing, {self.missing_cell}) 

95 self.assertEqual(self.cell_coadd.bbox, Box.factory[12900:13500, 9600:10050]) 

96 compare_cell_coadd_to_legacy( 

97 self, 

98 self.cell_coadd, 

99 self.legacy_cell_coadd, 

100 tract_bbox=Box.from_legacy(self.skymap[DP2_COADD_DATA_ID["tract"]].getBBox()), 

101 plane_map=self.plane_map, 

102 psf_points=self.psf_points, 

103 ) 

104 

105 def test_roundtrip(self) -> None: 

106 """Test serializing a CellCoadd and reading it back in, including 

107 subimage and component reads. 

108 """ 

109 with RoundtripFits(self, self.cell_coadd, "CellCoadd") as roundtrip: 

110 # Check a subimage read. The subbox only overlaps (but does not 

111 # fully cover) the middle 2 (of 4) cells in y, while covering 

112 # exactly the last column of cells in x. It does not cover the 

113 # missing cell. 

114 subbox = Box.factory[ 

115 self.cell_coadd.bbox.y.start + 252 : self.cell_coadd.bbox.y.stop - 175, 

116 self.cell_coadd.bbox.x.stop - 150 : self.cell_coadd.bbox.x.stop, 

117 ] 

118 subimage = roundtrip.get(bbox=subbox) 

119 assert_masked_images_equal(self, subimage, self.cell_coadd[subbox], expect_view=False) 

120 alternates: dict[str, Any] = {} 

121 with self.subTest(): 

122 subpsf = roundtrip.get("psf", bbox=subbox) 

123 self.assertEqual( 

124 subpsf.bounds.bbox, 

125 Box( 

126 y=Interval.factory[ 

127 self.cell_coadd.bbox.y.start + 150 : self.cell_coadd.bbox.y.stop - 150 

128 ], 

129 x=subbox.x, 

130 ), 

131 ) 

132 assert_psfs_equal(self, subpsf, self.cell_coadd.psf, points=self.make_psf_points(subbox)) 

133 self.assertEqual(roundtrip.get("bbox"), self.cell_coadd.bbox) 

134 alternates = { 

135 k: roundtrip.get(k) 

136 for k in [ 

137 "sky_projection", 

138 "image", 

139 "mask", 

140 "variance", 

141 "psf", 

142 "aperture_corrections", 

143 "provenance", 

144 ] 

145 } 

146 with self.subTest(): 

147 backgrounds = roundtrip.get("backgrounds") 

148 self.assertEqual(backgrounds.keys(), set()) 

149 self.assertIsNone(backgrounds.subtracted) 

150 with roundtrip.inspect() as fits: 

151 for extname in ["IMAGE", "MASK", "VARIANCE", "MASK_FRACTIONS/REJECTED"] + [ 

152 f"NOISE_REALIZATIONS/{n}" for n in range(len(self.cell_coadd.noise_realizations)) 

153 ]: 

154 self.assertEqual(fits[extname].header["ZTILE1"], self.cell_coadd.grid.cell_shape.x) 

155 self.assertEqual(fits[extname].header["ZTILE2"], self.cell_coadd.grid.cell_shape.y) 

156 # Fixture self-consistency: bbox and missing-cell set are what setUp 

157 # claims they are. 

158 self.assertEqual(self.cell_coadd.bounds.missing, {self.missing_cell}) 

159 self.assertEqual(self.cell_coadd.bbox, Box.factory[12900:13500, 9600:10050]) 

160 # Full round-trip fidelity, including background contents. 

161 assert_cell_coadds_equal(self, roundtrip.result, self.cell_coadd, expect_view=False) 

162 compare_cell_coadd_to_legacy( 

163 self, 

164 roundtrip.result, 

165 self.legacy_cell_coadd, 

166 tract_bbox=Box.from_legacy(self.skymap[DP2_COADD_DATA_ID["tract"]].getBBox()), 

167 plane_map=self.plane_map, 

168 alternates=alternates, 

169 psf_points=self.psf_points, 

170 ) 

171 

172 def test_fits_compression(self) -> None: 

173 """Test writing with quantized FITS compression.""" 

174 with RoundtripFits( 

175 self, 

176 self.cell_coadd, 

177 storage_class="CellCoadd", 

178 recipe="lossy16", 

179 compression_options={ 

180 "image": FitsCompressionOptions.LOSSY, 

181 "variance": FitsCompressionOptions.LOSSY, 

182 }, 

183 ) as roundtrip: 

184 with roundtrip.inspect() as fits: 

185 for extname in ["IMAGE", "MASK", "VARIANCE", "MASK_FRACTIONS/REJECTED"] + [ 

186 f"NOISE_REALIZATIONS/{n}" for n in range(len(self.cell_coadd.noise_realizations)) 

187 ]: 

188 with self.subTest(extname=extname): 

189 self.assertEqual(fits[extname].header["ZTILE1"], self.cell_coadd.grid.cell_shape.x) 

190 self.assertEqual(fits[extname].header["ZTILE2"], self.cell_coadd.grid.cell_shape.y) 

191 if extname == "MASK" or extname.startswith("MASK_FRACTIONS"): 

192 self.assertEqual(fits[extname].header["ZCMPTYPE"], "GZIP_2") 

193 else: 

194 self.assertEqual(fits[extname].header["ZCMPTYPE"], "RICE_1") 

195 self.assertEqual(fits[extname].header["ZQUANTIZ"], "SUBTRACTIVE_DITHER_2") 

196 

197 def test_fits_json_consistency(self) -> None: 

198 """FITS and JSON backends produce equal CellCoadds on round-trip.""" 

199 with ( 

200 RoundtripFits(self, self.cell_coadd) as fits_rt, 

201 RoundtripJson(self, self.cell_coadd) as json_rt, 

202 ): 

203 assert_cell_coadds_equal(self, self.cell_coadd, fits_rt.result, expect_view=False) 

204 assert_cell_coadds_equal(self, self.cell_coadd, json_rt.result, expect_view=False) 

205 assert_cell_coadds_equal(self, fits_rt.result, json_rt.result, expect_view=False) 

206 

207 def test_to_legacy(self) -> None: 

208 """Test converting a CellCoadd back into a legacy MultipleCellCoadd.""" 

209 legacy_cell_coadd = self.cell_coadd.to_legacy() 

210 compare_cell_coadd_to_legacy( 

211 self, 

212 self.cell_coadd, 

213 legacy_cell_coadd, 

214 tract_bbox=Box.from_legacy(self.skymap[DP2_COADD_DATA_ID["tract"]].getBBox()), 

215 plane_map=self.plane_map, 

216 psf_points=self.psf_points, 

217 ) 

218 

219 def test_to_legacy_exposure(self) -> None: 

220 """Test converting a CellCoadd back into a legacy Exposure.""" 

221 legacy_exposure = self.cell_coadd.to_legacy_exposure() 

222 

223 self.assertEqual(legacy_exposure.getFilter().bandLabel, self.cell_coadd.band) 

224 self.assertEqual(Box.from_legacy(legacy_exposure.getBBox()), self.cell_coadd.bbox) 

225 compare_masked_image_to_legacy( 

226 self, self.cell_coadd, legacy_exposure.maskedImage, plane_map=self.plane_map, expect_view=True 

227 ) 

228 compare_psf_to_legacy( 

229 self, 

230 self.cell_coadd.psf, 

231 legacy_exposure.getPsf(), 

232 points=self.psf_points, 

233 expect_legacy_raise_on_out_of_bounds=True, 

234 ) 

235 compare_sky_projection_to_legacy_wcs( 

236 self, 

237 self.cell_coadd.sky_projection, 

238 legacy_exposure.getWcs(), 

239 self.cell_coadd.sky_projection.pixel_frame, 

240 subimage_bbox=self.cell_coadd.bbox, 

241 is_fits=True, 

242 ) 

243 

244 

245if __name__ == "__main__": 

246 unittest.main()