Coverage for tests/test_cell_coadd.py: 27%

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 RoundtripNdf, 

30 assert_cell_coadds_equal, 

31 assert_images_equal, 

32 assert_masked_images_equal, 

33 assert_psfs_equal, 

34 compare_cell_coadd_to_legacy, 

35 compare_masked_image_to_legacy, 

36 compare_psf_to_legacy, 

37 compare_sky_projection_to_legacy_wcs, 

38) 

39 

40try: 

41 import h5py # noqa: F401 

42 

43 HAVE_H5PY = True 

44except ImportError: 

45 HAVE_H5PY = False 

46 

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

48 

49 

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

51class CellCoaddTestCase(unittest.TestCase): 

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

53 

54 @classmethod 

55 def setUpClass(cls) -> None: 

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

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

58 cls.plane_map = get_legacy_deep_coadd_mask_planes() 

59 cls.missing_cell = CellIJ(**DP2_COADD_MISSING_CELL) 

60 try: 

61 from lsst.cell_coadds import MultipleCellCoadd 

62 

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

64 except ImportError: 

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

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

67 cls.skymap = pickle.load(stream) 

68 cls.cell_coadd = CellCoadd.from_legacy( 

69 cls.legacy_cell_coadd, 

70 plane_map=cls.plane_map, 

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

72 ) 

73 

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

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

76 to be snapped to the cell_coadd grid. 

77 """ 

78 xc, yc = np.meshgrid( 

79 np.arange( 

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

81 bbox.x.stop, 

82 self.cell_coadd.grid.cell_shape.x, 

83 ), 

84 np.arange( 

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

86 bbox.y.stop, 

87 self.cell_coadd.grid.cell_shape.y, 

88 ), 

89 ) 

90 return YX( 

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

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

93 ) 

94 

95 def setUp(self) -> None: 

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

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

98 

99 def test_from_legacy(self) -> None: 

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

101 lsst.cell_coadds.MultipleCellCoadd. 

102 """ 

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

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

105 compare_cell_coadd_to_legacy( 

106 self, 

107 self.cell_coadd, 

108 self.legacy_cell_coadd, 

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

110 plane_map=self.plane_map, 

111 psf_points=self.psf_points, 

112 ) 

113 

114 def test_roundtrip(self) -> None: 

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

116 subimage and component reads. 

117 """ 

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

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

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

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

122 # missing cell. 

123 subbox = Box.factory[ 

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

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

126 ] 

127 subimage = roundtrip.get(bbox=subbox) 

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

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

130 with self.subTest(): 

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

132 self.assertEqual( 

133 subpsf.bounds.bbox, 

134 Box( 

135 y=Interval.factory[ 

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

137 ], 

138 x=subbox.x, 

139 ), 

140 ) 

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

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

143 alternates = { 

144 k: roundtrip.get(k) 

145 for k in [ 

146 "sky_projection", 

147 "image", 

148 "mask", 

149 "variance", 

150 "masked_image", 

151 "psf", 

152 "aperture_corrections", 

153 "provenance", 

154 "backgrounds", 

155 "bbox", 

156 ] 

157 } 

158 # Read all the components at once. 

159 all_components = roundtrip.get("components") 

160 self.assertEqual(set(all_components), set(alternates) - {"masked_image"}) 

161 self.assertEqual(all_components["bbox"], alternates["bbox"]) 

162 assert_psfs_equal(self, all_components["psf"], alternates["psf"]) 

163 assert_images_equal(self, all_components["image"], alternates["image"]) 

164 

165 with self.subTest(): 

166 backgrounds = roundtrip.get("backgrounds") 

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

168 self.assertIsNone(backgrounds.subtracted) 

169 with roundtrip.inspect() as fits: 

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

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

172 ]: 

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

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

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

176 # claims they are. 

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

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

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

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

181 compare_cell_coadd_to_legacy( 

182 self, 

183 roundtrip.result, 

184 self.legacy_cell_coadd, 

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

186 plane_map=self.plane_map, 

187 alternates=alternates, 

188 psf_points=self.psf_points, 

189 ) 

190 

191 def test_fits_compression(self) -> None: 

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

193 with RoundtripFits( 

194 self, 

195 self.cell_coadd, 

196 storage_class="CellCoadd", 

197 recipe="lossy16", 

198 compression_options={ 

199 "image": FitsCompressionOptions.LOSSY, 

200 "variance": FitsCompressionOptions.LOSSY, 

201 }, 

202 ) as roundtrip: 

203 with roundtrip.inspect() as fits: 

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

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

206 ]: 

207 with self.subTest(extname=extname): 

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

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

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

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

212 else: 

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

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

215 

216 def test_fits_json_consistency(self) -> None: 

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

218 with ( 

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

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

221 ): 

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

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

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

225 

226 def test_to_legacy(self) -> None: 

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

228 legacy_cell_coadd = self.cell_coadd.to_legacy() 

229 compare_cell_coadd_to_legacy( 

230 self, 

231 self.cell_coadd, 

232 legacy_cell_coadd, 

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

234 plane_map=self.plane_map, 

235 psf_points=self.psf_points, 

236 ) 

237 

238 def test_to_legacy_exposure(self) -> None: 

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

240 legacy_exposure = self.cell_coadd.to_legacy_exposure() 

241 

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

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

244 compare_masked_image_to_legacy( 

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

246 ) 

247 compare_psf_to_legacy( 

248 self, 

249 self.cell_coadd.psf, 

250 legacy_exposure.getPsf(), 

251 points=self.psf_points, 

252 expect_legacy_raise_on_out_of_bounds=True, 

253 ) 

254 compare_sky_projection_to_legacy_wcs( 

255 self, 

256 self.cell_coadd.sky_projection, 

257 legacy_exposure.getWcs(), 

258 self.cell_coadd.sky_projection.pixel_frame, 

259 subimage_bbox=self.cell_coadd.bbox, 

260 is_fits=True, 

261 ) 

262 

263 @unittest.skipUnless(HAVE_H5PY, "h5py is not installed") 

264 def test_round_trip_ndf(self) -> None: 

265 """NDF round-trip for CellCoadd, exercising hoisted long-named arrays. 

266 

267 This test covers the HDS name-shrinker fix for noise_realizations. 

268 """ 

269 with RoundtripNdf(self, self.cell_coadd, "CellCoadd") as roundtrip: 

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

271 

272 @unittest.skipUnless(HAVE_H5PY, "h5py is not installed") 

273 def test_fits_ndf_consistency(self) -> None: 

274 """FITS and NDF backends produce equal CellCoadds on round-trip.""" 

275 with ( 

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

277 RoundtripNdf(self, self.cell_coadd) as ndf_rt, 

278 ): 

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

280 assert_cell_coadds_equal(self, self.cell_coadd, ndf_rt.result, expect_view=False) 

281 assert_cell_coadds_equal(self, fits_rt.result, ndf_rt.result, expect_view=False) 

282 

283 

284if __name__ == "__main__": 

285 unittest.main()