Coverage for tests/test_ndf_layout.py: 96%

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 

12"""Layout sanity tests for NdfOutputArchive. 

13 

14Opens files written by NdfOutputArchive with raw h5py and verifies 

15the on-disk layout matches the HDS-on-HDF5 / NDF spec. 

16 

17Notes on mask routing 

18--------------------- 

19NDF serialization stores ``Mask`` arrays as a 3-D ``uint8`` DATA primitive 

20whose HDS axes are ``(x, y, mask-byte)``. The HDF5 dataset shape is reversed 

21from that, following hds-v5 convention. It also writes a 2-D ``QUALITY`` 

22view: single-byte masks are copied directly, while wider masks collapse to 

230/1 values. 

24""" 

25 

26from __future__ import annotations 

27 

28import numpy as np 

29import pytest 

30 

31from lsst.images import Box, Image, MaskedImage, MaskPlane, MaskSchema 

32from lsst.images.tests import RoundtripNdf 

33 

34try: 

35 import h5py 

36 

37 from lsst.images.ndf import _hds 

38 

39 HAVE_H5PY = True 

40except ImportError: 

41 HAVE_H5PY = False 

42 

43 

44skip_no_h5py = pytest.mark.skipif(not HAVE_H5PY, reason="h5py is not installed") 

45 

46 

47def _cls(node: h5py.Group) -> str: 

48 """Return the HDS type (CLASS attribute) of an h5py group as a 

49 Python str. 

50 """ 

51 val = node.attrs.get(_hds.ATTR_CLASS) 

52 if val is None: 52 ↛ 54line 52 didn't jump to line 54 because the condition on line 52 was never true

53 # Legacy fallback used by older HDS variants. 

54 val = node.attrs.get("HDSTYPE") 

55 if isinstance(val, bytes): 55 ↛ 57line 55 didn't jump to line 57 because the condition on line 55 was always true

56 return val.decode("ascii") 

57 return str(val) 

58 

59 

60def _hds_type(dataset: h5py.Dataset) -> str: 

61 """Return the HDS primitive type string inferred from a dataset's 

62 numpy dtype. 

63 """ 

64 dataset_type = dataset.id.get_type() 

65 if dataset_type.get_class() == h5py.h5t.BITFIELD: 

66 return "_LOGICAL" 

67 return _hds.hds_type_for_dtype(dataset.dtype) 

68 

69 

70def _hds_shape(dataset: h5py.Dataset) -> tuple[int, ...]: 

71 """Return the dataset shape in HDS/Fortran axis order.""" 

72 return tuple(reversed(dataset.shape)) 

73 

74 

75@skip_no_h5py 

76def test_image_layout() -> None: 

77 """Verify the on-disk layout produced by ``ndf.write()`` for a plain 

78 Image. 

79 """ 

80 image = Image( 

81 np.arange(20, dtype=np.float32).reshape(4, 5), 

82 bbox=Box.factory[10:14, 20:25], 

83 ) 

84 with RoundtripNdf(image) as roundtrip: 

85 f = roundtrip.inspect() 

86 # Root group carries CLASS="NDF". 

87 assert _cls(f["/"]) == "NDF" 

88 

89 # DATA_ARRAY is an ARRAY structure. 

90 assert "DATA_ARRAY" in f 

91 assert _cls(f["/DATA_ARRAY"]) == "ARRAY" 

92 

93 # DATA is a 2-D _REAL primitive whose shape matches the image. 

94 assert "DATA" in f["/DATA_ARRAY"] 

95 ds = f["/DATA_ARRAY/DATA"] 

96 assert _hds_type(ds) == "_REAL" 

97 assert ds.ndim == 2 

98 assert ds.shape == image.array.shape 

99 

100 # ORIGIN stores bbox lower bounds as int64 in (x_min, y_min) order. 

101 assert "ORIGIN" in f["/DATA_ARRAY"] 

102 origin = f["/DATA_ARRAY/ORIGIN"][()] 

103 assert origin.dtype == np.int64 

104 assert int(origin[0]) == 20 # x_min from Box.factory[10:14, 20:25] 

105 assert int(origin[1]) == 10 # y_min 

106 

107 # /MORE is the standard NDF extension container (EXT) and 

108 # /MORE/LSST carries the type "LSST" matching its name. 

109 assert "MORE" in f 

110 assert _cls(f["/MORE"]) == "EXT" 

111 assert "LSST" in f["/MORE"] 

