Coverage for tests/test_ndf_output_archive.py: 99%

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 json 

15import tempfile 

16import unittest 

17from unittest import mock 

18 

19import astropy.io.fits 

20import astropy.table 

21import astropy.units as u 

22import numpy as np 

23import pydantic 

24 

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

26from lsst.images._transforms import FrameLookupError, FrameSet, Transform 

27from lsst.images._transforms._frames import DetectorFrame, Frame 

28from lsst.images.fits import ExtensionKey, FitsOpaqueMetadata 

29from lsst.images.serialization import ArrayReferenceModel, InlineArrayModel, read 

30from lsst.images.serialization import open as open_archive 

31from lsst.images.tests import make_random_sky_projection 

32 

33try: 

34 import h5py 

35 

36 from lsst.images.ndf import ( 

37 NdfInputArchive, 

38 NdfOutputArchive, 

39 _hds, 

40 write, 

41 ) 

42 from lsst.images.ndf._hds import DAT__SZNAM 

43 

44 HAVE_H5PY = True 

45except ImportError: 

46 HAVE_H5PY = False 

47 

48 

49class TinyFrameSet(FrameSet): 

50 """Minimal concrete frame-set for archive bookkeeping tests.""" 

51 

52 def __contains__(self, frame: Frame) -> bool: 

53 return False 

54 

55 def __getitem__[I: Frame, O: Frame](self, key: tuple[I, O]) -> Transform[I, O]: 

56 raise FrameLookupError(key) 

57 

58 

59class TinyTree(pydantic.BaseModel): 

60 """A trivial Pydantic model used as a serialization stand-in.""" 

61 

62 name: str 

63 

64 

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

66class NdfOutputArchiveBasicsTestCase(unittest.TestCase): 

67 """Tests for `NdfOutputArchive` constructor and `serialize_direct`.""" 

68 

69 def test_serialize_direct_calls_serializer_with_nested_archive(self): 

70 with tempfile.NamedTemporaryFile(suffix=".sdf") as tmp: 

71 with h5py.File(tmp.name, "w") as f: 

72 arch = NdfOutputArchive(f) 

73 tree = arch.serialize_direct("top", lambda nested: TinyTree(name="hello")) 

74 self.assertEqual(tree.name, "hello") 

75 

76 def test_constructor_marks_root_as_ndf(self): 

77 """The constructor should set CLASS=NDF on the root group so that 

78 Starlink tools recognise the file as an NDF. 

79 """ 

80 with tempfile.NamedTemporaryFile(suffix=".sdf") as tmp: 

81 with h5py.File(tmp.name, "w") as f: 

82 NdfOutputArchive(f) 

83 with h5py.File(tmp.name, "r") as f: 

84 self.assertEqual(f["/"].attrs["CLASS"], b"NDF") 

85 

86 

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

88class NdfOutputArchiveAddArrayTestCase(unittest.TestCase): 

89 """Tests for `NdfOutputArchive.add_array` routing.""" 

90 

91 def test_top_level_image_routes_to_data_array(self): 

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

93 with tempfile.NamedTemporaryFile(suffix=".sdf") as tmp: 

94 with h5py.File(tmp.name, "w") as f: 

95 arch = NdfOutputArchive(f) 

96 ref = arch.add_array(data, name="image") 

97 self.assertEqual(ref.source, "ndf:/DATA_ARRAY/DATA") 

98 with h5py.File(tmp.name, "r") as f: 

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

100 self.assertEqual(ds.dtype, np.float32) 

101 np.testing.assert_array_equal(ds[()], data) 

102 self.assertEqual(f["/DATA_ARRAY"].attrs["CLASS"], b"ARRAY") 

103 origin = f["/DATA_ARRAY/ORIGIN"] 

104 self.assertEqual(origin.dtype, np.int64) 

105 self.assertEqual(origin.shape, (2,)) 

106 

107 def test_top_level_variance_routes_to_variance(self): 

108 data = np.full((3, 3), 0.5, dtype=np.float64) 

109 with tempfile.NamedTemporaryFile(suffix=".sdf") as tmp: 

110 with h5py.File(tmp.name, "w") as f: 

