Coverage for python/lsst/images/ndf/_hds.py: 89%
132 statements
« prev ^ index » next coverage.py v7.14.2, created at 2026-06-21 08:57 +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.
12"""HDS-on-HDF5 read/write helpers.
14These follow the conventions used by the canonical Starlink ``hds-v5``
15library (see ``reference/hds-v5/dat1.h`` and ``dat1New.c``):
17* HDS structures are HDF5 groups with a ``CLASS`` attribute holding the
18 HDS type string (``"NDF"``, ``"WCS"``, ``"EXT"``, ``"ARRAY"``, ...).
19* Arrays of structures additionally carry an ``HDS_STRUCTURE_DIMS``
20 attribute (deferred from v1; we only handle scalar structures).
21* The file root group, when it represents a top-level HDS structure,
22 carries an ``HDS_ROOT_NAME`` attribute giving the HDS object name.
23* HDS primitives are bare HDF5 datasets with no HDS-specific attributes;
24 the HDS type is inferred from the HDF5 dtype.
25"""
27from __future__ import annotations
29from collections.abc import Iterator, Sequence
30from typing import Any
32import h5py
33import numpy as np
35__all__ = (
36 "ATTR_CLASS",
37 "ATTR_ROOT_NAME",
38 "ATTR_STRUCTURE_DIMS",
39 "HDS_TO_NUMPY",
40 "NUMPY_TO_HDS",
41 "create_structure",
42 "decode_ndf_ast_data",
43 "encode_ndf_ast_data",
44 "hds_type_for_dtype",
45 "iter_children",
46 "open_structure",
47 "read_array",
48 "read_char_array",
49 "set_ascii_attr",
50 "set_root_name",
51 "write_array",
52 "write_char_array",
53)
56# Canonical attribute names used by hds-v5 (see reference/hds-v5/dat1.h).
57ATTR_CLASS = "CLASS"
58ATTR_STRUCTURE_DIMS = "HDS_STRUCTURE_DIMS"
59ATTR_ROOT_NAME = "HDS_ROOT_NAME"
62HDS_TO_NUMPY: dict[str, np.dtype] = {
63 "_LOGICAL": np.dtype(np.bool_),
64 "_REAL": np.dtype(np.float32),
65 "_DOUBLE": np.dtype(np.float64),
66 "_UBYTE": np.dtype(np.uint8),
67 "_WORD": np.dtype(np.int16),
68 "_INTEGER": np.dtype(np.int32),
69 "_INT64": np.dtype(np.int64),
70}
72NUMPY_TO_HDS: dict[np.dtype, str] = {
73 np.dtype(np.bool_): "_LOGICAL",
74 np.dtype(np.float32): "_REAL",
75 np.dtype(np.float64): "_DOUBLE",
76 np.dtype(np.uint8): "_UBYTE",
77 np.dtype(np.int16): "_WORD",
78 np.dtype(np.int32): "_INTEGER",
79 np.dtype(np.int64): "_INT64",
80}
83NDF_AST_DATA_WIDTH = 32
84NDF_AST_DATA_MIN_WIDTH = 16
86# HDS object-name length limit, from the Starlink DAT_PAR include
87# (``dat_par.h``): ``DAT__SZNAM 15`` ("Size of object name").
88# The sibling limits ``DAT__SZGRP`` (group name) and ``DAT__SZTYP`` (type
89# string) are also 15 today; only the object-name limit is enforced here,
90# because HDS type tags are derived from already-shrunk component names.
91DAT__SZNAM = 15
94def hds_type_for_dtype(dtype: np.dtype) -> str:
95 """Return the HDS type string for a numpy dtype.
97 Parameters
98 ----------
99 dtype
100 Numpy dtype to map to an HDS primitive type.
102 Returns
103 -------
104 str
105 HDS primitive type string.
107 Raises
108 ------
109 NotImplementedError
110 Raised if ``dtype`` does not map to a supported HDS primitive type.
112 Notes
113 -----
114 Fixed-width byte strings ``|S<N>`` map to ``"_CHAR*<N>"``. Numeric
115 dtypes are looked up in `NUMPY_TO_HDS`. Anything else raises
116 ``NotImplementedError``.
117 """
118 if dtype.kind == "S":
119 return f"_CHAR*{dtype.itemsize}"
120 try:
121 return NUMPY_TO_HDS[dtype]
122 except KeyError:
123 raise NotImplementedError(f"No HDS type mapping for dtype {dtype!r}.") from None
126def write_array(
127 parent: h5py.Group,
128 name: str,
129 data: np.ndarray,
130 *,
131 compression: str | None = None,
132 compression_opts: int | None = None,
133) -> h5py.Dataset:
134 """Write a numpy C-order array as an HDS primitive.
136 Parameters
137 ----------
138 parent
139 HDF5 group to receive the dataset.
140 name
141 Name of the primitive dataset to create.
142 data
143 Array to write.
144 compression
145 Optional compression algorithm accepted by `h5py.Group.create_dataset`.
146 compression_opts
147 Optional compression options accepted by `h5py.Group.create_dataset`.
149 Returns
150 -------
151 h5py.Dataset
152 Newly-created dataset.
154 Notes
155 -----
156 The HDF5 dataset carries no HDS-specific attributes; the HDS type
157 is inferred on read from the HDF5 dtype. Refuses dtypes that don't
158 map to a supported HDS primitive type.
160 The HDF5 dataset has the array's natural shape (C-order). Combined
161 with HDF5's native byte ordering, this matches the Fortran-on-disk
162 layout required by HDS for an NDF whose Fortran-order shape is the
163 reverse of ``data.shape``.
164 """
165 # Validate the dtype is supported up front so callers get a clear error.
166 hds_type_for_dtype(data.dtype)
167 if data.dtype == np.dtype(np.bool_):
168 return _write_logical_array(parent, name, data, compression=compression)
169 # h5py rejects compression options if no compression algorithm is set.
170 kwargs: dict[str, Any] = {}
171 if compression is not None: 171 ↛ 172line 171 didn't jump to line 172 because the condition on line 171 was never true
172 kwargs["compression"] = compression
173 if compression_opts is not None: 173 ↛ 174line 173 didn't jump to line 174 because the condition on line 173 was never true
174 kwargs["compression_opts"] = compression_opts
175 return parent.create_dataset(name, data=data, **kwargs)
178def _write_logical_array(
179 parent: h5py.Group,
180 name: str,
181 data: np.ndarray,
182 *,
183 compression: str | None = None,
184) -> h5py.Dataset:
185 """Write an HDS ``_LOGICAL`` primitive using the HDF5 bitfield type.
187 Parameters
188 ----------
189 parent
190 HDF5 group to receive the dataset.
191 name
192 Name of the primitive dataset to create.
193 data
194 Boolean array to write.
195 compression
196 Compression is not supported for this low-level bitfield writer.
198 Returns
199 -------
200 h5py.Dataset
201 Newly-created dataset.
203 Notes
204 -----
205 High-level h5py writes numpy bool data as an HDF5 enum, but hds-v5
206 identifies ``_LOGICAL`` primitives by the HDF5 bitfield class.
207 """
208 if compression is not None:
209 raise NotImplementedError("Compression is not implemented for HDS _LOGICAL arrays.")
210 logical_data = np.asarray(data, dtype=np.uint8)
211 if logical_data.shape:
212 space = h5py.h5s.create_simple(logical_data.shape)
213 else:
214 space = h5py.h5s.create(h5py.h5s.SCALAR)
215 dataset_id = h5py.h5d.create(
216 parent.id,
217 name.encode("ascii"),
218 h5py.h5t.STD_B8LE,
219 space,
220 )
221 dataset_id.write(
222 h5py.h5s.ALL,
223 h5py.h5s.ALL,
224 logical_data,
225 mtype=h5py.h5t.NATIVE_B8,
226 )
227 dataset_id.close()
228 return parent[name]
231def read_array(dataset: h5py.Dataset) -> np.ndarray:
232 """Read an HDS primitive into a C-order numpy array.
234 Parameters
235 ----------
236 dataset
237 HDF5 dataset to read.
239 Returns
240 -------
241 numpy.ndarray
242 Data read from ``dataset``.
244 Notes
245 -----
246 The HDS type is inferred from the HDF5 dtype. Raises
247 ``NotImplementedError`` if the dtype is not a supported numeric HDS
248 primitive type. Use `read_char_array` for ``_CHAR*N`` datasets.
249 """
250 if dataset.dtype.kind == "S":
251 raise ValueError(f"Dataset {dataset.name!r} is _CHAR*N; use read_char_array instead.")
252 dataset_type = dataset.id.get_type()
253 if dataset_type.get_class() == h5py.h5t.BITFIELD:
254 if dataset_type.get_size() not in {1, 4}:
255 raise NotImplementedError(
256 f"Dataset {dataset.name!r} has bitfield size {dataset_type.get_size()} "
257 "which does not map to HDS _LOGICAL."
258 )
259 data = dataset[()] != 0
260 if isinstance(data, np.ndarray):
261 return data.astype(np.bool_)
262 return np.atleast_1d(np.bool_(data))
263 if dataset.dtype not in NUMPY_TO_HDS:
264 raise NotImplementedError(
265 f"Dataset {dataset.name!r} has dtype {dataset.dtype} which does not "
266 f"map to a supported HDS primitive type."
267 )
268 return dataset[()]
271def write_char_array(
272 parent: h5py.Group,
273 name: str,
274 lines: Sequence[str],
275 *,
276 width: int = 80,
277) -> h5py.Dataset:
278 """Write a sequence of strings as a 1D HDS ``_CHAR*N`` primitive.
280 Parameters
281 ----------
282 parent
283 HDF5 group to receive the dataset.
284 name
285 Name of the primitive dataset to create.
286 lines
287 Strings to write. Each string must be ASCII and no longer than
288 ``width`` characters.
289 width
290 Fixed width of each HDS character-array element.
292 Returns
293 -------
294 h5py.Dataset
295 Newly-created dataset.
297 Raises
298 ------
299 ValueError
300 Raised if any string is not ASCII or is longer than ``width``.
302 Notes
303 -----
304 Each string is padded to ``width`` with trailing spaces (HDS
305 convention). The HDF5 dataset has dtype ``|S<width>``; no HDS-specific
306 attributes are written.
307 """
308 encoded_lines: list[bytes] = []
309 for n, line in enumerate(lines):
310 try:
311 encoded_line = line.encode("ascii")
312 except UnicodeEncodeError as err:
313 raise ValueError(f"Line {n} for {name!r} is not ASCII.") from err
314 if len(encoded_line) > width:
315 raise ValueError(
316 f"Line {n} for {name!r} is {len(encoded_line)} bytes, longer than width={width}."
317 )
318 encoded_lines.append(encoded_line.ljust(width))
319 encoded = np.array(
320 encoded_lines,
321 dtype=f"|S{width}",
322 )
323 return parent.create_dataset(name, data=encoded)
326def encode_ndf_ast_data(text: str, *, width: int = NDF_AST_DATA_WIDTH) -> list[str]:
327 """Encode AST Channel text for an NDF ``WCS.DATA`` component.
329 Parameters
330 ----------
331 text
332 Text emitted by an AST Channel.
333 width
334 Fixed width of the target NDF ``WCS.DATA`` records.
336 Returns
337 -------
338 list [str]
339 HDS character-array records.
341 Notes
342 -----
343 Starlink NDF stores each AST text line in one or more fixed-width
344 ``_CHAR*32`` records. The first character of each record is a flag:
345 a space starts a new AST line and ``+`` continues the previous one.
346 The payload is the AST text line with leading indentation removed.
347 """
348 if width < NDF_AST_DATA_MIN_WIDTH: 348 ↛ 349line 348 didn't jump to line 349 because the condition on line 348 was never true
349 raise ValueError(
350 f"NDF AST DATA record width {width} is too short; minimum is {NDF_AST_DATA_MIN_WIDTH}."
351 )
353 records: list[str] = []
354 payload_width = width - 1
355 for raw_line in text.splitlines():
356 line = raw_line.lstrip(" ").rstrip(" ")
357 if not line: 357 ↛ 358line 357 didn't jump to line 358 because the condition on line 357 was never true
358 continue
359 for start in range(0, len(line), payload_width):
360 flag = " " if start == 0 else "+"
361 records.append(f"{flag}{line[start : start + payload_width]}")
362 return records
365def decode_ndf_ast_data(records: Sequence[str]) -> str:
366 """Decode an NDF ``WCS.DATA`` component into AST Channel text.
368 Parameters
369 ----------
370 records
371 HDS character-array records from ``WCS.DATA``.
373 Returns
374 -------
375 str
376 AST Channel text.
378 Notes
379 -----
380 This reverses `encode_ndf_ast_data`. If the input does not look like
381 NDF AST records, it is treated as plain AST Channel text for backward
382 compatibility with earlier non-canonical files.
383 """
384 if not records: 384 ↛ 385line 384 didn't jump to line 385 because the condition on line 384 was never true
385 return ""
386 if any(record and record[0] not in {" ", "+"} for record in records): 386 ↛ 387line 386 didn't jump to line 387 because the condition on line 386 was never true
387 return "\n".join(records) + "\n"
389 lines: list[str] = []
390 current: list[str] = []
391 for record in records:
392 if not record: 392 ↛ 393line 392 didn't jump to line 393 because the condition on line 392 was never true
393 continue
394 flag = record[0]
395 payload = record[1:]
396 if flag == "+":
397 if current: 397 ↛ 400line 397 didn't jump to line 400 because the condition on line 397 was always true
398 current.append(payload)
399 else:
400 current = [payload]
401 else:
402 if current:
403 lines.append("".join(current).rstrip(" "))
404 current = [payload]
405 if current: 405 ↛ 407line 405 didn't jump to line 407 because the condition on line 405 was always true
406 lines.append("".join(current).rstrip(" "))
407 return "\n".join(lines) + ("\n" if lines else "")
410def read_char_array(dataset: h5py.Dataset) -> list[str]:
411 """Read an HDS ``_CHAR*N`` 1D primitive as a list of stripped strings.
413 Parameters
414 ----------
415 dataset
416 HDF5 dataset to read.
418 Returns
419 -------
420 list [str]
421 Dataset elements decoded from ASCII with trailing spaces stripped.
423 Notes
424 -----
425 Validates the dataset has a fixed-width byte-string dtype (``|S<N>``).
426 """
427 if dataset.dtype.kind != "S":
428 raise ValueError(f"Dataset {dataset.name!r} is not _CHAR*N (dtype {dataset.dtype}).")
429 if dataset.ndim == 0: 429 ↛ 430line 429 didn't jump to line 430 because the condition on line 429 was never true
430 raise ValueError(f"Dataset {dataset.name!r} is a scalar _CHAR*N; only 1-D arrays are supported.")
431 raw = dataset[()]
432 return [item.decode("ascii").rstrip(" ") for item in raw]
435def set_ascii_attr(target: h5py.Group | h5py.Dataset, name: str, value: str) -> None:
436 """Write a fixed-length ASCII byte attribute.
438 Parameters
439 ----------
440 target
441 HDF5 object whose attribute should be set.
442 name
443 Attribute name.
444 value
445 ASCII value to write.
447 Notes
448 -----
449 Canonical ``hds-v5`` stores ``CLASS`` and ``HDS_ROOT_NAME`` as
450 fixed-length ASCII byte strings (e.g. ``|S5`` for ``"ARRAY"``).
451 h5py's default for Python ``str`` is variable-length UTF-8, which
452 Starlink tools (KAPPA, ``hdstrace``) can't decode — they show
453 garbage in the type-tag column. Writing as fixed-length bytes
454 matches the canonical layout.
455 """
456 encoded = value.encode("ascii")
457 if name in target.attrs:
458 del target.attrs[name]
459 target.attrs.create(name, encoded, dtype=f"|S{len(encoded)}")
462def create_structure(parent: h5py.Group, name: str, hds_type: str) -> h5py.Group:
463 """Create a named HDS structure (h5py group with ``CLASS`` attribute).
465 Parameters
466 ----------
467 parent
468 Group to create the new structure under.
469 name
470 Component name (HDS rules apply: uppercase letters/digits/underscores,
471 max ``DAT__SZNAM`` characters; not enforced here).
472 hds_type
473 HDS type string for the new structure (e.g. ``"NDF"``, ``"WCS"``,
474 ``"ARRAY"``, ``"EXT"``).
475 """
476 group = parent.create_group(name)
477 set_ascii_attr(group, ATTR_CLASS, hds_type)
478 return group
481def set_root_name(file: h5py.File, hds_name: str, hds_type: str) -> None:
482 """Mark a file's root group as a top-level HDS structure.
484 Parameters
485 ----------
486 file
487 HDF5 file whose root group represents the HDS root object.
488 hds_name
489 HDS root object name.
490 hds_type
491 HDS type string for the root object.
493 Notes
494 -----
495 Sets ``HDS_ROOT_NAME`` (the HDS object name) and ``CLASS`` (the HDS
496 type) on the root group, matching what ``hds-v5`` writes for a root
497 structure created via :c:func:`dat1New`.
498 """
499 set_ascii_attr(file["/"], ATTR_ROOT_NAME, hds_name)
500 set_ascii_attr(file["/"], ATTR_CLASS, hds_type)
503def open_structure(parent: h5py.Group, name: str) -> tuple[h5py.Group, str]:
504 """Open a child structure by name. Returns ``(group, hds_type)``.
506 Parameters
507 ----------
508 parent
509 HDF5 group containing the structure.
510 name
511 Name of the child structure.
513 Returns
514 -------
515 group : h5py.Group
516 Opened HDF5 group.
517 hds_type : str
518 HDS type string from the structure's ``CLASS`` attribute.
520 Notes
521 -----
522 Raises ``ValueError`` if the child is not a group, or has no
523 ``CLASS`` attribute. Accepts the legacy ``HDSTYPE`` attribute name
524 as a fallback so files written by older HDS variants can still be
525 inspected.
526 """
527 obj = parent[name]
528 if not isinstance(obj, h5py.Group): 528 ↛ 529line 528 didn't jump to line 529 because the condition on line 528 was never true
529 raise ValueError(f"{parent.name}/{name} is a dataset, not a structure.")
530 hds_type = obj.attrs.get(ATTR_CLASS)
531 if hds_type is None:
532 # Legacy fallback for older HDS-on-HDF5 variants.
533 hds_type = obj.attrs.get("HDSTYPE")
534 if isinstance(hds_type, bytes):
535 hds_type = hds_type.decode("ascii")
536 if not isinstance(hds_type, str):
537 raise ValueError(f"Group {obj.name!r} has no {ATTR_CLASS!r} (or legacy HDSTYPE) attribute.")
538 return obj, hds_type
541def iter_children(group: h5py.Group) -> Iterator[tuple[str, h5py.Group | h5py.Dataset]]:
542 """Iterate over a structure's direct children.
544 Parameters
545 ----------
546 group
547 HDF5 group to inspect.
549 Yields ``(name, child)`` pairs where ``child`` is a group or dataset.
550 """
551 yield from group.items()