112 assert _cls(f["/MORE/LSST"]) == "LSST" 

113 

114 # Main JSON serialisation tree is present. 

115 assert "JSON" in f["/MORE/LSST"] 

116 

117 

118@skip_no_h5py 

119def test_masked_image_compatible_mask_layout() -> None: 

120 """Verify the on-disk layout for a MaskedImage whose mask fits in a 

121 single byte. 

122 

123 Even though the mask schema has only 2 planes (which would fit in a single 

124 NDF QUALITY byte), MaskedImage writes the native 3-D uint8 backing array 

125 in ``/MORE/LSST/MASK`` and a direct 2-D copy in ``/QUALITY``. 

126 """ 

127 planes = [MaskPlane("BAD", "Bad pixel"), MaskPlane("SAT", "Saturated")] 

128 schema = MaskSchema(planes) # default dtype=uint8, mask_size=1 

129 image = Image( 

130 np.arange(20, dtype=np.float32).reshape(4, 5), 

131 bbox=Box.factory[10:14, 20:25], 

132 ) 

133 # Pass an explicit float64 Image as variance so we can verify _DOUBLE 

134 # on disk (the default variance is float32, matching the image dtype). 

135 variance = Image(np.ones((4, 5), dtype=np.float64), bbox=image.bbox) 

136 masked = MaskedImage(image, mask_schema=schema, variance=variance) 

137 masked.mask.set("BAD", image.array % 2 == 0) 

138 masked.mask.set("SAT", image.array > 10) 

139 

140 with RoundtripNdf(masked) as roundtrip: 

141 f = roundtrip.inspect() 

142 assert "QUALITY" in f 

143 assert _cls(f["/QUALITY"]) == "QUALITY" 

144 assert _cls(f["/QUALITY/QUALITY"]) == "ARRAY" 

145 quality_ds = f["/QUALITY/QUALITY/DATA"] 

146 assert _hds_type(quality_ds) == "_UBYTE" 

147 assert quality_ds.shape == image.array.shape 

148 assert _hds_shape(quality_ds) == (image.array.shape[1], image.array.shape[0]) 

149 np.testing.assert_array_equal(quality_ds[()], masked.mask.array[:, :, 0]) 

150 quality_origin = f["/QUALITY/QUALITY/ORIGIN"] 

151 assert _hds_type(quality_origin) == "_INTEGER" 

152 assert list(quality_origin[()]) == [20, 10] 

153 bad_pixel = f["/QUALITY/QUALITY/BAD_PIXEL"] 

154 assert _hds_type(bad_pixel) == "_LOGICAL" 

155 assert not bad_pixel[()] 

156 assert f["/QUALITY/BADBITS"][()] == 255 

157 

158 # /MORE/LSST/MASK is a sub-NDF (CLASS="NDF") with a 

159 # canonical DATA_ARRAY structure containing DATA + ORIGIN. 

160 assert "MORE" in f 

161 assert "LSST" in f["/MORE"] 

162 assert "MASK" in f["/MORE/LSST"] 

163 assert _cls(f["/MORE/LSST/MASK"]) == "NDF" 

164 assert _cls(f["/MORE/LSST/MASK/DATA_ARRAY"]) == "ARRAY" 

165 mask_ds = f["/MORE/LSST/MASK/DATA_ARRAY/DATA"] 

166 assert _hds_type(mask_ds) == "_UBYTE" 

167 assert mask_ds.ndim == 3 

168 assert mask_ds.shape == (1, 4, 5) 

169 assert _hds_shape(mask_ds) == (5, 4, 1) 

170 origin = f["/MORE/LSST/MASK/DATA_ARRAY/ORIGIN"] 

171 assert origin.dtype == np.int64 

172 # The mask shares the parent image's bbox; the trailing mask 

173 # byte axis keeps a zero origin. 

174 assert list(origin[()]) == [20, 10, 0] 

175 bad_pixel = f["/MORE/LSST/MASK/DATA_ARRAY/BAD_PIXEL"] 

176 assert _hds_type(bad_pixel) == "_LOGICAL" 

177 assert not bad_pixel[()] 

178 

179 # VARIANCE is an ARRAY structure whose DATA is _DOUBLE (float64). 

180 assert "VARIANCE" in f 

181 assert _cls(f["/VARIANCE"]) == "ARRAY" 

182 assert "DATA" in f["/VARIANCE"] 

183 assert _hds_type(f["/VARIANCE/DATA"]) == "_DOUBLE" 