111 arch = NdfOutputArchive(f) 

112 ref = arch.add_array(data, name="variance") 

113 self.assertEqual(ref.source, "ndf:/VARIANCE/DATA") 

114 with h5py.File(tmp.name, "r") as f: 

115 self.assertEqual(f["/VARIANCE"].attrs["CLASS"], b"ARRAY") 

116 self.assertEqual(f["/VARIANCE/DATA"].dtype, np.float64) 

117 

118 def test_top_level_compatible_mask_routes_to_quality(self): 

119 data = np.array([[0, 1, 2], [3, 4, 5]], dtype=np.uint8) 

120 with tempfile.NamedTemporaryFile(suffix=".sdf") as tmp: 

121 with h5py.File(tmp.name, "w") as f: 

122 arch = NdfOutputArchive(f) 

123 ref = arch.add_array(data, name="mask") 

124 self.assertEqual(ref.source, "ndf:/QUALITY/QUALITY/DATA") 

125 with h5py.File(tmp.name, "r") as f: 

126 self.assertEqual(f["/QUALITY"].attrs["CLASS"], b"QUALITY") 

127 self.assertEqual(f["/QUALITY/QUALITY"].attrs["CLASS"], b"ARRAY") 

128 self.assertEqual(f["/QUALITY/QUALITY/DATA"].dtype, np.uint8) 

129 np.testing.assert_array_equal(f["/QUALITY/QUALITY/DATA"][()], data) 

130 self.assertEqual(f["/QUALITY/QUALITY/ORIGIN"].dtype, np.int32) 

131 self.assertEqual(f["/QUALITY/QUALITY/ORIGIN"].shape, (2,)) 

132 self.assertEqual(f["/QUALITY/QUALITY/BAD_PIXEL"].id.get_type().get_class(), h5py.h5t.BITFIELD) 

133 self.assertFalse(_hds.read_array(f["/QUALITY/QUALITY/BAD_PIXEL"])) 

134 self.assertEqual(f["/QUALITY/BADBITS"][()], 255) 

135 

136 def test_top_level_incompatible_mask_routes_to_more_lsst(self): 

137 # 3D mask array in NDF storage order (mask-byte, y, x) is hoisted 

138 # as a sub-NDF inside /MORE/LSST/MASK, with a compressed 2D view 

139 # exposed as /QUALITY/QUALITY for standard NDF applications. 

140 data = np.zeros((2, 3, 4), dtype=np.uint8) 

141 data[0, 1, 2] = 4 

142 data[1, 2, 3] = 8 

143 expected_quality = np.any(data != 0, axis=0).astype(np.uint8) 

144 with tempfile.NamedTemporaryFile(suffix=".sdf") as tmp: 

145 with h5py.File(tmp.name, "w") as f: 

146 arch = NdfOutputArchive(f) 

147 ref = arch.add_array(data, name="mask") 

148 self.assertEqual(ref.source, "ndf:/MORE/LSST/MASK/DATA_ARRAY/DATA") 

149 with h5py.File(tmp.name, "r") as f: 

150 # /MORE/LSST/MASK is a real NDF: top-level CLASS="NDF" 

151 # containing a DATA_ARRAY structure with DATA + ORIGIN. 

152 self.assertEqual(f["/MORE/LSST/MASK"].attrs["CLASS"], b"NDF") 

153 self.assertEqual(f["/MORE/LSST/MASK/DATA_ARRAY"].attrs["CLASS"], b"ARRAY") 

154 self.assertEqual(f["/MORE/LSST/MASK/DATA_ARRAY/DATA"].shape, data.shape) 

155 self.assertEqual(f["/QUALITY/QUALITY"].attrs["CLASS"], b"ARRAY") 

156 np.testing.assert_array_equal(f["/QUALITY/QUALITY/DATA"][()], expected_quality) 

157 self.assertEqual(f["/QUALITY/BADBITS"][()], 255) 

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

159 self.assertEqual(origin.dtype, np.int64) 

160 self.assertEqual(origin.shape, (3,)) 

161 

162 def test_long_hoisted_component_is_shrunk(self): 

