Coverage for tests/test_imageCutoutsBackendV2.py: 96%

1# This file is part of dax_images_cutout. 

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# This program is free software: you can redistribute it and/or modify 

10# it under the terms of the GNU General Public License as published by 

11# the Free Software Foundation, either version 3 of the License, or 

12# (at your option) any later version. 

13# 

14# This program is distributed in the hope that it will be useful, 

15# but WITHOUT ANY WARRANTY; without even the implied warranty of 

16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

17# GNU General Public License for more details. 

18# 

19# You should have received a copy of the GNU General Public License 

20# along with this program. If not, see <https://www.gnu.org/licenses/>. 

21 

22import os.path 

23import tempfile 

24import unittest 

25 

26import astropy.io.fits 

27import astropy.units as u 

28 

29import lsst.images 

30import lsst.images.serialization 

31import lsst.sphgeom 

32import lsst.utils.tests 

33from lsst.daf.butler import Butler 

34from lsst.dax.images.cutout import CutoutMode, ImageCutoutFactory, projection_finders, stencils 

35 

36TESTDIR = os.path.abspath(os.path.dirname(__file__)) 

37 

38 

39class TestImageCutoutsBackendV2(lsst.utils.tests.TestCase): 

40 """Tests for ImageCutoutsBackend using lsst.images data models.""" 

41 

42 def setUp(self): 

43 collection = "LSSTCam/runs/DRP/DP2/v30_0_8/DM-55060/deep_coadd_rewrite/20260612T200929Z" 

44 self.butler = Butler.from_config(os.path.join(TESTDIR, "butler"), collections=collection) 

45 self.enterContext(self.butler) 

46 

47 # Try: RA = 0:01:01.7 Dec = -3:02:13 

48 point = lsst.sphgeom.LonLat.fromDegrees(0.25708, -3.03694) 

49 radius = lsst.sphgeom.Angle((3 * u.arcsec).to_value(u.rad)) 

50 self.stencil = stencils.SkyCircle(point, radius) 

51 

52 # Projection finders are irrelevant in V2 but the constructor 

53 # still requires at least one. 

54 self.projectionFinders = ( 

55 projection_finders.ReadComponents(), 

56 projection_finders.ReadComponentsAstropyFits(), 

57 ) 

58 

59 self.dataId = {"patch": 76, "tract": 5428, "band": "g", "skymap": "lsst_cells_v2"} 

60 ref = self.butler.find_dataset("deep_coadd", data_id=self.dataId) 

61 assert ref is not None 

62 self.ref = ref 

63 

64 def test_extract_ref(self): 

65 """Test that extract_ref produces a reasonable cutout for all modes.""" 

66 for cutout_mode in CutoutMode: 

67 # Projection finder has no effect for V2 and tempdir is 

68 # not relevant for this test. 

69 proj_finder = self.projectionFinders[0] 

70 cutoutBackend = ImageCutoutFactory(self.butler, proj_finder, ".") 

71 result = cutoutBackend.extract_ref(self.stencil, self.ref, cutout_mode=cutout_mode) 

72 match result.cutout: 

73 case lsst.images.MaskedImage(): 

74 box = result.cutout.bbox 

75 array = result.cutout.image.array 

76 case lsst.images.Image(): 76 ↛ 79line 76 didn't jump to line 79 because the pattern on line 76 always matched

77 array = result.cutout.array 

78 box = result.cutout.bbox 

79 case _: 

80 raise RuntimeError(f"Unexpected cutout type: {type(result.cutout)}") 

81 

82 with self.subTest(cutout_mode=str(cutout_mode)): 

83 self.assertEqual(box.x.size, 29) 

84 self.assertEqual(box.y.size, 28) 

85 

86 # We are reading these values from the fuzzed data that 

87 # has no scientific content but we can test that each 

88 # cutout is returning the same value. These numbers have 

89 # not been validated by external tooling and will change 

90 # on redoing the fuzzing. 

91 self.assertFloatsAlmostEqual(array[14, 15], 2.788926362991333) 

92 self.assertFloatsAlmostEqual(array[1, 1], 7.183607578277588) 

93 self.assertFloatsAlmostEqual(array[27, 28], -11.042882919311523) 

94 

95 def test_outside_stencil_mask(self) -> None: 

96 """The ``OUTSIDE_STENCIL`` plane masks exactly the pixels outside the 

97 circular stencil for every mask-bearing cutout mode. 

98 

99 The image-only modes carry no mask and are exempt. The expected count 

100 is fixed by the stencil geometry and the test data WCS: it is the 

101 circle's complement within the ``29 x 28`` bounding box, so it is a 

102 non-empty, non-full subset (which also guards against an inverted 

103 mask). It will change if the stencil or the fuzzed WCS changes. 

104 """ 

105 image_only_modes = {CutoutMode.IMAGE_ONLY, CutoutMode.ASTROPY_IMAGE} 

106 expected_outside = 195 

107 for cutout_mode in CutoutMode: 

108 with self.subTest(cutout_mode=str(cutout_mode)): 

109 cutoutBackend = ImageCutoutFactory(self.butler, self.projectionFinders[0], ".") 

110 result = cutoutBackend.extract_ref(self.stencil, self.ref, cutout_mode=cutout_mode) 