184 

185 

186@skip_no_h5py 

187def test_masked_image_incompatible_mask_layout() -> None: 

188 """Verify the on-disk layout for a MaskedImage with more than 8 

189 mask planes. 

190 

191 A 12-plane uint8 mask has ``mask_size=2`` (two bytes per pixel), and the 

192 on-disk HDS axes are ``(x, y, mask-byte)``. 

193 """ 

194 planes = [MaskPlane(f"P{i}", f"Plane {i}") for i in range(12)] 

195 schema = MaskSchema(planes) # default uint8; mask_size = ceil(12/8) = 2 

196 image = Image( 

197 np.arange(20, dtype=np.float32).reshape(4, 5), 

198 bbox=Box.factory[10:14, 20:25], 

199 ) 

200 masked = MaskedImage(image, mask_schema=schema) 

201 masked.mask.set("P0", image.array % 2 == 0) 

202 masked.mask.set("P11", image.array > 10) 

203 expected_quality = np.any(masked.mask.array != 0, axis=2).astype(np.uint8) 

204 

205 with RoundtripNdf(masked) as roundtrip: 

206 f = roundtrip.inspect() 

207 assert "QUALITY" in f 

208 assert _cls(f["/QUALITY/QUALITY"]) == "ARRAY" 

209 quality_ds = f["/QUALITY/QUALITY/DATA"] 

210 assert _hds_type(quality_ds) == "_UBYTE" 

211 assert quality_ds.shape == image.array.shape 

212 np.testing.assert_array_equal(quality_ds[()], expected_quality) 

213 assert f["/QUALITY/BADBITS"][()] == 255 

214 

215 # /MORE/LSST/MASK is a sub-NDF. 

216 assert "MORE" in f 

217 assert "LSST" in f["/MORE"] 

218 assert "MASK" in f["/MORE/LSST"] 

219 assert _cls(f["/MORE/LSST/MASK"]) == "NDF" 

220 assert _cls(f["/MORE/LSST/MASK/DATA_ARRAY"]) == "ARRAY" 

221 

222 ds = f["/MORE/LSST/MASK/DATA_ARRAY/DATA"] 

223 assert _hds_type(ds) == "_UBYTE" 

224 assert ds.ndim == 3 

225 rows, cols = image.array.shape 

226 assert ds.shape == (2, rows, cols) 

227 assert _hds_shape(ds) == (cols, rows, 2) 

228 bad_pixel = f["/MORE/LSST/MASK/DATA_ARRAY/BAD_PIXEL"] 

229 assert _hds_type(bad_pixel) == "_LOGICAL" 

230 assert not bad_pixel[()] 

231 

232 

233@skip_no_h5py 

234def test_masked_image_many_plane_mask_layout() -> None: 

235 """Verify the on-disk layout for a MaskedImage with more than 31 planes.""" 

236 planes = [MaskPlane(f"P{i}", f"Plane {i}") for i in range(40)] 

237 schema = MaskSchema(planes) 

238 image = Image( 

239 np.arange(20, dtype=np.float32).reshape(4, 5), 

240 bbox=Box.factory[10:14, 20:25], 

241 ) 

242 masked = MaskedImage(image, mask_schema=schema) 

243 masked.mask.set("P0", image.array % 2 == 0) 

244 masked.mask.set("P17", image.array > 10) 

245 masked.mask.set("P39", image.array == 19) 

246 expected_quality = np.any(masked.mask.array != 0, axis=2).astype(np.uint8) 

247 

248 with RoundtripNdf(masked) as roundtrip: 

249 f = roundtrip.inspect() 

250 assert "QUALITY" in f 

251 assert _cls(f["/QUALITY/QUALITY"]) == "ARRAY" 

252 quality_ds = f["/QUALITY/QUALITY/DATA"] 

253 assert _hds_type(quality_ds) == "_UBYTE" 

254 assert quality_ds.shape == image.array.shape 

255 np.testing.assert_array_equal(quality_ds[()], expected_quality) 

256 assert f["/QUALITY/BADBITS"][()] == 255 

257 ds = f["/MORE/LSST/MASK/DATA_ARRAY/DATA"] 

258 assert _hds_type(ds) == "_UBYTE" 

259 assert ds.ndim == 3 

260 rows, cols = image.array.shape 

261 assert ds.shape == (5, rows, cols) 

262 assert _hds_shape(ds) == (cols, rows, 5)