163 # Regression for the cell_coadd failure: the /noise_realizations/0 

164 # archive path contains an 18-character component. 

165 data = np.array([[1.0, 2.0]], dtype=np.float32) 

166 with tempfile.NamedTemporaryFile(suffix=".sdf") as tmp: 

167 with h5py.File(tmp.name, "w") as f: 

168 arch = NdfOutputArchive(f) 

169 ref = arch.add_array(data, name="noise_realizations/0") 

170 # The reported path is exactly what is stored in the JSON. 

171 self.assertTrue(ref.source.startswith("ndf:/MORE/LSST/")) 

172 self.assertTrue(ref.source.endswith("/DATA_ARRAY/DATA")) 

173 with h5py.File(tmp.name, "r") as f: 

174 # Every HDS component is within the limit. 

175 hdf5_path = ref.source[len("ndf:") :] 

176 for component in hdf5_path.strip("/").split("/"): 

177 self.assertLessEqual(len(component), DAT__SZNAM) 

178 # The node the JSON points at actually exists. 

179 self.assertIn(hdf5_path, f) 

180 

181 def test_long_name_round_trips_through_input_archive(self): 

182 from lsst.images.ndf import NdfInputArchive 

183 

184 data = np.arange(6, dtype=np.float32).reshape(2, 3) 

185 with tempfile.NamedTemporaryFile(suffix=".sdf") as tmp: 

186 with h5py.File(tmp.name, "w") as f: 

187 arch = NdfOutputArchive(f) 

188 ref = arch.add_array(data, name="noise_realizations/0") 

189 with NdfInputArchive.open(tmp.name) as inp: 

190 read_back = inp.get_array(ref) 

191 np.testing.assert_array_equal(read_back, data) 

192 

193 def test_repeated_long_name_gets_distinct_versioned_paths(self): 

194 data = np.array([[1.0]], dtype=np.float32) 

195 with tempfile.NamedTemporaryFile(suffix=".sdf") as tmp: 

196 with h5py.File(tmp.name, "w") as f: 

197 arch = NdfOutputArchive(f) 

198 first = arch.add_array(data, name="noise_realizations_value") 

199 second = arch.add_array(data, name="noise_realizations_value") 

200 self.assertNotEqual(first.source, second.source) 

201 # The second occurrence keeps a visible _2 version suffix. 

202 second_leaf = second.source[len("ndf:") :].split("/")[-3] 

203 self.assertTrue(second_leaf.endswith("_2")) 

204 with h5py.File(tmp.name, "r") as f: 

205 self.assertIn(first.source[len("ndf:") :], f) 

206 self.assertIn(second.source[len("ndf:") :], f) 

207 

208 def test_nested_array_hoists_as_sub_ndf(self): 

209 # Hoisted numeric arrays land under /MORE/LSST as hierarchical 

210 # sub-NDFs (CLASS="NDF" with DATA_ARRAY/DATA + ORIGIN inside) so 

211 # Starlink tools can inspect them as ordinary NDFs while each HDS 

212 # component stays short. 

213 data = np.array([[1.0, 2.0]], dtype=np.float32) 

214 with tempfile.NamedTemporaryFile(suffix=".sdf") as tmp: 

215 with h5py.File(tmp.name, "w") as f: 

216 arch = NdfOutputArchive(f) 

217 ref = arch.add_array(data, name="psf/coefficients") 

218 self.assertEqual(ref.source, "ndf:/MORE/LSST/PSF/COEFFICIENTS/DATA_ARRAY/DATA") 

219 with h5py.File(tmp.name, "r") as f: 

220 self.assertIn("MORE", f) 

221 self.assertIn("LSST", f["/MORE"]) 

222 self.assertIn("PSF", f["/MORE/LSST"]) 

223 self.assertIn("COEFFICIENTS", f["/MORE/LSST/PSF"]) 

224 sub = f["/MORE/LSST/PSF/COEFFICIENTS"] 

225 self.assertEqual(sub.attrs["CLASS"], b"NDF") 

226 self.assertEqual(sub["DATA_ARRAY"].attrs["CLASS"], b"ARRAY") 