111 result.mask() 

112 if cutout_mode in image_only_modes: 

113 self.assertIsInstance(result.cutout, lsst.images.Image) 

114 continue 

115 mask = result.cutout.mask 

116 self.assertIn("OUTSIDE_STENCIL", mask.schema.names) 

117 outside = mask.get("OUTSIDE_STENCIL") 

118 total = result.cutout.bbox.x.size * result.cutout.bbox.y.size 

119 n_outside = int(outside.sum()) 

120 # A real circle complement: neither empty nor the full box. 

121 self.assertGreater(n_outside, 0) 

122 self.assertLess(n_outside, total) 

123 self.assertEqual(n_outside, expected_outside) 

124 

125 def test_off_edge_cutout(self) -> None: 

126 """Test that we get a truncated cutout at the edge of the image.""" 

127 # Shift the default position slightly so we fall partly off the edge. 

128 # Shift Y such that the bounding box is [-17:12] vs [0:64] 

129 point = lsst.sphgeom.LonLat.fromDegrees(0.25708, -3.03894) 

130 radius = lsst.sphgeom.Angle((3 * u.arcsec).to_value(u.rad)) 

131 self.stencil = stencils.SkyCircle(point, radius, clip=True) 

132 

133 proj_finder = self.projectionFinders[0] 

134 # Should not write any file out so tempdir is irrelevant. 

135 cutoutBackend = ImageCutoutFactory(self.butler, proj_finder, ".") 

136 result = cutoutBackend.extract_ref(self.stencil, self.ref, cutout_mode=CutoutMode.MASKED_IMAGE) 

137 box = result.cutout.bbox 

138 self.assertEqual(box.x.size, 29) 

139 self.assertEqual(box.y.size, 12) 

140 

141 def test_process_ref(self) -> None: 

142 with tempfile.TemporaryDirectory() as tempdir: 

143 for cutout_mode in CutoutMode: 

144 proj_finder = self.projectionFinders[0] 

145 cutoutBackend = ImageCutoutFactory(self.butler, proj_finder, tempdir) 

146 

147 output = cutoutBackend.process_ref(self.stencil, self.ref, cutout_mode=cutout_mode) 

148 self.assertTrue(output.exists()) 

149 

150 # We should be able to read this back in using generic reader. 

151 result = lsst.images.serialization.read_archive(output) 

152 self.assertIsInstance(result, lsst.images.GeneralizedImage) 

153 

154 def test_process_uuid(self) -> None: 

155 with tempfile.TemporaryDirectory() as tempdir: 

156 for cutout_mode in CutoutMode: 

157 proj_finder = self.projectionFinders[0] 

158 cutoutBackend = ImageCutoutFactory(self.butler, proj_finder, tempdir) 

159 

160 output = cutoutBackend.process_uuid(self.stencil, self.ref.id, cutout_mode=cutout_mode) 

161 self.assertTrue(output.exists()) 

162 

163 # We should be able to read this back in using generic reader. 

164 result = lsst.images.serialization.read_archive(output) 

165 self.assertIsInstance(result, lsst.images.GeneralizedImage) 

166 

167 def test_provenance_in_primary_header(self): 

168 """Cutout provenance must land in the primary FITS header for every 

169 cutout mode, including the native ``lsst.images`` container modes. 

170 """ 

171 provenance_keys = ("CUTVERS",) 

172 # Card comments that fit within the 80-column FITS card limit and so 

173 # survive the round trip, covering both the provenance keys (DATE-CUT) 

174 # and the stencil keys (ST_TYPE). Some cards can potentially have 

175 # truncated comments. 

176 provenance_comments = { 

177 "BTLRUUID": "Butler UUID of full image", 

178 "BTLRNAME": "Butler dataset type", 

179 "ST_TYPE": "Type of stencil used to create this cutout", 

180 "ST_RA": "[deg] Circle center Right Ascension", 

181 "DATE-CUT": "Time of cutout extraction", 

182 } 

183 

184 with tempfile.TemporaryDirectory() as tempdir: 

185 cutoutBackend = ImageCutoutFactory(self.butler, self.projectionFinders[0], tempdir) 

186 for cutout_mode in CutoutMode: 

187 with self.subTest(cutout_mode=str(cutout_mode)): 

188 output = cutoutBackend.process_ref(self.stencil, self.ref, cutout_mode=cutout_mode) 

189 with output.open("rb") as fh, astropy.io.fits.open(fh) as hdul: 

190 header = hdul[0].header 

191 for key in provenance_keys: 

192 self.assertIn( 

193 key, 

194 header, 

195 f"{key} missing from primary header for {cutout_mode}", 

196 ) 

197 self.assertEqual(header["BTLRUUID"], self.ref.id.hex) 

198 self.assertEqual(header["BTLRNAME"], self.ref.datasetType.name) 

199 for key, comment in provenance_comments.items(): 

200 self.assertEqual( 

201 header.comments[key], 

202 comment, 

203 f"{key} comment wrong in primary header for {cutout_mode}", 

204 ) 

205 

206 

207if __name__ == "__main__": 207 ↛ 208line 207 didn't jump to line 208 because the condition on line 207 was never true

208 unittest.main()