227 np.testing.assert_array_equal(sub["DATA_ARRAY/DATA"][()], data) 

228 origin = sub["DATA_ARRAY/ORIGIN"] 

229 self.assertEqual(origin.dtype, np.int64) 

230 self.assertEqual(origin.shape, (data.ndim,)) 

231 

232 def test_colliding_shrunk_names_raise(self): 

233 data = np.array([[1.0]], dtype=np.float32) 

234 with tempfile.NamedTemporaryFile(suffix=".sdf") as tmp: 

235 with h5py.File(tmp.name, "w") as f: 

236 arch = NdfOutputArchive(f) 

237 # Force both long names to shrink to the same HDS token. 

238 with mock.patch.object( 

239 arch._name_shrinker, 

240 "shrink", 

241 side_effect=lambda name, *a, **k: name.upper() if len(name) <= DAT__SZNAM else "CLASH", 

242 ): 

243 arch.add_array(data, name="long_component_name_one") 

244 with self.assertRaisesRegex(ValueError, "name collision"): 

245 arch.add_array(data, name="long_component_name_two") 

246 

247 

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

249class NdfOutputArchivePointerTestCase(unittest.TestCase): 

250 """Tests for `NdfOutputArchive.serialize_pointer` and 

251 `serialize_frame_set`. 

252 """ 

253 

254 def test_serialize_pointer_writes_subtree_and_returns_pointer(self): 

255 with tempfile.NamedTemporaryFile(suffix=".sdf") as tmp: 

256 with h5py.File(tmp.name, "w") as f: 

257 arch = NdfOutputArchive(f) 

258 ptr = arch.serialize_pointer( 

259 "psf", 

260 lambda nested: TinyTree(name="gaussian"), 

261 key=("psf", 1), 

262 ) 

263 self.assertEqual(ptr.path, "/MORE/LSST/PSF/JSON") 

264 with h5py.File(tmp.name, "r") as f: 

265 # The hoisted sub-tree is stored as a "JSON" _CHAR*N 

266 # child of the target structure. 

267 raw = f["/MORE/LSST/PSF/JSON"][()] 

268 joined = b"".join(raw).decode("ascii").rstrip(" ") 

269 self.assertIn('"name":"gaussian"', joined.replace(" ", "")) 

270 

271 def test_serialize_pointer_caches_by_key(self): 

272 with tempfile.NamedTemporaryFile(suffix=".sdf") as tmp: 

273 with h5py.File(tmp.name, "w") as f: 

274 arch = NdfOutputArchive(f) 

275 ptr1 = arch.serialize_pointer( 

276 "psf", 

277 lambda nested: TinyTree(name="first"), 

278 key=("psf", 1), 

279 ) 

280 # Same key -> returns cached pointer; serializer not re-run 

281 # (we'd otherwise overwrite the file content with "second"). 

282 ptr2 = arch.serialize_pointer( 

283 "psf", 

284 lambda nested: TinyTree(name="second"), 

285 key=("psf", 1), 

286 ) 

287 self.assertEqual(ptr1, ptr2) 

288 with h5py.File(tmp.name, "r") as f: 

289 raw = f["/MORE/LSST/PSF/JSON"][()] 

290 joined = b"".join(raw).decode("ascii").rstrip(" ") 

291 self.assertIn("first", joined) 

292 self.assertNotIn("second", joined) 

293 

294 def test_serialize_pointer_preserves_nested_arrays(self): 

295 # Regression test: a pointer target that writes a nested array via 

296 # the nested archive must round-trip with that array still in the 

297 # file. Previously the pointer JSON was written at the target path 

298 # itself, clobbering any nested data the serializer produced. 

299 class TreeWithArray(pydantic.BaseModel): 

300 name: str 

301 data: ArrayReferenceModel 

302 

303 payload = np.arange(6, dtype=np.float32).reshape(2, 3) 

304 with tempfile.NamedTemporaryFile(suffix=".sdf") as tmp: 

305 with h5py.File(tmp.name, "w") as f: 

306 arch = NdfOutputArchive(f) 

307 ptr = arch.serialize_pointer( 

308 "psf", 

309 lambda nested: TreeWithArray( 

310 name="gaussian", 

311 data=nested.add_array(payload, name="parameters"), 

312 ), 

313 key=("psf", 1), 

314 ) 

315 self.assertEqual(ptr.path, "/MORE/LSST/PSF/JSON") 

316 with h5py.File(tmp.name, "r") as f: 

317 # JSON stayed at <path>/JSON; the nested array is still 

318 # accessible at the sub-NDF path the serializer wrote. 

319 self.assertIn("/MORE/LSST/PSF/JSON", f) 

320 self.assertIn("/MORE/LSST/PSF/PARAMETERS/DATA_ARRAY/DATA", f) 

321 np.testing.assert_array_equal(f["/MORE/LSST/PSF/PARAMETERS/DATA_ARRAY/DATA"][()], payload) 

322 # The pointer-target structure is typed after its leaf 

323 # name, not the generic EXT. 

324 self.assertEqual(f["/MORE/LSST/PSF"].attrs["CLASS"], b"PSF") 

325 self.assertEqual(f["/MORE/LSST"].attrs["CLASS"], b"LSST") 

326 

327 def test_serialize_frame_set_records_for_iter(self): 

328 # serialize_frame_set is delegated to serialize_pointer plus 

329 # recording the (FrameSet, pointer) pair for iter_frame_sets, 

330 # mirroring how FITS and JSON archives behave. 

331 frame_set = TinyFrameSet() 

332 with tempfile.NamedTemporaryFile(suffix=".sdf") as tmp: 

333 with h5py.File(tmp.name, "w") as f: 

334 arch = NdfOutputArchive(f) 

335 ptr = arch.serialize_frame_set( 

336 "wcs/pixel_to_sky", 

337 frame_set, 

338 lambda nested: TinyTree(name="proj"), 

339 key=("frame_set", 1), 

340 ) 

341 self.assertEqual(ptr.path, "/MORE/LSST/WCS/PIXEL_TO_SKY/JSON") 

342 recorded = list(arch.iter_frame_sets()) 

343 self.assertEqual(len(recorded), 1) 

344 self.assertIs(recorded[0][0], frame_set) 

345 self.assertEqual(recorded[0][1].path, "/MORE/LSST/WCS/PIXEL_TO_SKY/JSON") 

346 

347 

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

349class NdfOutputArchiveAddTableTestCase(unittest.TestCase): 

350 """Tests for `NdfOutputArchive.add_table` and `add_structured_array`.""" 

351 

352 def test_add_table_returns_inline_table_model(self): 

353 t = astropy.table.Table({"a": [1, 2, 3], "b": [4.0, 5.0, 6.0]}) 

354 with tempfile.NamedTemporaryFile(suffix=".sdf") as tmp: 

355 with h5py.File(tmp.name, "w") as f: 

356 arch = NdfOutputArchive(f) 

357 model = arch.add_table(t, name="some_table") 

358 self.assertEqual(len(model.columns), 2) 

359 # v1 stores tables inline in the JSON tree. 

360 self.assertIsInstance(model.columns[0].data, InlineArrayModel) 

361 

362 def test_add_structured_array_writes_column_ndfs_with_units(self): 

363 rec = np.zeros(3, dtype=[("x", np.float64), ("y", np.int32)]) 

364 rec["x"] = [1.0, 2.0, 3.0] 

365 rec["y"] = [10, 20, 30] 

366 with tempfile.NamedTemporaryFile(suffix=".sdf") as tmp: 

367 with h5py.File(tmp.name, "w") as f: 

368 arch = NdfOutputArchive(f) 

369 model = arch.add_structured_array( 

370 rec, 

371 name="rec", 

372 units={"x": u.m}, 

373 descriptions={"y": "the y values"}, 

374 ) 

375 self.assertEqual(len(model.columns), 2) 

376 self.assertIsInstance(model.columns[0].data, ArrayReferenceModel) 

377 # Confirm units/descriptions were applied. 

378 col_x = next(c for c in model.columns if c.name == "x") 

379 col_y = next(c for c in model.columns if c.name == "y") 

380 self.assertEqual(col_x.unit, u.m) 

381 self.assertEqual(col_y.description, "the y values") 

382 self.assertEqual(col_x.data.source, "ndf:/MORE/LSST/REC/X/DATA_ARRAY/DATA") 

383 self.assertEqual(col_y.data.source, "ndf:/MORE/LSST/REC/Y/DATA_ARRAY/DATA") 

384 with h5py.File(tmp.name, "r") as f: 

385 self.assertEqual(f["/MORE/LSST/REC/X"].attrs["CLASS"], b"NDF") 

386 np.testing.assert_array_equal(f["/MORE/LSST/REC/X/DATA_ARRAY/DATA"][()], rec["x"]) 

387 self.assertEqual(f["/MORE/LSST/REC/Y"].attrs["CLASS"], b"NDF") 

388 np.testing.assert_array_equal(f["/MORE/LSST/REC/Y/DATA_ARRAY/DATA"][()], rec["y"]) 

389 with NdfInputArchive.open(tmp.name) as archive: 

390 recovered = archive.get_structured_array(model) 

391 np.testing.assert_array_equal(recovered, rec) 

392 

393 def test_add_single_column_structured_array_uses_table_name(self): 

394 rec = np.zeros(1, dtype=[("solution", np.float64, (4,))]) 

395 rec["solution"] = [[1.0, 2.0, 3.0, 4.0]] 

396 with tempfile.NamedTemporaryFile(suffix=".sdf") as tmp: 

397 with h5py.File(tmp.name, "w") as f: 

398 arch = NdfOutputArchive(f) 

399 model = arch.add_structured_array(rec, name="psf/piff/interp/solution") 

400 self.assertEqual(len(model.columns), 1) 

401 column = model.columns[0] 

402 self.assertIsInstance(column.data, ArrayReferenceModel) 

403 self.assertEqual( 

404 column.data.source, 

405 "ndf:/MORE/LSST/PSF/PIFF/INTERP/SOLUTION/DATA_ARRAY/DATA", 

406 ) 

407 self.assertEqual(column.data.shape, [4]) 

408 with h5py.File(tmp.name, "r") as f: 

409 self.assertIn("PSF", f["/MORE/LSST"]) 

410 self.assertIn("PIFF", f["/MORE/LSST/PSF"]) 

411 self.assertIn("INTERP", f["/MORE/LSST/PSF/PIFF"]) 

412 self.assertIn("SOLUTION", f["/MORE/LSST/PSF/PIFF/INTERP"]) 

413 np.testing.assert_array_equal( 

414 f["/MORE/LSST/PSF/PIFF/INTERP/SOLUTION/DATA_ARRAY/DATA"][()], 

415 rec["solution"], 

416 ) 

417 

418 def test_structured_array_long_name_is_shrunk_and_versioned(self): 

419 dtype = np.dtype([("alpha", "f8"), ("beta", "i4")]) 

420 arr = np.zeros(3, dtype=dtype) 

421 with tempfile.NamedTemporaryFile(suffix=".sdf") as tmp: 

422 with h5py.File(tmp.name, "w") as f: 

423 arch = NdfOutputArchive(f) 

424 first = arch.add_structured_array(arr, name="catalog_of_long_named_sources") 

425 second = arch.add_structured_array(arr, name="catalog_of_long_named_sources") 

426 for model in (first, second): 

427 for column in model.columns: 

428 token = column.data.source[len("ndf:") :] 

429 for component in token.strip("/").split("/"): 

430 self.assertLessEqual(len(component), DAT__SZNAM) 

431 # The two structured arrays land in distinct sub-trees. 

432 self.assertNotEqual( 

433 first.columns[0].data.source, 

434 second.columns[0].data.source, 

435 ) 

436 # The second structured array's parent token keeps a 

437 # visible _2 version suffix. 

438 second_parent = second.columns[0].data.source[len("ndf:") :].strip("/").split("/")[-4] 

439 self.assertTrue(second_parent.endswith("_2")) 

440 with h5py.File(tmp.name, "r") as f: 

441 for model in (first, second): 

442 for column in model.columns: 

443 self.assertIn(column.data.source[len("ndf:") :], f) 

444 

445 

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

447class NdfWriteWcsTestCase(unittest.TestCase): 

448 """Tests for /WCS/DATA serialization in ndf.write().""" 

449 

450 def test_write_with_projection_creates_wcs_component(self): 

451 rng = np.random.default_rng(42) 

452 det_frame = DetectorFrame(instrument="TestInst", detector=4, bbox=Box.factory[1:4096, 1:4096]) 

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

454 sky_projection = make_random_sky_projection(rng, det_frame, Box.factory[1:4096, 1:4096]) 

455 image = Image( 

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

457 bbox=bbox, 

458 sky_projection=sky_projection, 

459 ) 

460 with tempfile.NamedTemporaryFile(suffix=".sdf", delete_on_close=False) as tmp: 

461 tmp.close() 

462 write(image, tmp.name) 

463 with h5py.File(tmp.name, "r") as f: 

464 self.assertIn("WCS", f) 

465 self.assertEqual(f["/WCS"].attrs["CLASS"], b"WCS") 

466 wcs_data = f["/WCS/DATA"] 

467 self.assertEqual(wcs_data.dtype, np.dtype("|S32")) 

468 records = [s.decode("ascii").rstrip(" ") for s in wcs_data[()]] 

469 self.assertTrue(all(record[0] in {" ", "+"} for record in records)) 

470 self.assertFalse(any(record.startswith("#") for record in records)) 

471 text = _hds.decode_ndf_ast_data(records) 

472 stripped = [line.lstrip() for line in text.splitlines()] 

473 self.assertTrue(any(s.startswith("Begin FrameSet") for s in stripped)) 

474 self.assertTrue(any(s.startswith("End FrameSet") for s in stripped)) 

475 self.assertIn('Domain = "GRID"', stripped) 

476 self.assertIn('Domain = "PIXEL"', stripped) 

477 self.assertIn("Sft1 = -19", stripped) 

478 self.assertIn("Sft2 = -9", stripped) 

479 

480 def test_write_without_projection_omits_wcs_component(self): 

481 # Image with no sky_projection -> no /WCS in the file. 

482 image = Image(np.zeros((2, 2), dtype=np.float32)) 

483 with tempfile.NamedTemporaryFile(suffix=".sdf", delete_on_close=False) as tmp: 

484 tmp.close() 

485 write(image, tmp.name) 

486 with h5py.File(tmp.name, "r") as f: 

487 self.assertNotIn("WCS", f) 

488 

489 def test_mask_sub_ndf_gets_3d_wcs(self): 

490 # When an incompatible mask is hoisted to /MORE/LSST/MASK as a 

491 # sub-NDF, it should carry its own 3D /WCS. The first two axes 

492 # retain the parent image sky_projection while the third axis is 

493 # a generic mask-byte coordinate. 

494 rng = np.random.default_rng(42) 

495 det_frame = DetectorFrame(instrument="TestInst", detector=4, bbox=Box.factory[1:4096, 1:4096]) 

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

497 sky_projection = make_random_sky_projection(rng, det_frame, Box.factory[1:4096, 1:4096]) 

498 # 12-plane schema -> native 3D uint8 mask, hoisted to /MORE/LSST/MASK. 

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

500 image = Image( 

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

502 bbox=bbox, 

503 sky_projection=sky_projection, 

504 ) 

505 masked = MaskedImage(image, mask_schema=MaskSchema(planes)) 

506 with tempfile.NamedTemporaryFile(suffix=".sdf", delete_on_close=False) as tmp: 

507 tmp.close() 

508 write(masked, tmp.name) 

509 with h5py.File(tmp.name, "r") as f: 

510 # Top-level WCS is present (existing behaviour). 

511 self.assertIn("WCS", f) 

512 top_lines = [s.decode("ascii") for s in f["/WCS/DATA"][()]] 

513 self.assertIn("MASK", f["/MORE/LSST"]) 

514 self.assertIn("WCS", f["/MORE/LSST/MASK"]) 

515 self.assertEqual(f["/MORE/LSST/MASK/WCS"].attrs["CLASS"], b"WCS") 

516 mask_lines = [s.decode("ascii") for s in f["/MORE/LSST/MASK/WCS/DATA"][()]] 

517 self.assertNotEqual(top_lines, mask_lines) 

518 mask_text = _hds.decode_ndf_ast_data(mask_lines) 

519 stripped = [line.lstrip() for line in mask_text.splitlines()] 

520 self.assertIn("Naxes = 3", stripped) 

521 self.assertIn('Domain = "GRID"', stripped) 

522 self.assertIn('Domain = "PIXEL"', stripped) 

523 self.assertIn("Sft1 = -19", stripped) 

524 self.assertIn("Sft2 = -9", stripped) 

525 self.assertIn("Sft3 = 1", stripped) 

526 self.assertIn("Begin CmpFrame", stripped) 

527 self.assertIn("Begin SkyFrame", stripped) 

528 self.assertIn('Domain = "MASK"', stripped) 

529 self.assertIn("Begin CmpMap", stripped) 

530 self.assertIn("Series = 0", stripped) 

531 

532 def test_mask_sub_ndf_no_wcs_when_image_has_no_projection(self): 

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

534 masked = MaskedImage( 

535 Image(np.zeros((4, 5), dtype=np.float32)), 

536 mask_schema=MaskSchema(planes), 

537 ) 

538 with tempfile.NamedTemporaryFile(suffix=".sdf", delete_on_close=False) as tmp: 

539 tmp.close() 

540 write(masked, tmp.name) 

541 with h5py.File(tmp.name, "r") as f: 

542 self.assertNotIn("WCS", f) 

543 self.assertIn("MASK", f["/MORE/LSST"]) 

544 self.assertNotIn("WCS", f["/MORE/LSST/MASK"]) 

545 

546 

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

548class NdfWriteFunctionTestCase(unittest.TestCase): 

549 """End-to-end tests for the module-level `write()` function.""" 

550 

551 def test_write_image_produces_valid_layout(self): 

552 image = Image( 

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

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

555 ) 

556 with tempfile.NamedTemporaryFile(suffix=".sdf", delete_on_close=False) as tmp: 

557 tmp.close() 

558 tree = write(image, tmp.name) 

559 self.assertIsNotNone(tree) 

560 with h5py.File(tmp.name, "r") as f: 

561 # Root is an NDF with a name. 

562 self.assertEqual(f["/"].attrs["CLASS"], b"NDF") 

563 self.assertIn("HDS_ROOT_NAME", f["/"].attrs) 

564 # DATA_ARRAY uses the complex form (DATA + ORIGIN). 

565 self.assertEqual(f["/DATA_ARRAY"].attrs["CLASS"], b"ARRAY") 

566 np.testing.assert_array_equal(f["/DATA_ARRAY/DATA"][()], image.array) 

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

568 self.assertEqual(origin.dtype, np.int64) 

569 self.assertEqual(len(origin), 2) 

570 # ORIGIN encodes bbox lower bounds in Fortran order. The exact 

571 # values depend on Box's API; just verify it isn't the 

572 # all-zeros placeholder when the bbox is non-trivial. 

573 self.assertFalse((origin == 0).all()) 

574 # Main JSON tree at /MORE/LSST/JSON. 

575 self.assertIn("MORE", f) 

576 self.assertIn("LSST", f["/MORE"]) 

577 self.assertIn("JSON", f["/MORE/LSST"]) 

578 

579 def test_write_image_preserves_opaque_fits_metadata(self): 

580 image = Image(np.zeros((2, 2), dtype=np.float32)) 

581 # Attach an opaque-metadata primary header to the image. 

582 primary = astropy.io.fits.Header() 

583 primary["FOO"] = ("bar", "test card") 

584 long_value = "x" * 100 

585 primary["LONGSTR"] = (long_value, "long string value") 

586 opaque = FitsOpaqueMetadata() 

587 opaque.add_header(primary, name="", ver=1) 

588 image._opaque_metadata = opaque 

589 with tempfile.NamedTemporaryFile(suffix=".sdf", delete_on_close=False) as tmp: 

590 tmp.close() 

591 write(image, tmp.name) 

592 with h5py.File(tmp.name, "r") as f: