Merge branch 'feat/p2-python-bindings' into feat/p2-perf-coverage

# Conflicts:
#	CHANGELOG.md
#	README.md

.gitea/workflows/ci.yml

@@ -38,7 +38,9 @@ jobs:
           # (clawhdf5-tools) interop tests compare against.
           apt-get install -y --no-install-recommends python3 python3-venv cmake hdf5-tools
           python3 -m venv /opt/interop
-          /opt/interop/bin/pip install --no-cache-dir h5py numpy netCDF4 xarray hdf5plugin
+          # maturin + pytest: ci-test.sh builds the Python package
+          # (crates/clawhdf5-py) and runs its tests against h5py.
+          /opt/interop/bin/pip install --no-cache-dir h5py numpy netCDF4 xarray hdf5plugin maturin pytest
           echo "/opt/interop/bin" >> "$GITHUB_PATH"
       - name: Show interop library versions
         # h5dump's version too: the h5rs dump test requires its exact output

.gitignore

@@ -5,3 +5,5 @@ benchmarks/longmemeval/*.json
 # Local model weights (MiniLM etc.) — large, not committed
 weights/
 .venv
+__pycache__/
+.pytest_cache/

CHANGELOG.md

@@ -230,6 +230,119 @@
   are not byte-identical to earlier versions). Regression test:
   `crates/clawhdf5/tests/writer_groups_interop.rs`.
 
+### Python bindings (2026-09-26)
+- **Panic: selections of v4 implicit-index chunked datasets** (pre-existing,
+  facade `Dataset::read_selection`, Rust callers too). A hyperslab whose
+  bounding box covered more than half of a chunked dataset with the implicit
+  index (`libver='latest'`, early allocation, no filters) panicked with
+  "index out of bounds" in `generate_implicit_chunks`: the fallback in
+  `data_read::read_raw_data_selection` passed the layout's chunk dimensions,
+  element-size dimension included, and then decoded the whole dataset
+  anyway. That arm now decodes and extracts directly, for every chunk index.
+  `crates/clawhdf5/tests/v4_chunk_index_selection.rs` reads small and large
+  hyperslabs of all five v4 indexes (single chunk, implicit, fixed array,
+  extensible array, B-tree v2) and compares them with h5py; it panicked
+  before. The Python bindings made this easy to reach (`ds[0:3]` on
+  libhdf5's `h5fc_ext*.h5` test files).
+- **`pip install` / `maturin develop` now gives `import clawhdf5`.** The
+  distribution in `crates/clawhdf5-py/pyproject.toml` was still called
+  `rustyhdf5` while the extension module was `clawhdf5`, and the package's
+  tests imported `rustyhdf5`, so they failed at collection. Distribution,
+  module and tests now all say `clawhdf5`, and the module has
+  `__version__`.
+- **h5py-style reads that read the selection, not the dataset.** `ds[...]`
+  used to read the whole dataset and slice it in numpy, and knew six
+  dtypes. Now
+  integers (negative from the end), slices with positive steps, `...`, one
+  increasing list of integers per key and compound field names map onto the
+  facade's hyperslab selection (a list is read one group of neighbouring
+  chunks at a time and picked from in memory), with h5py's results (numpy scalar for an all-integer
+  key, 0-d array for `scalar[...]`) and h5py's errors for everything else
+  (negative steps, `None`, boolean masks, out-of-range indices).
+  `Dataset.dtype` is the numpy dtype h5py reports, for every integer and
+  IEEE float width (incl. `float16`) in either byte order, `bool`, enums
+  (base integer with `metadata['enum']`), complex (`r`/`i` compounds),
+  fixed strings (`S<n>`), variable-length strings (`object` of `bytes`, as
+  h5py), variable-length sequences (`object` of arrays), opaque (`V<n>`),
+  HDF5 array types and compounds (numpy structured, offsets and padding
+  kept, nested). The bytes the library returns become the numpy array's
+  buffer without a copy. Types the mapping cannot describe exactly
+  (references, bitfields, time, non-IEEE floats, integers with padding
+  bits, variable-length members inside compounds) raise `TypeError` rather
+  than return guessed data. Attributes come back as h5py returns them
+  (numpy scalars and arrays with the stored dtype, `str` for
+  variable-length strings, `numpy.bytes_` for fixed ones — **a change**:
+  string attributes written by this package are fixed-length and used to
+  come back as `str` — and `clawhdf5.Empty` for a null dataspace, which
+  datasets return too). `Group`/`File` gain `get`, `values`, `items`,
+  iteration, `len`, `name`, absolute and relative paths (`g['/a/b']`,
+  `g['c/d']`, `f['/']`); `Dataset` gains `ndim`, `size`, `maxshape`,
+  `name`, `len()` and `numpy.asarray(ds)`. File access and decoding run
+  with the GIL released, so Python threads read in parallel.
+  `crates/clawhdf5-py/tests/test_read_vs_h5py.py` compares every read with
+  h5py 3.16 (HDF5 2.0) on a file h5py writes. One difference is h5py's:
+  it returns variable-length sequences of big-endian floats unswapped; this
+  package returns the stored values.
+- **A panic in the library is an ordinary Python exception.** PyO3 turns a
+  Rust panic into `PanicException`, a `BaseException` that `except
+  Exception` does not catch. Every call from the bindings into the library
+  is now guarded and a panic becomes `clawhdf5.InternalError` (a
+  `RuntimeError`) naming the object; with the implicit-index panic above
+  restored, `ds[0:30]` raises it.
+- **Wrong data: uninitialised padding in compound results of index lists.**
+  `ds[[0, 3, 6]]` joined one read per run with `np.concatenate`, which
+  copies structured dtypes field by field into an `np.empty` result, so the
+  padding bytes held whatever was in memory (pointers were seen) and leaked
+  through `tobytes()`, hashes and write-backs. The runs' bytes are now joined
+  in Rust, whole elements at a time, so the result carries the bytes read
+  from the file (h5py's, zero for files it wrote) and stays zero-copy.
+  The h5py comparisons now also compare every byte of structured values
+  (`test_compound_padding_bytes_match_h5py` and `assert_same`).
+- **Index lists no longer decode the same chunks once per run.** A list
+  index was one uncached hyperslab read per run of consecutive indices, so
+  on a chunked, compressed dataset every run decoded its chunk again:
+  `d[list(range(0, 200000, 40))]` over 20 gzip chunks took 8 s (h5py:
+  0.014 s). The list is now read in groups — for a chunked dataset a group
+  ends only where a whole chunk holds no selected index, so each chunk is
+  decoded once; otherwise at a gap of more than 64 KiB — and the selected
+  rows are picked from each group in Rust. The same read now takes 3.8 ms
+  (h5py 4.1 ms; release build on tank, best of 5).
+  `test_a_long_index_list_decodes_each_chunk_once` compares 1-D, 2-D and
+  contiguous cases with h5py under a 2 s bound (5.8 s before, debug build).
+- **Groups and datasets remember where they are.** Every `ds[...]`, and
+  every `g[k]`, resolved its path from the root again (two or three times
+  per open), and in a large group each resolution scans the group's links,
+  so visiting a group was quadratic: 4000 scalar datasets in one group took
+  39 s (`libver='earliest'`) and 131 s (`'latest'`) to list, read and
+  re-read in `test_big_groups_are_not_quadratic`; now 0.3 s each (debug
+  build). A `Dataset` keeps its object's address, and a `Group` (and the
+  file's root) its address and, once listed, its link table. New facade
+  API: `File::dataset_at(address)` opens a dataset without resolving a
+  path. libhdf5's `h5stat_newgrat.h5` (35001 members in the root): listing
+  takes 0.03 s and 2000 opens 1 ms (h5py: 0.022 s).
+- **`ds[np.array(1)]` is an integer index**, as in h5py; a 0-d integer
+  array went down the index-list path and raised a confusing `TypeError`.
+  The h5py comparison keys now include 0-d arrays on every axis.
+- **Tests that would notice a held GIL, and our extra errors.**
+  `test_reads_release_the_gil` times a Python thread spinning while another
+  reads: with the read made to hold the GIL it stalls for the whole read
+  (0.062 s of a 0.064 s read) and the test fails; released, its longest
+  stall is about 3 ms. (The existing threads test only checked values.)
+  `test_errors_match_h5py` now also requires that every key h5py reads
+  reads here too, with the same result, and covers more keys (0-d arrays,
+  repeated and empty lists, `()`, `...`).
+- **Docs say when a selection reads more than itself.** The README and
+  the package README said `ds[...]` reads only the selected elements,
+  without condition. The library decodes the whole dataset when the
+  selection's bounding box covers more than half of it, and for compact,
+  virtual, unwritten and non-default-fill chunked datasets; the READMEs,
+  the facade's `read_selection` docs and `docs/known-issues.md` now say so.
+- **CI builds and tests the Python package.** It was excluded from CI.
+  `scripts/ci-test.sh` now lints `clawhdf5-py`, builds the wheel with
+  maturin, unpacks it under `target/` and runs the pytest suite; skipped
+  without maturin/pytest in `$CLAWHDF5_PYTHON`, a failure then under
+  `CLAWHDF5_REQUIRE_INTEROP=1`. The CI interop venv installs both.
+
 ### Plugin filters (2026-09-26)
 - **LZF, bitshuffle, bzip2 and Blosc read and write, in pure Rust.** Files
   written by h5py with `compression="lzf"`, or with hdf5plugin's

README.md

@@ -73,8 +73,9 @@ breaking change, are in [CHANGELOG.md](CHANGELOG.md).
 - Default fusion weights are now the measured 0.4 / 0.6 (v2.5.0). Re-ranking had
   been discarding the retrieval score, costing the Markdown backend 40.6pp of
   Hit@1; fixed in v2.6.0.
-- Selection reads decode only the chunks they touch (a 64×64 window: 105 ms to
-  0.39 ms), and full reads are 1.2–1.9× faster (v2.5.0).
+- Selection reads whose bounding box covers at most half the dataset decode
+  only the chunks they touch (a 64×64 window: 105 ms to 0.39 ms), and full
+  reads are 1.2–1.9× faster (v2.5.0).
 
 **Memory**
 - A loaded store holds ~30% less (embeddings stored once, v2.6.0), and the
@@ -432,6 +433,54 @@ b.write("groups.h5")?;
 A group holds at most 65 535 links; more is an error, as is a link over
 65 515 bytes (a very long soft-link target) in a group of more than 8 links.
 
+### Python
+
+`crates/clawhdf5-py` is a Python package (PyO3 + numpy) that reads HDF5 with
+an h5py-shaped API and no libhdf5. It is not on PyPI; build it with
+[maturin](https://www.maturin.rs) into a virtualenv:
+
+```bash
+python -m venv .venv && . .venv/bin/activate
+pip install maturin numpy
+maturin develop --release -m crates/clawhdf5-py/Cargo.toml
+python -c "import clawhdf5; print(clawhdf5.__version__)"
+```
+
+```python
+import numpy as np
+import clawhdf5
+
+with clawhdf5.File("data.h5", "r") as f:
+    print(list(f.keys()))              # sorted member names, like h5py
+    ds = f["group/temperatures"]       # relative or absolute ("/group/...") paths
+    print(ds.shape, ds.dtype)          # dtype is the numpy dtype h5py reports
+    block = ds[100:200, ::4]           # a small selection reads only its chunks
+    row = ds[-1]                       # integers drop the axis
+    picked = ds[[1, 5, 9], :]          # one increasing index list per key
+    units = ds.attrs["units"]          # attributes come back as h5py returns them
+    everything = np.asarray(ds)
+
+    records = f["table"]               # compound -> numpy structured array
+    ids = records["id"]                # one field
+```
+
+Reads cover integers and IEEE floats of every width in either byte order,
+`bool`, enums, complex, fixed and variable-length strings, variable-length
+sequences, opaque, HDF5 array types and compounds; other types (references,
+bitfields, ...) raise `TypeError` instead of returning guessed data. Keys
+follow h5py (negative steps, `None` and boolean masks are refused). The
+read itself runs with the GIL released, so Python threads read in parallel.
+A selection whose bounding box covers at most half the dataset decodes only
+the chunks (or contiguous rows) that box overlaps; a larger one — including
+a strided slice across the whole dataset — decodes the whole dataset, as
+do datasets that are compact, virtual, unwritten, or chunked with a
+non-default fill value (`docs/known-issues.md`). An index list is read one
+group of neighbouring chunks at a time.
+Writing (`File(path, "w")`, `create_dataset`, `create_group`, `attrs[...] =`)
+covers `float64`, `float32`, `int64`, `int32` and `uint8` arrays. The tests
+in `crates/clawhdf5-py/tests` compare every read with h5py; run them with
+`pip install pytest h5py && pytest crates/clawhdf5-py/tests`.
+
 ### Agent Memory
 
 ```rust

crates/clawhdf5-format/src/chunked_read.rs

@@ -513,6 +513,10 @@ fn collect_chunk_info_inner(
 ///
 /// Chunks are stored contiguously starting at `base_address`. No stored index;
 /// addresses are computed from the chunk position.
+///
+/// `chunk_dimensions` are the spatial chunk dimensions, one per entry of
+/// `dataset_dims` — not the layout message's list, which carries the element
+/// size as an extra last dimension.
 pub fn generate_implicit_chunks(
     base_address: u64,
     dataset_dims: &[u64],

crates/clawhdf5-format/src/data_read.rs

@@ -288,9 +288,11 @@ pub fn read_raw_data_indexed(
 
 /// Read raw bytes for only the selected elements of a dataset.
 ///
-/// For chunked layouts, only chunks that intersect the selection are read
-/// and decompressed. For compact/contiguous layouts, the full data is read
-/// and then the selection is extracted.
+/// When the selection's bounding box covers at most half the dataset, only
+/// that box is materialised — the overlapping rows of a contiguous dataset,
+/// the overlapping chunks of a chunked one, whatever its chunk index (see
+/// [`crate::partial_read`]). Otherwise, and for compact and virtual
+/// layouts, the whole dataset is decoded and the selection extracted.
 #[allow(clippy::too_many_arguments)]
 pub fn read_raw_data_selection(
     file_data: &[u8],
@@ -358,85 +360,17 @@ pub fn read_raw_data_selection(
         }
         DataLayout::Chunked {
             chunk_dimensions,
-            btree_address,
             version,
-            chunk_index_type,
             ..
         } => {
+            // `partial_read` declined (a bounding box covering most of the
+            // dataset, or a selection it doesn't box), so decode every chunk
+            // and pick the selection out, whatever the chunk index. This arm
+            // used to enumerate the chunks first — passing the layout's
+            // chunk dimensions, element-size dimension included, to the
+            // implicit-index generator, which then indexed past the rank and
+            // panicked — only to decode the full dataset anyway.
             crate::chunked_read::chunk_geometry(chunk_dimensions, *version, dataspace, elem_size)?;
-            // For chunked data, only read chunks that intersect the selection
-            let chunk_dims: Vec<u64> = chunk_dimensions.iter().map(|&d| d as u64).collect();
-            let rank = dims.len();
-
-            // Collect chunk info from B-tree
-            let chunks = if *version == 4 {
-                match chunk_index_type {
-                    Some(2) => {
-                        // Implicit index
-                        crate::chunked_read::generate_implicit_chunks(
-                            btree_address.unwrap_or(0),
-                            dims,
-                            chunk_dimensions,
-                            elem_size as u32,
-                        )
-                    }
-                    _ => {
-                        if let Some(_addr) = btree_address {
-                            // Use extensible array or fixed array
-                            // Fall back to full read for complex v4 index types
-                            let full_data = read_raw_data_full(
-                                file_data,
-                                layout,
-                                dataspace,
-                                datatype,
-                                pipeline,
-                                offset_size,
-                                length_size,
-                            )?;
-                            return extract_selection_from_buffer(
-                                &full_data, dims, elem_size, selection,
-                            );
-                        } else {
-                            return Ok(Vec::new());
-                        }
-                    }
-                }
-            } else {
-                // v3: B-tree v1
-                if let Some(addr) = btree_address {
-                    crate::chunked_read::collect_chunk_info_checked(
-                        file_data,
-                        *addr,
-                        chunk_dimensions,
-                        offset_size,
-                        length_size,
-                    )?
-                } else {
-                    return Ok(Vec::new());
-                }
-            };
-
-            // Filter chunks to only those that intersect the selection
-            let intersecting: Vec<_> = chunks
-                .iter()
-                .filter(|ci| {
-                    let offsets: Vec<u64> = ci.offsets.iter().take(rank).copied().collect();
-                    selection.intersects_chunk(&offsets, &chunk_dims[..rank])
-                })
-                .collect();
-
-            if intersecting.is_empty() {
-                return Ok(Vec::new());
-            }
-
-            // Decompress only the intersecting chunks
-            let _chunk_total_bytes: usize =
-                chunk_dims.iter().map(|&d| d as usize).product::<usize>() * elem_size;
-            let _element_size_u32 = elem_size as u32;
-
-            // First, assemble only the intersecting chunks into a partial buffer,
-            // then extract the selection. For simplicity, we assemble into a full
-            // dataset buffer and extract (same as contiguous path).
             let full_data = read_raw_data_full(
                 file_data,
                 layout,

crates/clawhdf5-py/Cargo.toml

@@ -3,7 +3,7 @@ name = "clawhdf5-py"
 version = "2.7.0"
 edition = "2024"
 rust-version.workspace = true
-description = "Python bindings for rustyhdf5 — a pure-Rust HDF5 library"
+description = "Python bindings for clawhdf5 — a pure-Rust HDF5 library"
 license = "MIT"
 repository = "https://git.redclaw.dev/quantumclaw/clawhdf5"
 readme = "README.md"

crates/clawhdf5-py/README.md

@@ -3,23 +3,81 @@
 [![crates.io](https://img.shields.io/crates/v/clawhdf5-py.svg)](https://crates.io/crates/clawhdf5-py)
 [![docs.rs](https://docs.rs/clawhdf5-py/badge.svg)](https://docs.rs/clawhdf5-py)
 
-Python bindings for clawhdf5 — a pure-Rust HDF5 library.
+Python bindings for clawhdf5 — a pure-Rust HDF5 library. The package is
+`clawhdf5` (`import clawhdf5`); it needs numpy and no libhdf5.
 
-## Features
+## Install
 
-- h5py-compatible API (`File`, `Group`, `Dataset`)
-- NumPy array integration
-- Read and write HDF5 files from Python with no C dependencies
+Not on PyPI yet. Build it into a virtualenv with [maturin](https://www.maturin.rs):
 
-## Usage
+```bash
+pip install maturin numpy
+cd crates/clawhdf5-py
+maturin develop --release
+python -c "import clawhdf5; print(clawhdf5.__version__)"
+```
+
+## Reading
+
+The read API follows h5py:
 
 ```python
+import numpy as np
 import clawhdf5
 
-with clawhdf5.File('data.h5', 'r') as f:
-    data = f['/dataset'][:]
+with clawhdf5.File("data.h5", "r") as f:
+    f.keys(), f["group"].items(), "group/data" in f
+    ds = f["group/data"]           # or f["/group/data"], f["group"]["data"]
+    ds.shape, ds.dtype, ds.attrs["units"]
+    ds[10:20, ::2]                 # a small selection reads only its chunks
+    ds[-1], ds[..., 0], ds[[1, 4, 7]]
+    np.asarray(ds)
+    f["table"]["id"]               # a compound field
 ```
 
+- `Dataset.dtype` is the numpy dtype h5py reports: integers and IEEE floats
+  of every width in either byte order, `bool`, enums (with
+  `dtype.metadata['enum']`), complex, `S<n>` fixed strings, `object` for
+  variable-length strings (`bytes` values) and sequences (array values),
+  `V<n>` opaque, array types, and compounds as structured dtypes.
+  Other types raise `TypeError`.
+- Keys are h5py's: integers, slices with a positive step, `...`, one
+  increasing list of integers, compound field names. Each maps onto a
+  hyperslab selection. `None`, negative steps and boolean masks are refused
+  with h5py's errors.
+- What is read from the file: a selection whose bounding box covers at
+  most half the dataset decodes only the chunks (or contiguous rows) the box
+  overlaps. The library decodes the whole dataset for a larger box
+  (including a strided slice such as `ds[::100]` across a chunked dataset),
+  and for compact, virtual and unwritten datasets and chunked ones with a
+  non-default fill value. An index list is read one group of neighbouring
+  chunks at a time (a new group only past a chunk with no selected index),
+  so each chunk is decoded once. `ds[()]`, `ds[...]` and `np.asarray(ds)`
+  use the file's chunk cache; other selections do not.
+- The bytes the library reads become the numpy array's buffer without a
+  copy, and the read runs with the GIL released, so threads read in
+  parallel. A bug in the library (a Rust panic) raises
+  `clawhdf5.InternalError`, a `RuntimeError`.
+- Attributes return what h5py returns; `clawhdf5.Empty` stands for a null
+  dataspace (h5py's `Empty`).
+
+## Writing
+
+`clawhdf5.File(path, "w")` with `create_dataset(name, data=array,
+chunks=..., compression="gzip")`, `create_group` and `attrs[...] = ...`
+writes `float64`, `float32`, `int64`, `int32` and `uint8` arrays; the file is
+written on `close()`.
+
+## Tests
+
+```bash
+pip install pytest h5py
+pytest crates/clawhdf5-py/tests
+```
+
+`tests/test_read_vs_h5py.py` compares every read with h5py on a file h5py
+writes. `scripts/ci-test.sh` builds the wheel and runs these in CI.
+
 ## License
 
 MIT

crates/clawhdf5-py/pyproject.toml

@@ -3,12 +3,15 @@ requires = ["maturin>=1.0,<2.0"]
 build-backend = "maturin"
 
 [project]
-name = "rustyhdf5"
+name = "clawhdf5"
 version = "2.7.0"
-description = "Python bindings for rustyhdf5 — a pure-Rust HDF5 library"
+description = "Python bindings for clawhdf5 — a pure-Rust HDF5 library"
 requires-python = ">=3.8"
 license = { text = "MIT" }
 dependencies = ["numpy"]
 
 [tool.maturin]
 features = ["extension-module"]
+# The extension module is `clawhdf5` (the cdylib's [lib] name): the
+# distribution, the import name and the #[pymodule] all agree.
+module-name = "clawhdf5"

crates/clawhdf5-py/src/attrs.rs

@@ -1,24 +1,31 @@
 //! PyAttrs — dict-like access to HDF5 attributes.
 
-use std::collections::HashMap;
 use std::sync::{Arc, Mutex};
 
+use clawhdf5_format::attribute::AttributeMessage;
+use pyo3::exceptions::{PyKeyError, PyTypeError, PyValueError};
 use pyo3::prelude::*;
-use pyo3::types::PyList;
+use pyo3::types::{PyList, PyTuple};
 
-use crate::{OwnedAttrValue, attr_value_to_py, py_to_attr_value};
+use crate::convert::{Converter, Elements, resolve_vl};
+use crate::{OwnedAttrValue, PyEmpty, attr_value_to_py, node, py_to_attr_value};
 
 /// Backing storage for attributes.
 enum AttrsInner {
-    /// Read-only attributes from an existing HDF5 object.
-    Read(HashMap<String, clawhdf5_rs::AttrValue>),
+    /// Attributes of an object in a file opened for reading, sorted by name.
+    Read {
+        file: Arc<clawhdf5_rs::File>,
+        attrs: Vec<AttributeMessage>,
+    },
     /// Writable attribute list shared with a parent (PyFile or PyGroup).
     Write(Arc<Mutex<Vec<(String, OwnedAttrValue)>>>),
 }
 
 /// Dict-like access to HDF5 attributes.
 ///
-/// In read mode, provides immutable access to attribute key/value pairs.
+/// In read mode, values are what h5py returns: numpy scalars for scalar
+/// attributes, numpy arrays otherwise, `str` for variable-length strings,
+/// `numpy.bytes_` for fixed-length ones, and `Empty` for a null dataspace.
 /// In write mode, attributes set here are accumulated and written when
 /// the parent file is closed.
 #[pyclass(name = "Attrs")]
@@ -27,11 +34,13 @@ pub struct PyAttrs {
 }
 
 impl PyAttrs {
-    /// Create a read-only attrs from an existing attribute map.
-    pub(crate) fn from_read(map: HashMap<String, clawhdf5_rs::AttrValue>) -> Self {
-        Self {
-            inner: AttrsInner::Read(map),
-        }
+    /// The attributes of the object at `addr` (whose path is `path`) in a
+    /// file opened for reading.
+    pub(crate) fn read(file: Arc<clawhdf5_rs::File>, addr: u64, path: &str) -> PyResult<Self> {
+        let attrs = node::attributes(&file, addr, path)?;
+        Ok(Self {
+            inner: AttrsInner::Read { file, attrs },
+        })
     }
 
     /// Create a writable attrs that shares storage with a parent object.
@@ -46,11 +55,11 @@ impl PyAttrs {
 impl PyAttrs {
     fn __getitem__(&self, py: Python<'_>, key: &str) -> PyResult<Py<PyAny>> {
         match &self.inner {
-            AttrsInner::Read(map) => match map.get(key) {
-                Some(val) => Ok(attr_value_to_py(py, val)),
-                None => Err(PyErr::new::<pyo3::exceptions::PyKeyError, _>(
-                    key.to_string(),
-                )),
+            AttrsInner::Read { file, attrs } => match attrs.iter().find(|a| a.name == key) {
+                Some(attr) => Ok(attr_to_py(py, file, attr)?.unbind()),
+                None => Err(PyKeyError::new_err(format!(
+                    "Can't open attribute (can't locate attribute: '{key}')"
+                ))),
             },
             AttrsInner::Write(store) => {
                 let guard = store.lock().unwrap();
@@ -60,16 +69,14 @@ impl PyAttrs {
                         return Ok(attr_value_to_py(py, &attr_val));
                     }
                 }
-                Err(PyErr::new::<pyo3::exceptions::PyKeyError, _>(
-                    key.to_string(),
-                ))
+                Err(PyKeyError::new_err(key.to_string()))
             }
         }
     }
 
     fn __setitem__(&self, key: &str, value: &Bound<'_, PyAny>) -> PyResult<()> {
         match &self.inner {
-            AttrsInner::Read(_) => Err(PyErr::new::<pyo3::exceptions::PyIOError, _>(
+            AttrsInner::Read { .. } => Err(PyErr::new::<pyo3::exceptions::PyIOError, _>(
                 "cannot set attributes on a read-only file",
             )),
             AttrsInner::Write(store) => {
@@ -88,14 +95,14 @@ impl PyAttrs {
 
     fn __len__(&self) -> usize {
         match &self.inner {
-            AttrsInner::Read(map) => map.len(),
+            AttrsInner::Read { attrs, .. } => attrs.len(),
             AttrsInner::Write(store) => store.lock().unwrap().len(),
         }
     }
 
     fn __contains__(&self, key: &str) -> bool {
         match &self.inner {
-            AttrsInner::Read(map) => map.contains_key(key),
+            AttrsInner::Read { attrs, .. } => attrs.iter().any(|a| a.name == key),
             AttrsInner::Write(store) => store.lock().unwrap().iter().any(|(k, _)| k == key),
         }
     }
@@ -111,10 +118,20 @@ impl PyAttrs {
         format!("<HDF5 Attrs ({n} members)>")
     }
 
+    /// The value of `key`, or `default` if there is no such attribute.
+    #[pyo3(signature = (key, default=None))]
+    fn get(&self, py: Python<'_>, key: &str, default: Option<Py<PyAny>>) -> PyResult<Py<PyAny>> {
+        if self.__contains__(key) {
+            self.__getitem__(py, key)
+        } else {
+            Ok(default.unwrap_or_else(|| py.None()))
+        }
+    }
+
     /// Return attribute names as a list.
     fn keys(&self, py: Python<'_>) -> PyResult<Py<PyAny>> {
         let names: Vec<String> = match &self.inner {
-            AttrsInner::Read(map) => map.keys().cloned().collect(),
+            AttrsInner::Read { attrs, .. } => attrs.iter().map(|a| a.name.clone()).collect(),
             AttrsInner::Write(store) => store
                 .lock()
                 .unwrap()
@@ -129,7 +146,10 @@ impl PyAttrs {
     /// Return attribute values as a list.
     fn values(&self, py: Python<'_>) -> PyResult<Py<PyAny>> {
         let vals: Vec<Py<PyAny>> = match &self.inner {
-            AttrsInner::Read(map) => map.values().map(|v| attr_value_to_py(py, v)).collect(),
+            AttrsInner::Read { file, attrs } => attrs
+                .iter()
+                .map(|a| attr_to_py(py, file, a).map(Bound::unbind))
+                .collect::<PyResult<_>>()?,
             AttrsInner::Write(store) => store
                 .lock()
                 .unwrap()
@@ -147,10 +167,10 @@ impl PyAttrs {
     /// Return attribute (key, value) pairs as a list of tuples.
     fn items(&self, py: Python<'_>) -> PyResult<Py<PyAny>> {
         let pairs: Vec<(String, Py<PyAny>)> = match &self.inner {
-            AttrsInner::Read(map) => map
+            AttrsInner::Read { file, attrs } => attrs
                 .iter()
-                .map(|(k, v)| (k.clone(), attr_value_to_py(py, v)))
-                .collect(),
+                .map(|a| Ok((a.name.clone(), attr_to_py(py, file, a)?.unbind())))
+                .collect::<PyResult<_>>()?,
             AttrsInner::Write(store) => store
                 .lock()
                 .unwrap()
@@ -166,28 +186,69 @@ impl PyAttrs {
     }
 }
 
+/// An attribute's value as h5py returns it.
+fn attr_to_py<'py>(
+    py: Python<'py>,
+    file: &clawhdf5_rs::File,
+    attr: &AttributeMessage,
+) -> PyResult<Bound<'py, PyAny>> {
+    crate::no_panic(|| {
+        let sb = file.superblock();
+        let conv = Converter::new(py, &attr.datatype, sb.offset_size)
+            .map_err(|e| prefix_err(py, &attr.name, e))?;
+        if node::is_null(&attr.dataspace) {
+            return Ok(PyEmpty::new(conv.dtype).into_pyobject(py)?.into_any());
+        }
+        let shape: Vec<usize> = attr
+            .dataspace
+            .dimensions
+            .iter()
+            .map(|&d| d as usize)
+            .collect();
+        let n: usize = shape.iter().product();
+        let data = if conv.is_vl() {
+            let want = n * conv.elem_size;
+            if attr.raw_data.len() < want {
+                return Err(PyValueError::new_err(format!(
+                    "attribute {}: {} bytes of variable-length references, expected {want}",
+                    attr.name,
+                    attr.raw_data.len(),
+                )));
+            }
+            let raw = &attr.raw_data[..want];
+            let file_data = file.as_bytes();
+            let (osz, lsz, unit) = (sb.offset_size, sb.length_size, conv.vl_unit);
+            Elements::Vl(
+                py.detach(|| resolve_vl(file_data, raw, n, osz, lsz, unit))
+                    .map_err(|e| PyValueError::new_err(format!("attribute {}: {e}", attr.name)))?,
+            )
+        } else {
+            Elements::Bytes(attr.raw_data.clone())
+        };
+        let arr = conv
+            .to_array(py, data, &shape, true)
+            .map_err(|e| prefix_err(py, &attr.name, e))?;
+        if shape.is_empty() {
+            // A scalar dataspace: h5py returns the element itself.
+            return arr.get_item(PyTuple::empty(py));
+        }
+        Ok(arr)
+    })
+}
+
+fn prefix_err(py: Python<'_>, name: &str, e: PyErr) -> PyErr {
+    let msg = format!("attribute {name}: {}", e.value(py));
+    if e.is_instance_of::<PyTypeError>(py) {
+        PyTypeError::new_err(msg)
+    } else {
+        PyValueError::new_err(msg)
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
 
-    #[test]
-    fn read_attrs_len() {
-        let mut map = HashMap::new();
-        map.insert("a".into(), clawhdf5_rs::AttrValue::I64(1));
-        map.insert("b".into(), clawhdf5_rs::AttrValue::F64(2.0));
-        let attrs = PyAttrs::from_read(map);
-        assert_eq!(attrs.__len__(), 2);
-    }
-
-    #[test]
-    fn read_attrs_contains() {
-        let mut map = HashMap::new();
-        map.insert("x".into(), clawhdf5_rs::AttrValue::String("hello".into()));
-        let attrs = PyAttrs::from_read(map);
-        assert!(attrs.__contains__("x"));
-        assert!(!attrs.__contains__("y"));
-    }
-
     #[test]
     fn write_attrs_len() {
         let store = Arc::new(Mutex::new(Vec::new()));

crates/clawhdf5-py/src/convert.rs

@@ -0,0 +1,600 @@
+//! HDF5 datatypes as numpy dtypes, and element bytes as numpy arrays.
+//!
+//! The dtype a file's datatype maps to is the one h5py reports for it
+//! (byte order kept, compound offsets and padding kept, `r`/`i` compounds as
+//! complex, the `FALSE`/`TRUE` enum as `bool`, fixed strings as `S<n>`,
+//! variable-length data as `object`). For every fixed-size type that dtype
+//! describes the file's element bytes exactly, so the bytes the library
+//! returns become the array's buffer as they are: the `Vec<u8>` is handed to
+//! numpy without a copy and viewed as the dtype.
+//!
+//! Anything this mapping cannot describe exactly — non-IEEE floats, integers
+//! with padding bits, VAX byte order, references, bitfields, time, and
+//! variable-length members inside compounds or arrays — is a `TypeError`,
+//! never a best-effort guess.
+
+use std::collections::HashMap;
+
+use clawhdf5_format::datatype::{CharacterSet, Datatype, DatatypeByteOrder};
+use clawhdf5_format::global_heap::GlobalHeapCollection;
+use numpy::PyArray1;
+use pyo3::exceptions::{PyTypeError, PyValueError};
+use pyo3::prelude::*;
+use pyo3::types::{PyBytes, PyDict, PyList, PyString, PyTuple};
+
+/// How the elements of a datatype become Python values.
+#[derive(Clone, Debug, PartialEq)]
+pub(crate) enum Layout {
+    /// Fixed-size elements numpy reads as they are.
+    Fixed,
+    /// A top-level HDF5 array type: elements are viewed as the base dtype and
+    /// the array's dimensions are appended to the shape (as h5py does).
+    Subarray(Vec<usize>),
+    /// Variable-length string: a global heap reference per element.
+    VlString { utf8: bool },
+    /// Variable-length sequence of a fixed-size base type.
+    VlSequence,
+}
+
+/// Everything needed to turn a dataset's or attribute's bytes into numpy.
+pub(crate) struct Converter {
+    /// The dtype reported to Python (`Dataset.dtype`).
+    pub dtype: Py<PyAny>,
+    /// The dtype the element bytes are viewed as: `dtype` itself, the base
+    /// of a subarray, or the base of a variable-length sequence.
+    pub view: Py<PyAny>,
+    pub layout: Layout,
+    /// Bytes per element in the raw buffer the library returns.
+    pub elem_size: usize,
+    /// For variable-length data, bytes per unit of an element's stored
+    /// length: 1 for strings, the base type's size for sequences.
+    pub vl_unit: usize,
+}
+
+fn unsupported(what: impl std::fmt::Display) -> PyErr {
+    PyTypeError::new_err(format!(
+        "clawhdf5 cannot read this datatype into numpy: {what}"
+    ))
+}
+
+fn byte_order_char(order: &DatatypeByteOrder, size: u32) -> PyResult<&'static str> {
+    if size == 1 {
+        return Ok("|");
+    }
+    match order {
+        DatatypeByteOrder::LittleEndian => Ok("<"),
+        DatatypeByteOrder::BigEndian => Ok(">"),
+        DatatypeByteOrder::Vax => Err(unsupported("VAX byte order")),
+    }
+}
+
+/// The numpy format string of an integer type, if it is a plain one.
+fn int_format(dt: &Datatype) -> PyResult<String> {
+    match dt {
+        Datatype::FixedPoint {
+            size,
+            byte_order,
+            signed,
+            bit_offset,
+            bit_precision,
+        } => {
+            if !matches!(size, 1 | 2 | 4 | 8) {
+                return Err(unsupported(format!("{size}-byte integer")));
+            }
+            if *bit_offset != 0 || u32::from(*bit_precision) != size * 8 {
+                return Err(unsupported(format!(
+                    "integer with {bit_precision} significant bits at offset {bit_offset} in {size} bytes"
+                )));
+            }
+            let kind = if *signed { 'i' } else { 'u' };
+            Ok(format!(
+                "{}{kind}{size}",
+                byte_order_char(byte_order, *size)?
+            ))
+        }
+        other => Err(unsupported(format!("{other:?} is not an integer"))),
+    }
+}
+
+/// The numpy format string of an IEEE 754 binary16/32/64 type.
+fn float_format(dt: &Datatype) -> PyResult<String> {
+    let Datatype::FloatingPoint {
+        size,
+        byte_order,
+        bit_offset,
+        bit_precision,
+        exponent_location,
+        exponent_size,
+        mantissa_location,
+        mantissa_size,
+        exponent_bias,
+    } = dt
+    else {
+        return Err(unsupported(format!("{dt:?} is not a float")));
+    };
+    // (exponent location, exponent size, mantissa size, bias) of IEEE 754.
+    let ieee = match size {
+        2 => (10, 5, 10, 15),
+        4 => (23, 8, 23, 127),
+        8 => (52, 11, 52, 1023),
+        _ => return Err(unsupported(format!("{size}-byte float"))),
+    };
+    let layout = (
+        *exponent_location,
+        *exponent_size,
+        *mantissa_size,
+        *exponent_bias,
+    );
+    if *bit_offset != 0
+        || u32::from(*bit_precision) != size * 8
+        || *mantissa_location != 0
+        || layout != ieee
+    {
+        return Err(unsupported(format!(
+            "non-IEEE {size}-byte float (exponent {exponent_size} bits at {exponent_location}, \
+             mantissa {mantissa_size} bits at {mantissa_location}, bias {exponent_bias})"
+        )));
+    }
+    Ok(format!("{}f{size}", byte_order_char(byte_order, *size)?))
+}
+
+/// `r`/`i` compounds of two identical IEEE floats are complex numbers in h5py.
+fn complex_format(
+    size: u32,
+    members: &[clawhdf5_format::datatype::CompoundMember],
+) -> Option<String> {
+    let [re, im] = members else { return None };
+    if re.name != "r" || im.name != "i" || re.datatype != im.datatype {
+        return None;
+    }
+    let Datatype::FloatingPoint {
+        size: fsize,
+        byte_order,
+        ..
+    } = &re.datatype
+    else {
+        return None;
+    };
+    if !matches!(fsize, 4 | 8)
+        || re.byte_offset != 0
+        || im.byte_offset != u64::from(*fsize)
+        || size != 2 * fsize
+    {
+        return None;
+    }
+    float_format(&re.datatype).ok()?;
+    let order = byte_order_char(byte_order, *fsize).ok()?;
+    Some(format!("{order}c{}", 2 * fsize))
+}
+
+/// The members of an enum as `{name: value}`.
+fn enum_members<'py>(
+    py: Python<'py>,
+    base: &Datatype,
+    members: &[clawhdf5_format::datatype::EnumMember],
+) -> PyResult<Bound<'py, PyDict>> {
+    let signed = matches!(base, Datatype::FixedPoint { signed: true, .. });
+    let dict = PyDict::new(py);
+    for m in members {
+        let value: Py<PyAny> = if signed {
+            let v = clawhdf5_format::data_read::read_as_i64(&m.value, base)
+                .map_err(|e| PyValueError::new_err(format!("enum member {}: {e}", m.name)))?;
+            let v = *v.first().ok_or_else(|| {
+                PyValueError::new_err(format!("enum member {} has no value", m.name))
+            })?;
+            v.into_pyobject(py)?.into_any().unbind()
+        } else {
+            let v = clawhdf5_format::data_read::read_as_u64(&m.value, base)
+                .map_err(|e| PyValueError::new_err(format!("enum member {}: {e}", m.name)))?;
+            let v = *v.first().ok_or_else(|| {
+                PyValueError::new_err(format!("enum member {} has no value", m.name))
+            })?;
+            v.into_pyobject(py)?.into_any().unbind()
+        };
+        dict.set_item(&m.name, value)?;
+    }
+    Ok(dict)
+}
+
+/// Whether an enum is h5py's boolean: a one-byte integer with exactly the
+/// members `FALSE` = 0 and `TRUE` = 1.
+fn is_h5py_bool(base: &Datatype, members: &[clawhdf5_format::datatype::EnumMember]) -> bool {
+    if base.type_size() != 1 || members.len() != 2 {
+        return false;
+    }
+    let value = |name: &str| {
+        members
+            .iter()
+            .find(|m| m.name == name)
+            .and_then(|m| m.value.first().copied())
+    };
+    value("FALSE") == Some(0) && value("TRUE") == Some(1)
+}
+
+fn np_dtype<'py>(py: Python<'py>, spec: impl IntoPyObject<'py>) -> PyResult<Bound<'py, PyAny>> {
+    py.import("numpy")?.getattr("dtype")?.call1((spec,))
+}
+
+fn np_dtype_with_metadata<'py>(
+    py: Python<'py>,
+    spec: impl IntoPyObject<'py>,
+    metadata: Bound<'py, PyDict>,
+) -> PyResult<Bound<'py, PyAny>> {
+    let kwargs = PyDict::new(py);
+    kwargs.set_item("metadata", metadata)?;
+    py.import("numpy")?
+        .getattr("dtype")?
+        .call((spec,), Some(&kwargs))
+}
+
+/// The numpy dtype of a fixed-size datatype, whose element bytes numpy can
+/// read as they are.
+pub(crate) fn fixed_dtype<'py>(py: Python<'py>, dt: &Datatype) -> PyResult<Bound<'py, PyAny>> {
+    match dt {
+        Datatype::FixedPoint { .. } => np_dtype(py, int_format(dt)?),
+        Datatype::FloatingPoint { .. } => np_dtype(py, float_format(dt)?),
+        Datatype::String { size, charset, .. } => {
+            if *size == 0 {
+                return Err(unsupported("zero-length fixed string"));
+            }
+            let meta = PyDict::new(py);
+            let enc = match charset {
+                CharacterSet::Ascii => "ascii",
+                CharacterSet::Utf8 => "utf-8",
+            };
+            meta.set_item("h5py_encoding", enc)?;
+            np_dtype_with_metadata(py, format!("S{size}"), meta)
+        }
+        Datatype::Opaque { size, .. } => {
+            if *size == 0 {
+                return Err(unsupported("zero-length opaque type"));
+            }
+            np_dtype(py, format!("V{size}"))
+        }
+        Datatype::Enumeration {
+            base_type, members, ..
+        } => {
+            let base = int_format(base_type)?;
+            if is_h5py_bool(base_type, members) {
+                return np_dtype(py, "?");
+            }
+            let meta = PyDict::new(py);
+            meta.set_item("enum", enum_members(py, base_type, members)?)?;
+            np_dtype_with_metadata(py, base, meta)
+        }
+        Datatype::Compound { size, members } => {
+            if let Some(c) = complex_format(*size, members) {
+                return np_dtype(py, c);
+            }
+            let names = PyList::empty(py);
+            let formats = PyList::empty(py);
+            let offsets = PyList::empty(py);
+            for m in members {
+                let end = m.byte_offset.checked_add(u64::from(m.datatype.type_size()));
+                if end.is_none_or(|end| end > u64::from(*size)) {
+                    return Err(PyValueError::new_err(format!(
+                        "compound member {} lies outside the {size}-byte compound",
+                        m.name
+                    )));
+                }
+                names.append(&m.name)?;
+                formats.append(fixed_dtype(py, &m.datatype).map_err(|e| {
+                    unsupported(format!("compound member {}: {}", m.name, e.value(py)))
+                })?)?;
+                offsets.append(m.byte_offset)?;
+            }
+            let spec = PyDict::new(py);
+            spec.set_item("names", names)?;
+            spec.set_item("formats", formats)?;
+            spec.set_item("offsets", offsets)?;
+            spec.set_item("itemsize", size)?;
+            np_dtype(py, spec)
+        }
+        Datatype::Array {
+            base_type,
+            dimensions,
+        } => {
+            let base = fixed_dtype(py, base_type)?;
+            let dims = PyTuple::new(py, dimensions)?;
+            np_dtype(py, (base, dims))
+        }
+        Datatype::VariableLength { is_string, .. } => Err(unsupported(if *is_string {
+            "variable-length string inside a compound or array type"
+        } else {
+            "variable-length sequence inside a compound or array type"
+        })),
+        Datatype::Reference { .. } => Err(unsupported("object/region references")),
+        Datatype::BitField { .. } => Err(unsupported("bitfield")),
+        Datatype::Time { .. } => Err(unsupported("time")),
+    }
+}
+
+impl Converter {
+    /// The converter for a dataset's or attribute's datatype.
+    pub(crate) fn new(py: Python<'_>, dt: &Datatype, offset_size: u8) -> PyResult<Self> {
+        match dt {
+            Datatype::VariableLength {
+                is_string: true,
+                charset,
+                ..
+            } => {
+                let utf8 = matches!(charset, Some(CharacterSet::Utf8));
+                let meta = PyDict::new(py);
+                if utf8 {
+                    meta.set_item("vlen", py.get_type::<PyString>())?;
+                } else {
+                    meta.set_item("vlen", py.get_type::<PyBytes>())?;
+                }
+                let dtype = np_dtype_with_metadata(py, "O", meta)?;
+                Ok(Self {
+                    view: dtype.clone().unbind(),
+                    dtype: dtype.unbind(),
+                    layout: Layout::VlString { utf8 },
+                    elem_size: vl_ref_size(offset_size)?,
+                    vl_unit: 1,
+                })
+            }
+            Datatype::VariableLength {
+                is_string: false,
+                base_type,
+                ..
+            } => {
+                let base = fixed_dtype(py, base_type)?;
+                let meta = PyDict::new(py);
+                meta.set_item("vlen", &base)?;
+                let dtype = np_dtype_with_metadata(py, "O", meta)?;
+                Ok(Self {
+                    dtype: dtype.unbind(),
+                    view: base.unbind(),
+                    layout: Layout::VlSequence,
+                    elem_size: vl_ref_size(offset_size)?,
+                    vl_unit: base_type.type_size() as usize,
+                })
+            }
+            Datatype::Array {
+                base_type,
+                dimensions,
+            } => {
+                let dtype = fixed_dtype(py, dt)?;
+                let base = fixed_dtype(py, base_type)?;
+                Ok(Self {
+                    dtype: dtype.unbind(),
+                    view: base.unbind(),
+                    layout: Layout::Subarray(dimensions.iter().map(|&d| d as usize).collect()),
+                    elem_size: dt.type_size() as usize,
+                    vl_unit: 0,
+                })
+            }
+            _ => {
+                let dtype = fixed_dtype(py, dt)?;
+                Ok(Self {
+                    view: dtype.clone().unbind(),
+                    dtype: dtype.unbind(),
+                    layout: Layout::Fixed,
+                    elem_size: dt.type_size() as usize,
+                    vl_unit: 0,
+                })
+            }
+        }
+    }
+
+    pub(crate) fn is_vl(&self) -> bool {
+        matches!(self.layout, Layout::VlString { .. } | Layout::VlSequence)
+    }
+
+    /// An empty array of `shape` (some dimension is zero).
+    pub(crate) fn empty<'py>(
+        &self,
+        py: Python<'py>,
+        shape: &[usize],
+    ) -> PyResult<Bound<'py, PyAny>> {
+        let np = py.import("numpy")?;
+        match &self.layout {
+            Layout::Subarray(dims) => {
+                let mut full = shape.to_vec();
+                full.extend_from_slice(dims);
+                np.call_method1("empty", (PyTuple::new(py, full)?, self.view.bind(py)))
+            }
+            _ => np.call_method1("empty", (PyTuple::new(py, shape)?, self.dtype.bind(py))),
+        }
+    }
+
+    /// Turn decoded element data into a numpy array of `shape`.
+    ///
+    /// `str_values` decodes variable-length strings to `str` (what h5py
+    /// does for attributes) instead of `bytes` (what it does for datasets).
+    pub(crate) fn to_array<'py>(
+        &self,
+        py: Python<'py>,
+        data: Elements,
+        shape: &[usize],
+        str_values: bool,
+    ) -> PyResult<Bound<'py, PyAny>> {
+        let n: usize = shape.iter().product();
+        match (data, &self.layout) {
+            (Elements::Bytes(bytes), Layout::Fixed) => {
+                bytes_as_array(py, bytes, self.view.bind(py), shape)
+            }
+            (Elements::Bytes(bytes), Layout::Subarray(dims)) => {
+                let mut full = shape.to_vec();
+                full.extend_from_slice(dims);
+                bytes_as_array(py, bytes, self.view.bind(py), &full)
+            }
+            (Elements::Vl(items), Layout::VlString { .. }) => {
+                check_count(items.len(), n)?;
+                let mut objs: Vec<Py<PyAny>> = Vec::with_capacity(items.len());
+                for item in items {
+                    let obj = if str_values {
+                        PyBytes::new(py, &item)
+                            .call_method1("decode", ("utf-8", "surrogateescape"))?
+                            .unbind()
+                    } else {
+                        PyBytes::new(py, &item).into_any().unbind()
+                    };
+                    objs.push(obj);
+                }
+                object_array(py, objs, shape)
+            }
+            (Elements::Vl(items), Layout::VlSequence) => {
+                check_count(items.len(), n)?;
+                let base = self.view.bind(py);
+                let itemsize: usize = base.getattr("itemsize")?.extract()?;
+                let mut objs: Vec<Py<PyAny>> = Vec::with_capacity(items.len());
+                for item in items {
+                    if item.len() % itemsize != 0 {
+                        return Err(PyValueError::new_err(format!(
+                            "variable-length element of {} bytes is not a whole number of {itemsize}-byte values",
+                            item.len()
+                        )));
+                    }
+                    let len = item.len() / itemsize;
+                    objs.push(bytes_as_array(py, item, base, &[len])?.unbind());
+                }
+                object_array(py, objs, shape)
+            }
+            _ => Err(PyValueError::new_err(
+                "internal error: element data does not match the datatype",
+            )),
+        }
+    }
+}
+
+/// Element data as read, before it becomes numpy.
+pub(crate) enum Elements {
+    /// The elements' bytes, back to back.
+    Bytes(Vec<u8>),
+    /// Each variable-length element's bytes, resolved from the global heap.
+    Vl(Vec<Vec<u8>>),
+}
+
+fn vl_ref_size(offset_size: u8) -> PyResult<usize> {
+    // The library sizes a variable-length element as 16 bytes (a length, an
+    // 8-byte heap address and an index) whatever the file's offset size.
+    // Refuse the other sizes rather than read misaligned references.
+    if offset_size != 8 {
+        return Err(unsupported(format!(
+            "variable-length data in a file with {offset_size}-byte offsets"
+        )));
+    }
+    Ok(4 + usize::from(offset_size) + 4)
+}
+
+fn check_count(got: usize, want: usize) -> PyResult<()> {
+    if got != want {
+        return Err(PyValueError::new_err(format!(
+            "read {got} elements, expected {want}"
+        )));
+    }
+    Ok(())
+}
+
+/// A numpy array over `bytes` without copying them: the `Vec` becomes the
+/// array's buffer and is viewed as `dtype` with `shape`.
+pub(crate) fn bytes_as_array<'py>(
+    py: Python<'py>,
+    bytes: Vec<u8>,
+    dtype: &Bound<'py, PyAny>,
+    shape: &[usize],
+) -> PyResult<Bound<'py, PyAny>> {
+    let itemsize: usize = dtype.getattr("itemsize")?.extract()?;
+    let n: usize = shape.iter().product();
+    if n.checked_mul(itemsize) != Some(bytes.len()) {
+        return Err(PyValueError::new_err(format!(
+            "read {} bytes, expected {n} elements of {itemsize} bytes",
+            bytes.len()
+        )));
+    }
+    let shape = PyTuple::new(py, shape)?;
+    if n == 0 {
+        return py.import("numpy")?.call_method1("empty", (shape, dtype));
+    }
+    let raw = PyArray1::from_vec(py, bytes);
+    let arr = raw
+        .call_method1("view", (dtype,))?
+        .call_method1("reshape", (shape,))?;
+    // A `Vec<u8>` carries no alignment promise. numpy copes with unaligned
+    // arrays, but slowly and not in every routine, so hand out an aligned
+    // copy in the (allocator-dependent, rare) case the buffer is not.
+    if !arr
+        .getattr("flags")?
+        .getattr("aligned")?
+        .extract::<bool>()?
+    {
+        return arr.call_method0("copy");
+    }
+    Ok(arr)
+}
+
+fn object_array<'py>(
+    py: Python<'py>,
+    objs: Vec<Py<PyAny>>,
+    shape: &[usize],
+) -> PyResult<Bound<'py, PyAny>> {
+    let arr = PyArray1::from_vec(py, objs);
+    arr.call_method1("reshape", (PyTuple::new(py, shape)?,))
+}
+
+/// Resolve variable-length elements (global heap references in `raw`) to
+/// their bytes: each element's stored length times `unit` (1 for strings,
+/// the base type's size for sequences). Pure Rust, so it runs without the
+/// GIL.
+pub(crate) fn resolve_vl(
+    file_data: &[u8],
+    raw: &[u8],
+    count: usize,
+    offset_size: u8,
+    length_size: u8,
+    unit: usize,
+) -> Result<Vec<Vec<u8>>, String> {
+    let refs = clawhdf5_format::vl_data::parse_vl_references(raw, count as u64, offset_size)
+        .map_err(|e| e.to_string())?;
+    let undefined = match offset_size {
+        2 => 0xFFFF,
+        4 => 0xFFFF_FFFF,
+        _ => u64::MAX,
+    };
+    let mut collections: HashMap<u64, GlobalHeapCollection> = HashMap::new();
+    let mut out = Vec::with_capacity(refs.len());
+    for vl in &refs {
+        if vl.collection_address == 0 || vl.collection_address == undefined {
+            if vl.length != 0 {
+                return Err(format!(
+                    "variable-length element of length {} has no heap address",
+                    vl.length
+                ));
+            }
+            out.push(Vec::new());
+            continue;
+        }
+        let coll = match collections.entry(vl.collection_address) {
+            std::collections::hash_map::Entry::Occupied(e) => e.into_mut(),
+            std::collections::hash_map::Entry::Vacant(e) => {
+                let addr = usize::try_from(vl.collection_address)
+                    .map_err(|_| "global heap address out of range".to_string())?;
+                e.insert(
+                    GlobalHeapCollection::parse(file_data, addr, length_size)
+                        .map_err(|e| e.to_string())?,
+                )
+            }
+        };
+        let index = u16::try_from(vl.object_index)
+            .map_err(|_| format!("global heap object index {} out of range", vl.object_index))?;
+        let obj = coll.get_object(index).ok_or_else(|| {
+            format!(
+                "global heap object {index} not found in the collection at {}",
+                vl.collection_address
+            )
+        })?;
+        let need = (vl.length as usize)
+            .checked_mul(unit)
+            .ok_or("variable-length element too long")?;
+        if need > obj.data.len() {
+            return Err(format!(
+                "variable-length element of {need} bytes in a {}-byte heap object",
+                obj.data.len()
+            ));
+        }
+        out.push(obj.data[..need].to_vec());
+    }
+    Ok(out)
+}

crates/clawhdf5-py/src/dataset.rs

@@ -1,241 +1,373 @@
-//! PyDataset — read access to HDF5 datasets with numpy integration.
+//! PyDataset — h5py-style read access to HDF5 datasets.
+//!
+//! `ds[key]` parses the key into hyperslab selections (see `select`) and
+//! reads them through the facade's `read_selection`, which decodes only the
+//! chunks a small selection touches (see its docs for when it decodes the
+//! whole dataset instead); the
+//! bytes it returns become the numpy array's buffer without a copy (see
+//! `convert`). All file access and decoding runs with the GIL released, so
+//! Python threads reading the same or different datasets run in parallel.
 
 use std::sync::Arc;
 
-use numpy::PyArrayDyn;
-use numpy::ndarray::{ArrayD, IxDyn};
+use clawhdf5_format::datatype::Datatype;
+use clawhdf5_format::object_header::ObjectHeader;
+use pyo3::exceptions::{PyTypeError, PyValueError};
 use pyo3::prelude::*;
-use pyo3::types::PyList;
-
-use clawhdf5_rs::DType;
+use pyo3::types::{PyList, PyTuple};
 
 use crate::attrs::PyAttrs;
-use crate::to_py_err;
+use crate::convert::{Converter, Elements, resolve_vl};
+use crate::select::{self, Plan};
+use crate::{PyEmpty, node, to_py_err};
 
-/// A handle to an HDF5 dataset (read mode).
+/// A dataset in a file opened for reading.
 ///
-/// Supports numpy-style indexing via `__getitem__`:
 /// ```python
-/// ds = f['dataset_name']
-/// data = ds[:]           # read all data as numpy array
-/// shape = ds.shape
-/// dtype = ds.dtype
+/// ds = f['group/dataset']
+/// ds.shape, ds.dtype, ds.attrs['units']
+/// block = ds[10:20, ::2]      # a small selection reads only its chunks
 /// ```
 #[pyclass(name = "Dataset")]
 pub struct PyDataset {
     file: Arc<clawhdf5_rs::File>,
     path: String,
-    cached_shape: Vec<u64>,
-    cached_dtype: DType,
+    /// Where the dataset's object header is: reads open it from here rather
+    /// than resolve `path` again.
+    addr: u64,
+    /// `None` for a dataset with a null dataspace (h5py's `Empty`).
+    shape: Option<Vec<u64>>,
+    /// The chunk shape, for a chunked dataset.
+    chunks: Option<Vec<u64>>,
+    datatype: Datatype,
+    /// Why the datatype cannot be read into numpy, if it cannot.
+    conv: Result<Converter, String>,
 }
 
 impl PyDataset {
-    pub fn new(file: Arc<clawhdf5_rs::File>, path: String) -> PyResult<Self> {
-        let ds = file.dataset(&path).map_err(to_py_err)?;
-        let cached_shape = ds.shape().map_err(to_py_err)?;
-        let cached_dtype = ds.dtype().map_err(to_py_err)?;
+    pub(crate) fn open(
+        py: Python<'_>,
+        file: Arc<clawhdf5_rs::File>,
+        path: String,
+        addr: u64,
+        hdr: &ObjectHeader,
+    ) -> PyResult<Self> {
+        crate::no_panic(|| {
+            let null = node::is_null(&node::dataspace(&file, hdr)?);
+            let (shape, datatype) = {
+                let ds = file.dataset_at(addr).map_err(to_py_err)?;
+                let shape = if null {
+                    None
+                } else {
+                    Some(ds.shape().map_err(to_py_err)?)
+                };
+                (shape, ds.raw_datatype().map_err(to_py_err)?)
+            };
+            let conv = Converter::new(py, &datatype, file.superblock().offset_size)
+                .map_err(|e| e.value(py).to_string());
+            let chunks = shape
+                .as_ref()
+                .and_then(|s| node::chunk_shape(&file, hdr, s.len()));
             Ok(Self {
                 file,
                 path,
-            cached_shape,
-            cached_dtype,
+                addr,
+                shape,
+                chunks,
+                datatype,
+                conv,
             })
+        })
+    }
+
+    fn converter(&self) -> PyResult<&Converter> {
+        self.conv
+            .as_ref()
+            .map_err(|msg| PyTypeError::new_err(format!("{}: {msg}", node::name(&self.path))))
+    }
+
+    /// Read the selection described by `plan` into a numpy array.
+    fn read_plan<'py>(&self, py: Python<'py>, plan: &Plan) -> PyResult<Bound<'py, PyAny>> {
+        let conv = self.converter()?;
+        let dims = self.shape.as_deref().unwrap_or(&[]);
+        let out_shape = plan.out_shape();
+
+        let arr = if plan.is_empty() {
+            conv.empty(py, &out_shape)?
+        } else {
+            let (vl, elem_size, unit) = (conv.is_vl(), conv.elem_size, conv.vl_unit);
+            let list_axis = plan.list_axis();
+            let chunk_len = match (&self.chunks, list_axis) {
+                (Some(c), Some(a)) => c.get(a).copied(),
+                _ => None,
+            };
+            let (reads, list_axis) = plan.reads(dims, chunk_len, elem_size);
+            let read_shape = plan.read_shape();
+            let file = &*self.file;
+            let addr = self.addr;
+            // Everything below touches only Rust data: release the GIL.
+            let read = || -> Result<Elements, ReadError> {
+                let ds = file.dataset_at(addr)?;
+                let mut blocks = Vec::with_capacity(reads.len());
+                for read in reads {
+                    let raw = ds.read_selection(&read.sel)?;
+                    let mut shape = read.shape;
+                    let want = shape.iter().product::<usize>() * elem_size;
+                    if raw.len() != want {
+                        return Err(ReadError::Other(format!(
+                            "read {} bytes, expected {want}",
+                            raw.len()
+                        )));
+                    }
+                    let raw = match (&read.pick, list_axis) {
+                        (Some(pick), Some(axis)) => {
+                            let kept = select::gather_along(&raw, &shape, axis, pick, elem_size);
+                            shape[axis] = pick.len();
+                            kept
+                        }
+                        _ => raw,
+                    };
+                    blocks.push((raw, shape));
+                }
+                // Several blocks only for a list index: join their bytes
+                // (every byte of every element, padding included) along
+                // that axis before anything becomes numpy.
+                let raw = match (blocks.len(), list_axis) {
+                    (1, _) => blocks.pop().expect("one block").0,
+                    (_, Some(axis)) => select::join_along(&blocks, axis, elem_size),
+                    _ => {
+                        return Err(ReadError::Other(
+                            "several reads without an index list".into(),
+                        ));
+                    }
+                };
+                if !vl {
+                    return Ok(Elements::Bytes(raw));
+                }
+                let sb = file.superblock();
+                let n = read_shape.iter().product();
+                resolve_vl(
+                    file.as_bytes(),
+                    &raw,
+                    n,
+                    sb.offset_size,
+                    sb.length_size,
+                    unit,
+                )
+                .map(Elements::Vl)
+                .map_err(ReadError::Other)
+            };
+            let data = py
+                .detach(|| {
+                    std::panic::catch_unwind(std::panic::AssertUnwindSafe(read))
+                        .unwrap_or_else(|p| Err(ReadError::Panic(crate::panic_text(&*p))))
+                })
+                .map_err(|e| e.into_py(&self.path))?;
+            let joined = conv.to_array(py, data, &read_shape, false)?;
+            // Drop the axes indexed by an integer (length 1 in the blocks).
+            let mut shape = out_shape.clone();
+            if let crate::convert::Layout::Subarray(sub) = &conv.layout {
+                shape.extend_from_slice(sub);
+            }
+            joined.call_method1("reshape", (PyTuple::new(py, shape)?,))?
+        };
+
+        let arr = select_fields(py, arr, &plan.fields)?;
+        if plan.scalar {
+            return arr.get_item(PyTuple::empty(py));
+        }
+        Ok(arr)
     }
 }
 
-/// Map a `DType` to a numpy dtype string.
-fn dtype_to_numpy_str(dt: &DType) -> &'static str {
-    match dt {
-        DType::F64 => "float64",
-        DType::F32 => "float32",
-        DType::I64 => "int64",
-        DType::I32 => "int32",
-        DType::I16 => "int16",
-        DType::I8 => "int8",
-        DType::U64 => "uint64",
-        DType::U32 => "uint32",
-        DType::U16 => "uint16",
-        DType::U8 => "uint8",
-        DType::String | DType::VariableLengthString => "object",
-        _ => "object",
+/// An error from the read closure, turned into a Python error with the GIL.
+enum ReadError {
+    Lib(clawhdf5_rs::Error),
+    Other(String),
+    Panic(String),
 }
+
+impl From<clawhdf5_rs::Error> for ReadError {
+    fn from(e: clawhdf5_rs::Error) -> Self {
+        ReadError::Lib(e)
+    }
+}
+
+impl ReadError {
+    fn into_py(self, path: &str) -> PyErr {
+        match self {
+            ReadError::Lib(e) => to_py_err(e),
+            ReadError::Other(msg) => PyValueError::new_err(format!("{}: {msg}", node::name(path))),
+            ReadError::Panic(msg) => crate::InternalError::new_err(format!(
+                "{}: clawhdf5 internal error (please report it): {msg}",
+                node::name(path)
+            )),
+        }
+    }
+}
+
+/// Keep only the named compound fields, as h5py's `ds['x']` / `ds['x', 'y']`.
+fn select_fields<'py>(
+    py: Python<'py>,
+    arr: Bound<'py, PyAny>,
+    fields: &[String],
+) -> PyResult<Bound<'py, PyAny>> {
+    if fields.is_empty() {
+        return Ok(arr);
+    }
+    let names = arr.getattr("dtype")?.getattr("names")?;
+    if names.is_none() {
+        return Err(PyValueError::new_err(
+            "Field names only allowed for compound types",
+        ));
+    }
+    let names: Vec<String> = names.extract()?;
+    for f in fields {
+        if !names.contains(f) {
+            return Err(PyValueError::new_err(format!(
+                "Field {f} does not appear in this type."
+            )));
+        }
+    }
+    let np = py.import("numpy")?;
+    if let [one] = fields {
+        return np.call_method1("ascontiguousarray", (arr.get_item(one)?,));
+    }
+    let picked = arr.get_item(PyList::new(py, fields)?)?;
+    py.import("numpy.lib.recfunctions")?
+        .call_method1("repack_fields", (picked,))
 }
 
 #[pymethods]
 impl PyDataset {
-    /// The shape of the dataset as a tuple.
+    /// The shape of the dataset (`None` for an empty/null dataspace).
     #[getter]
-    fn shape(&self, py: Python<'_>) -> PyResult<Py<PyAny>> {
-        let tuple = pyo3::types::PyTuple::new(py, self.cached_shape.iter().map(|&d| d as usize))?;
-        Ok(tuple.into_any().unbind())
+    fn shape<'py>(&self, py: Python<'py>) -> PyResult<Bound<'py, PyAny>> {
+        match &self.shape {
+            Some(s) => Ok(PyTuple::new(py, s)?.into_any()),
+            None => Ok(py.None().into_bound(py)),
+        }
     }
 
-    /// The numpy dtype string of the dataset.
+    /// The maximum shape (`None` per unlimited dimension), like h5py.
     #[getter]
-    fn dtype(&self) -> &'static str {
-        dtype_to_numpy_str(&self.cached_dtype)
+    fn maxshape<'py>(&self, py: Python<'py>) -> PyResult<Bound<'py, PyAny>> {
+        crate::no_panic(|| {
+            let Some(shape) = &self.shape else {
+                return Ok(py.None().into_bound(py));
+            };
+            let max = self
+                .file
+                .dataset_at(self.addr)
+                .and_then(|ds| ds.max_dimensions())
+                .map_err(to_py_err)?
+                .unwrap_or_else(|| shape.clone());
+            let items: Vec<Option<u64>> = max
+                .into_iter()
+                .map(|d| (d != u64::MAX).then_some(d))
+                .collect();
+            Ok(PyTuple::new(py, items)?.into_any())
+        })
     }
 
-    /// Attribute access (read-only).
+    /// The dataset's numpy dtype, as h5py reports it.
+    #[getter]
+    fn dtype<'py>(&self, py: Python<'py>) -> PyResult<Bound<'py, PyAny>> {
+        Ok(self.converter()?.dtype.bind(py).clone())
+    }
+
+    #[getter]
+    fn ndim(&self) -> usize {
+        self.shape.as_ref().map_or(0, Vec::len)
+    }
+
+    /// Number of elements (`None` for an empty/null dataspace, as h5py).
+    #[getter]
+    fn size(&self) -> Option<u64> {
+        self.shape.as_ref().map(|s| s.iter().product())
+    }
+
+    /// The dataset's full name, e.g. `/group/data`.
+    #[getter]
+    fn name(&self) -> String {
+        node::name(&self.path)
+    }
+
+    /// The dataset's attributes (read-only, dict-like).
     #[getter]
     fn attrs(&self) -> PyResult<PyAttrs> {
-        let ds = self.file.dataset(&self.path).map_err(to_py_err)?;
-        let map = ds.attrs().map_err(to_py_err)?;
-        Ok(PyAttrs::from_read(map))
+        PyAttrs::read(Arc::clone(&self.file), self.addr, &self.path)
     }
 
-    /// Read data via indexing. Supports `ds[:]`, `ds[0]`, `ds[0:5]`, etc.
-    ///
-    /// The full dataset is always read from the underlying file; the index
-    /// is then applied on the resulting numpy array.
-    fn __getitem__<'py>(&self, py: Python<'py>, key: &Bound<'py, PyAny>) -> PyResult<Py<PyAny>> {
-        let arr = self.read_as_numpy(py)?;
-        let indexed = arr.get_item(key)?;
-        Ok(indexed.unbind())
+    /// Read with h5py indexing: integers, slices with positive steps,
+    /// `...`, one increasing list of integers, and compound field names.
+    /// A selection whose bounding box covers at most half the dataset reads
+    /// only the chunks (or contiguous rows) it overlaps.
+    fn __getitem__<'py>(
+        &self,
+        py: Python<'py>,
+        key: &Bound<'py, PyAny>,
+    ) -> PyResult<Bound<'py, PyAny>> {
+        let Some(dims) = &self.shape else {
+            let is_empty_tuple = key.cast::<PyTuple>().is_ok_and(|t| t.is_empty());
+            let is_ellipsis = key.is_instance_of::<pyo3::types::PyEllipsis>();
+            if is_empty_tuple || is_ellipsis {
+                let empty = PyEmpty::new(self.converter()?.dtype.clone_ref(py));
+                return Ok(empty.into_pyobject(py)?.into_any());
+            }
+            return Err(PyValueError::new_err("Empty datasets cannot be sliced"));
+        };
+        let plan = select::parse(key, dims)?;
+        self.read_plan(py, &plan)
     }
 
-    fn __repr__(&self) -> String {
+    /// `numpy.asarray(ds)` reads the whole dataset.
+    #[pyo3(signature = (dtype=None, copy=None))]
+    fn __array__<'py>(
+        &self,
+        py: Python<'py>,
+        dtype: Option<&Bound<'py, PyAny>>,
+        copy: Option<bool>,
+    ) -> PyResult<Bound<'py, PyAny>> {
+        let _ = copy; // every read is a fresh array
+        let Some(dims) = &self.shape else {
+            return Err(PyValueError::new_err("an empty dataset has no array value"));
+        };
+        let ellipsis = pyo3::types::PyEllipsis::get(py).to_owned().into_any();
+        let plan = select::parse(&ellipsis, dims)?;
+        let arr = self.read_plan(py, &plan)?;
+        match dtype {
+            Some(dt) => arr.call_method1("astype", (dt,)),
+            None => Ok(arr),
+        }
+    }
+
+    fn __len__(&self) -> PyResult<usize> {
+        match self.shape.as_deref() {
+            Some([first, ..]) => Ok(*first as usize),
+            _ => Err(PyTypeError::new_err(
+                "Attempt to take len() of scalar dataset",
+            )),
+        }
+    }
+
+    fn __repr__(&self, py: Python<'_>) -> String {
+        let dtype = match &self.conv {
+            Ok(c) => c
+                .dtype
+                .bind(py)
+                .str()
+                .map(|s| s.to_string())
+                .unwrap_or_default(),
+            Err(_) => format!("{:?}", self.datatype),
+        };
+        let shape = match &self.shape {
+            Some(s) => format!("{s:?}"),
+            None => "None".to_string(),
+        };
         format!(
-            "<HDF5 Dataset \"{}\": shape {:?}, dtype {}>",
-            self.path,
-            self.cached_shape,
-            dtype_to_numpy_str(&self.cached_dtype),
+            "<HDF5 dataset \"{}\": shape {shape}, type \"{dtype}\">",
+            node::name(&self.path)
         )
     }
-
-    fn __len__(&self) -> usize {
-        self.cached_shape.first().copied().unwrap_or(0) as usize
-    }
-}
-
-impl PyDataset {
-    /// Read the full dataset and return it as a numpy array (or list for strings).
-    ///
-    /// For numeric types, the Rust I/O (file reading + decompression) is
-    /// performed inside `py.detach()` so that the GIL is released
-    /// during the potentially expensive operation.  The numpy array
-    /// construction still happens with the GIL held.
-    fn read_as_numpy<'py>(&self, py: Python<'py>) -> PyResult<Bound<'py, PyAny>> {
-        let file = &self.file;
-        let path = &self.path;
-        let shape: Vec<usize> = self.cached_shape.iter().map(|&d| d as usize).collect();
-
-        match &self.cached_dtype {
-            DType::F64 => {
-                let data = py
-                    .detach(|| file.dataset(path).and_then(|ds| ds.read_f64()))
-                    .map_err(to_py_err)?;
-                let nd = ArrayD::from_shape_vec(IxDyn(&shape), data)
-                    .map_err(|e| PyErr::new::<pyo3::exceptions::PyValueError, _>(e.to_string()))?;
-                let arr = PyArrayDyn::from_owned_array(py, nd);
-                Ok(arr.into_any())
-            }
-            DType::F32 => {
-                let data = py
-                    .detach(|| file.dataset(path).and_then(|ds| ds.read_f32()))
-                    .map_err(to_py_err)?;
-                let nd = ArrayD::from_shape_vec(IxDyn(&shape), data)
-                    .map_err(|e| PyErr::new::<pyo3::exceptions::PyValueError, _>(e.to_string()))?;
-                let arr = PyArrayDyn::from_owned_array(py, nd);
-                Ok(arr.into_any())
-            }
-            DType::I32 => {
-                let data = py
-                    .detach(|| file.dataset(path).and_then(|ds| ds.read_i32()))
-                    .map_err(to_py_err)?;
-                let nd = ArrayD::from_shape_vec(IxDyn(&shape), data)
-                    .map_err(|e| PyErr::new::<pyo3::exceptions::PyValueError, _>(e.to_string()))?;
-                let arr = PyArrayDyn::from_owned_array(py, nd);
-                Ok(arr.into_any())
-            }
-            DType::I64 => {
-                let data = py
-                    .detach(|| file.dataset(path).and_then(|ds| ds.read_i64()))
-                    .map_err(to_py_err)?;
-                let nd = ArrayD::from_shape_vec(IxDyn(&shape), data)
-                    .map_err(|e| PyErr::new::<pyo3::exceptions::PyValueError, _>(e.to_string()))?;
-                let arr = PyArrayDyn::from_owned_array(py, nd);
-                Ok(arr.into_any())
-            }
-            DType::U8 => {
-                // Try zero-copy first (contiguous layout), fall back to
-                // read_u64 + cast for chunked/compact datasets.
-                let data: Vec<u8> = py
-                    .detach(|| {
-                        let ds = file.dataset(path)?;
-                        match ds.read_u8_zerocopy() {
-                            Ok(slice) => Ok(slice.to_vec()),
-                            Err(_) => {
-                                let raw = ds.read_u64()?;
-                                Ok(raw.iter().map(|&v| v as u8).collect())
-                            }
-                        }
-                    })
-                    .map_err(to_py_err)?;
-                let nd = ArrayD::from_shape_vec(IxDyn(&shape), data)
-                    .map_err(|e| PyErr::new::<pyo3::exceptions::PyValueError, _>(e.to_string()))?;
-                let arr = PyArrayDyn::from_owned_array(py, nd);
-                Ok(arr.into_any())
-            }
-            DType::U64 => {
-                let data = py
-                    .detach(|| file.dataset(path).and_then(|ds| ds.read_u64()))
-                    .map_err(to_py_err)?;
-                let nd = ArrayD::from_shape_vec(IxDyn(&shape), data)
-                    .map_err(|e| PyErr::new::<pyo3::exceptions::PyValueError, _>(e.to_string()))?;
-                let arr = PyArrayDyn::from_owned_array(py, nd);
-                Ok(arr.into_any())
-            }
-            DType::String | DType::VariableLengthString => {
-                // String reads need the GIL for PyList construction, but we
-                // release it during the Rust I/O portion.
-                let data = py
-                    .detach(|| file.dataset(path).and_then(|ds| ds.read_string()))
-                    .map_err(to_py_err)?;
-                let list = PyList::new(py, &data)?;
-                Ok(list.into_any())
-            }
-            other => Err(PyErr::new::<pyo3::exceptions::PyTypeError, _>(format!(
-                "unsupported dataset dtype for reading: {other}"
-            ))),
-        }
-    }
-}
-
-#[cfg(test)]
-mod tests {
-    use super::*;
-
-    #[test]
-    fn dtype_mapping() {
-        assert_eq!(dtype_to_numpy_str(&DType::F64), "float64");
-        assert_eq!(dtype_to_numpy_str(&DType::F32), "float32");
-        assert_eq!(dtype_to_numpy_str(&DType::I32), "int32");
-        assert_eq!(dtype_to_numpy_str(&DType::I64), "int64");
-        assert_eq!(dtype_to_numpy_str(&DType::U8), "uint8");
-        assert_eq!(dtype_to_numpy_str(&DType::String), "object");
-    }
-
-    #[test]
-    fn dataset_from_file() {
-        let mut b = clawhdf5_rs::FileBuilder::new();
-        b.create_dataset("vals").with_f64_data(&[1.0, 2.0, 3.0]);
-        let bytes = b.finish().unwrap();
-        let file = Arc::new(clawhdf5_rs::File::from_bytes(bytes).unwrap());
-        let ds = PyDataset::new(file, "vals".into()).unwrap();
-        assert_eq!(ds.cached_shape, vec![3]);
-        assert_eq!(ds.cached_dtype, DType::F64);
-    }
-
-    #[test]
-    fn dataset_len() {
-        let mut b = clawhdf5_rs::FileBuilder::new();
-        b.create_dataset("data")
-            .with_i32_data(&[10, 20, 30, 40])
-            .with_shape(&[2, 2]);
-        let bytes = b.finish().unwrap();
-        let file = Arc::new(clawhdf5_rs::File::from_bytes(bytes).unwrap());
-        let ds = PyDataset::new(file, "data".into()).unwrap();
-        assert_eq!(ds.__len__(), 2);
-    }
 }

crates/clawhdf5-py/src/file.rs

@@ -4,10 +4,10 @@ use std::path::PathBuf;
 use std::sync::{Arc, Mutex};
 
 use pyo3::prelude::*;
+use pyo3::types::PyList;
 
 use crate::attrs::PyAttrs;
-use crate::dataset::PyDataset;
-use crate::group::{PyGroup, WriteGroupState, finalize_write_group};
+use crate::group::{PyGroup, ReadGroup, WriteGroupState, finalize_write_group};
 use crate::{DatasetSpec, OwnedAttrValue, apply_dataset_spec, extract_numpy_data, to_py_err};
 
 /// Internal state for write mode.
@@ -35,10 +35,12 @@ struct WriteState {
 #[pyclass(name = "File")]
 pub struct PyFile {
     inner: Option<FileInner>,
+    filename: String,
 }
 
 enum FileInner {
-    Read(Arc<clawhdf5_rs::File>),
+    /// The root group; it holds the file.
+    Read(ReadGroup),
     Write(WriteState),
 }
 
@@ -51,15 +53,20 @@ impl PyFile {
     ///   mode: 'r' for read (default), 'w' for write
     #[new]
     #[pyo3(signature = (path, mode="r"))]
-    fn new(path: &str, mode: &str) -> PyResult<Self> {
+    fn new(py: Python<'_>, path: &str, mode: &str) -> PyResult<Self> {
+        let filename = path.to_string();
         match mode {
             "r" => {
-                let file = clawhdf5_rs::File::open(path).map_err(to_py_err)?;
+                let file = py.detach(|| {
+                    crate::no_panic(|| clawhdf5_rs::File::open(path).map_err(to_py_err))
+                })?;
                 Ok(Self {
-                    inner: Some(FileInner::Read(Arc::new(file))),
+                    inner: Some(FileInner::Read(root_group(Arc::new(file)))),
+                    filename,
                 })
             }
             "w" => Ok(Self {
+                filename,
                 inner: Some(FileInner::Write(WriteState {
                     path: PathBuf::from(path),
                     root_datasets: Vec::new(),
@@ -101,44 +108,51 @@ impl PyFile {
         Ok(false) // don't suppress exceptions
     }
 
-    /// Get a child object (dataset or group) by path.
+    /// Get a child object (dataset or group) by path; `f['/']` is the root.
     fn __getitem__(&self, py: Python<'_>, key: &str) -> PyResult<Py<PyAny>> {
-        let file = self.read_file()?;
-        // Try dataset first
-        match file.dataset(key) {
-            Ok(_) => {
-                let ds = PyDataset::new(Arc::clone(file), key.to_string())?;
-                Ok(ds.into_pyobject(py)?.into_any().unbind())
-            }
-            Err(clawhdf5_rs::Error::NotADataset(_)) => {
-                let grp = PyGroup::from_read(Arc::clone(file), key.to_string());
-                Ok(grp.into_pyobject(py)?.into_any().unbind())
-            }
-            Err(_) => {
-                // Could be a group (no DataLayout message, no error)
-                match file.group(key) {
-                    Ok(_) => {
-                        let grp = PyGroup::from_read(Arc::clone(file), key.to_string());
-                        Ok(grp.into_pyobject(py)?.into_any().unbind())
-                    }
-                    Err(e) => Err(PyErr::new::<pyo3::exceptions::PyKeyError, _>(format!(
-                        "{key}: {e}"
-                    ))),
-                }
-            }
+        self.read_file()?.get_item(py, key)
     }
+
+    /// `f.get(key, default=None)`.
+    #[pyo3(signature = (key, default=None))]
+    fn get(&self, py: Python<'_>, key: &str, default: Option<Py<PyAny>>) -> PyResult<Py<PyAny>> {
+        self.read_file()?.get(py, key, default)
     }
 
     /// List the names of all children in the root group.
     fn keys(&self, py: Python<'_>) -> PyResult<Py<PyAny>> {
-        let file = self.read_file()?;
-        let root = file.root();
-        let mut names = root.datasets().map_err(to_py_err)?;
-        let groups = root.groups().map_err(to_py_err)?;
-        names.extend(groups);
-        names.sort();
-        let list = pyo3::types::PyList::new(py, &names)?;
-        Ok(list.into_any().unbind())
+        let names = self.read_file()?.member_names()?;
+        Ok(PyList::new(py, names)?.into_any().unbind())
+    }
+
+    fn values(&self, py: Python<'_>) -> PyResult<Py<PyAny>> {
+        let vals = self.read_file()?.values(py)?;
+        Ok(PyList::new(py, vals)?.into_any().unbind())
+    }
+
+    fn items(&self, py: Python<'_>) -> PyResult<Py<PyAny>> {
+        let items = self.read_file()?.items(py)?;
+        Ok(PyList::new(py, items)?.into_any().unbind())
+    }
+
+    fn __iter__(&self, py: Python<'_>) -> PyResult<Py<PyAny>> {
+        self.keys(py)?.call_method0(py, "__iter__")
+    }
+
+    fn __len__(&self) -> PyResult<usize> {
+        Ok(self.read_file()?.member_names()?.len())
+    }
+
+    /// The root group's name, `/`.
+    #[getter]
+    fn name(&self) -> &'static str {
+        "/"
+    }
+
+    /// The path the file was opened with.
+    #[getter]
+    fn filename(&self) -> &str {
+        &self.filename
     }
 
     /// Create a dataset in the root group (write mode only).
@@ -192,10 +206,7 @@ impl PyFile {
     #[getter]
     fn attrs(&self) -> PyResult<PyAttrs> {
         match self.inner.as_ref() {
-            Some(FileInner::Read(file)) => {
-                let map = file.root().attrs().map_err(to_py_err)?;
-                Ok(PyAttrs::from_read(map))
-            }
+            Some(FileInner::Read(root)) => root.attrs(),
             Some(FileInner::Write(state)) => Ok(PyAttrs::from_write(Arc::clone(&state.root_attrs))),
             None => Err(PyErr::new::<pyo3::exceptions::PyIOError, _>(
                 "file is closed",
@@ -205,8 +216,8 @@ impl PyFile {
 
     fn __repr__(&self) -> String {
         match &self.inner {
-            Some(FileInner::Read(f)) => {
-                format!("<HDF5 File (read, {} bytes)>", f.as_bytes().len())
+            Some(FileInner::Read(root)) => {
+                format!("<HDF5 File (read, {} bytes)>", root.file.as_bytes().len())
             }
             Some(FileInner::Write(s)) => {
                 format!("<HDF5 File (write, \"{}\")>", s.path.display())
@@ -216,13 +227,13 @@ impl PyFile {
     }
 
     fn __contains__(&self, key: &str) -> PyResult<bool> {
-        let file = self.read_file()?;
-        Ok(file.dataset(key).is_ok() || file.group(key).is_ok())
+        Ok(self.read_file()?.contains(key))
     }
 }
 
 impl PyFile {
-    fn read_file(&self) -> PyResult<&Arc<clawhdf5_rs::File>> {
+    /// The root group of a file opened for reading.
+    fn read_file(&self) -> PyResult<&ReadGroup> {
         match &self.inner {
             Some(FileInner::Read(f)) => Ok(f),
             Some(FileInner::Write(_)) => Err(PyErr::new::<pyo3::exceptions::PyIOError, _>(
@@ -260,8 +271,14 @@ fn parse_compression(
     }
 }
 
+fn root_group(file: Arc<clawhdf5_rs::File>) -> ReadGroup {
+    let root = file.superblock().root_group_address;
+    ReadGroup::new(file, String::new(), root)
+}
+
 /// Build and write the HDF5 file from accumulated write state.
 fn finalize_write(state: WriteState) -> PyResult<()> {
+    crate::no_panic(|| {
         let mut builder = clawhdf5_rs::FileBuilder::new();
 
         // Root attributes
@@ -285,6 +302,7 @@ fn finalize_write(state: WriteState) -> PyResult<()> {
 
         builder.write(&state.path).map_err(to_py_err)?;
         Ok(())
+    })
 }
 
 #[cfg(test)]

crates/clawhdf5-py/src/group.rs

@@ -1,13 +1,14 @@
 //! PyGroup — navigable HDF5 group with read and write support.
 
-use std::sync::{Arc, Mutex};
+use std::collections::HashMap;
+use std::sync::{Arc, Mutex, OnceLock};
 
+use pyo3::exceptions::{PyIOError, PyKeyError, PyValueError};
 use pyo3::prelude::*;
 use pyo3::types::PyList;
 
 use crate::attrs::PyAttrs;
-use crate::dataset::PyDataset;
-use crate::{DatasetSpec, OwnedAttrValue, apply_dataset_spec, extract_numpy_data, to_py_err};
+use crate::{DatasetSpec, OwnedAttrValue, apply_dataset_spec, extract_numpy_data, node};
 
 /// Shared state for a group being written.
 pub(crate) struct WriteGroupState {
@@ -18,32 +19,24 @@ pub(crate) struct WriteGroupState {
 
 /// An HDF5 group.
 ///
-/// In read mode, provides `__getitem__` navigation and child listing.
-/// In write mode, supports `create_dataset` and `create_group` and
-/// attribute setting.
-///
-/// ```python
-/// grp = f['group_name']
-/// grp.keys()
-/// ds = grp['dataset']
-/// ```
+/// In read mode it behaves like an h5py group: `grp['name']`,
+/// `grp['sub/path']` and `grp['/absolute/path']`, `keys()`, `values()`,
+/// `items()`, iteration, `len()`, `in`, `get()`, `name` and `attrs`.
+/// In write mode, supports `create_dataset` and attribute setting.
 #[pyclass(name = "Group")]
 pub struct PyGroup {
     inner: GroupInner,
 }
 
 enum GroupInner {
-    Read {
-        file: Arc<clawhdf5_rs::File>,
-        path: String,
-    },
+    Read(ReadGroup),
     Write(Arc<Mutex<WriteGroupState>>),
 }
 
 impl PyGroup {
-    pub(crate) fn from_read(file: Arc<clawhdf5_rs::File>, path: String) -> Self {
+    pub(crate) fn from_read(file: Arc<clawhdf5_rs::File>, path: String, addr: u64) -> Self {
         Self {
-            inner: GroupInner::Read { file, path },
+            inner: GroupInner::Read(ReadGroup::new(file, path, addr)),
         }
     }
 
@@ -52,63 +45,172 @@ impl PyGroup {
             inner: GroupInner::Write(state),
         }
     }
+
+    fn read_group(&self, what: &str) -> PyResult<&ReadGroup> {
+        match &self.inner {
+            GroupInner::Read(g) => Ok(g),
+            GroupInner::Write(_) => Err(PyIOError::new_err(format!(
+                "cannot {what} a group opened for writing"
+            ))),
+        }
+    }
+}
+
+/// A group in a file opened for reading (a file is its root group, as in
+/// h5py). It keeps its own address and, once listed, its links, so looking
+/// up a child neither resolves the path from the root nor scans the group's
+/// links again: visiting every member of a large group is linear, not
+/// quadratic.
+pub(crate) struct ReadGroup {
+    pub file: Arc<clawhdf5_rs::File>,
+    pub path: String,
+    pub addr: u64,
+    /// Link name -> object address (soft links resolved), filled on first use.
+    links: OnceLock<HashMap<String, u64>>,
+    /// Names of the datasets and subgroups, sorted (h5py's order).
+    members: OnceLock<Vec<String>>,
+}
+
+impl ReadGroup {
+    pub(crate) fn new(file: Arc<clawhdf5_rs::File>, path: String, addr: u64) -> Self {
+        Self {
+            file,
+            path,
+            addr,
+            links: OnceLock::new(),
+            members: OnceLock::new(),
+        }
+    }
+
+    fn links(&self) -> PyResult<&HashMap<String, u64>> {
+        if let Some(links) = self.links.get() {
+            return Ok(links);
+        }
+        let entries = crate::no_panic(|| {
+            clawhdf5_format::group_v2::resolve_group_children(
+                self.file.as_bytes(),
+                self.file.superblock(),
+                self.addr,
+            )
+            .map_err(|e| PyValueError::new_err(format!("{}: {e}", node::name(&self.path))))
+        })?;
+        let map = entries
+            .into_iter()
+            .map(|e| (e.name, e.object_header_address))
+            .collect();
+        Ok(self.links.get_or_init(|| map))
+    }
+
+    /// The path and address of `key` (a name, a relative or an absolute path).
+    fn locate(&self, key: &str) -> PyResult<(String, u64)> {
+        let path = node::join(&self.path, key);
+        let rel = if self.path.is_empty() {
+            Some(path.as_str())
+        } else if path == self.path {
+            Some("")
+        } else {
+            path.strip_prefix(self.path.as_str())
+                .and_then(|r| r.strip_prefix('/'))
+        };
+        let addr = match rel {
+            // A direct child: the link table, when it has the name.
+            Some(name) if !name.is_empty() && !name.contains('/') => {
+                match self.links()?.get(name) {
+                    Some(&a) => a,
+                    None => node::resolve_from(&self.file, self.addr, name, &path)?,
+                }
+            }
+            Some(rel) => node::resolve_from(&self.file, self.addr, rel, &path)?,
+            None => node::address(&self.file, &path)?,
+        };
+        Ok((path, addr))
+    }
+
+    /// `group[key]`.
+    pub(crate) fn get_item(&self, py: Python<'_>, key: &str) -> PyResult<Py<PyAny>> {
+        let (path, addr) = self.locate(key)?;
+        node::open(py, &self.file, path, addr)
+    }
+
+    /// `group.get(key, default)`.
+    pub(crate) fn get(
+        &self,
+        py: Python<'_>,
+        key: &str,
+        default: Option<Py<PyAny>>,
+    ) -> PyResult<Py<PyAny>> {
+        match self.get_item(py, key) {
+            Err(e) if e.is_instance_of::<PyKeyError>(py) => {
+                Ok(default.unwrap_or_else(|| py.None()))
+            }
+            other => other,
+        }
+    }
+
+    /// Names of the group's datasets and subgroups, sorted (h5py's order).
+    pub(crate) fn member_names(&self) -> PyResult<&[String]> {
+        if let Some(m) = self.members.get() {
+            return Ok(m);
+        }
+        let mut names = Vec::new();
+        for (name, &addr) in self.links()? {
+            let hdr = node::header_at(&self.file, addr, &node::join(&self.path, name))?;
+            if matches!(
+                node::kind(&hdr),
+                Some(node::Kind::Dataset | node::Kind::Group)
+            ) {
+                names.push(name.clone());
+            }
+        }
+        names.sort_by(|a, b| a.as_bytes().cmp(b.as_bytes()));
+        Ok(self.members.get_or_init(|| names))
+    }
+
+    pub(crate) fn contains(&self, key: &str) -> bool {
+        self.locate(key)
+            .and_then(|(path, addr)| node::header_at(&self.file, addr, &path))
+            .ok()
+            .and_then(|h| node::kind(&h))
+            .is_some_and(|k| k != node::Kind::Datatype)
+    }
+
+    pub(crate) fn values(&self, py: Python<'_>) -> PyResult<Vec<Py<PyAny>>> {
+        self.member_names()?
+            .iter()
+            .map(|n| self.get_item(py, n))
+            .collect()
+    }
+
+    pub(crate) fn items(&self, py: Python<'_>) -> PyResult<Vec<(String, Py<PyAny>)>> {
+        self.member_names()?
+            .iter()
+            .map(|n| Ok((n.clone(), self.get_item(py, n)?)))
+            .collect()
+    }
+
+    pub(crate) fn attrs(&self) -> PyResult<PyAttrs> {
+        PyAttrs::read(Arc::clone(&self.file), self.addr, &self.path)
+    }
 }
 
 #[pymethods]
 impl PyGroup {
     /// Get a child object (dataset or subgroup) by name or path.
     fn __getitem__(&self, py: Python<'_>, key: &str) -> PyResult<Py<PyAny>> {
-        match &self.inner {
-            GroupInner::Read { file, path } => {
-                let full_path = if path.is_empty() {
-                    key.to_string()
-                } else {
-                    format!("{path}/{key}")
-                };
-                // Try dataset first
-                match file.dataset(&full_path) {
-                    Ok(_) => {
-                        let ds = PyDataset::new(Arc::clone(file), full_path)?;
-                        Ok(ds.into_pyobject(py)?.into_any().unbind())
-                    }
-                    Err(clawhdf5_rs::Error::NotADataset(_)) => {
-                        let grp = PyGroup::from_read(Arc::clone(file), full_path);
-                        Ok(grp.into_pyobject(py)?.into_any().unbind())
-                    }
-                    Err(e) => {
-                        // Could be a group without a DataLayout message
-                        match file.group(&full_path) {
-                            Ok(_) => {
-                                let grp = PyGroup::from_read(Arc::clone(file), full_path);
-                                Ok(grp.into_pyobject(py)?.into_any().unbind())
-                            }
-                            Err(_) => Err(PyErr::new::<pyo3::exceptions::PyKeyError, _>(format!(
-                                "{key}: {e}"
-                            ))),
-                        }
-                    }
-                }
-            }
-            GroupInner::Write(_) => Err(PyErr::new::<pyo3::exceptions::PyIOError, _>(
-                "cannot read children from a group opened for writing",
-            )),
+        self.read_group("read children from")?.get_item(py, key)
     }
+
+    /// `group.get(key, default=None)`.
+    #[pyo3(signature = (key, default=None))]
+    fn get(&self, py: Python<'_>, key: &str, default: Option<Py<PyAny>>) -> PyResult<Py<PyAny>> {
+        self.read_group("read children from")?.get(py, key, default)
     }
 
     /// List the names of all children (datasets and subgroups).
     fn keys(&self, py: Python<'_>) -> PyResult<Py<PyAny>> {
         match &self.inner {
-            GroupInner::Read { file, path } => {
-                let group = if path.is_empty() {
-                    file.root()
-                } else {
-                    file.group(path).map_err(to_py_err)?
-                };
-                let mut names = group.datasets().map_err(to_py_err)?;
-                let groups = group.groups().map_err(to_py_err)?;
-                names.extend(groups);
-                names.sort();
-                let list = PyList::new(py, &names)?;
+            GroupInner::Read(g) => {
+                let list = PyList::new(py, g.member_names()?)?;
                 Ok(list.into_any().unbind())
             }
             GroupInner::Write(state) => {
@@ -120,6 +222,36 @@ impl PyGroup {
         }
     }
 
+    fn values(&self, py: Python<'_>) -> PyResult<Py<PyAny>> {
+        let g = self.read_group("read children from")?;
+        Ok(PyList::new(py, g.values(py)?)?.into_any().unbind())
+    }
+
+    fn items(&self, py: Python<'_>) -> PyResult<Py<PyAny>> {
+        let g = self.read_group("read children from")?;
+        Ok(PyList::new(py, g.items(py)?)?.into_any().unbind())
+    }
+
+    fn __iter__(&self, py: Python<'_>) -> PyResult<Py<PyAny>> {
+        self.keys(py)?.call_method0(py, "__iter__")
+    }
+
+    fn __len__(&self) -> PyResult<usize> {
+        match &self.inner {
+            GroupInner::Read(g) => Ok(g.member_names()?.len()),
+            GroupInner::Write(state) => Ok(state.lock().unwrap().datasets.len()),
+        }
+    }
+
+    /// The group's full name, e.g. `/sensors`.
+    #[getter]
+    fn name(&self) -> String {
+        match &self.inner {
+            GroupInner::Read(g) => node::name(&g.path),
+            GroupInner::Write(state) => node::name(&state.lock().unwrap().name),
+        }
+    }
+
     /// Create a dataset inside this group (write mode only).
     ///
     /// Parameters:
@@ -161,7 +293,7 @@ impl PyGroup {
                 state.lock().unwrap().datasets.push(spec);
                 Ok(())
             }
-            GroupInner::Read { .. } => Err(PyErr::new::<pyo3::exceptions::PyIOError, _>(
+            GroupInner::Read { .. } => Err(PyIOError::new_err(
                 "cannot create datasets on a read-only group",
             )),
         }
@@ -171,15 +303,7 @@ impl PyGroup {
     #[getter]
     fn attrs(&self) -> PyResult<PyAttrs> {
         match &self.inner {
-            GroupInner::Read { file, path } => {
-                let group = if path.is_empty() {
-                    file.root()
-                } else {
-                    file.group(path).map_err(to_py_err)?
-                };
-                let map = group.attrs().map_err(to_py_err)?;
-                Ok(PyAttrs::from_read(map))
-            }
+            GroupInner::Read(g) => g.attrs(),
             GroupInner::Write(state) => {
                 let store = Arc::clone(&state.lock().unwrap().attrs);
                 Ok(PyAttrs::from_write(store))
@@ -189,12 +313,9 @@ impl PyGroup {
 
     fn __repr__(&self) -> String {
         match &self.inner {
-            GroupInner::Read { path, .. } => {
-                if path.is_empty() {
-                    "<HDF5 Group \"/\" (root)>".to_string()
-                } else {
-                    format!("<HDF5 Group \"/{path}\">")
-                }
+            GroupInner::Read(g) => {
+                let n = g.member_names().map_or(0, |m| m.len());
+                format!("<HDF5 group \"{}\" ({n} members)>", node::name(&g.path))
             }
             GroupInner::Write(state) => {
                 let name = &state.lock().unwrap().name;
@@ -205,14 +326,7 @@ impl PyGroup {
 
     fn __contains__(&self, key: &str) -> PyResult<bool> {
         match &self.inner {
-            GroupInner::Read { file, path } => {
-                let full_path = if path.is_empty() {
-                    key.to_string()
-                } else {
-                    format!("{path}/{key}")
-                };
-                Ok(file.dataset(&full_path).is_ok() || file.group(&full_path).is_ok())
-            }
+            GroupInner::Read(g) => Ok(g.contains(key)),
             GroupInner::Write(state) => {
                 let guard = state.lock().unwrap();
                 Ok(guard.datasets.iter().any(|d| d.name == key))
@@ -244,26 +358,28 @@ mod tests {
     use super::*;
 
     #[test]
-    fn read_group_construction() {
+    fn member_names_are_sorted() {
         let mut b = clawhdf5_rs::FileBuilder::new();
-        let mut g = b.create_group("grp");
+        b.create_dataset("zeta").with_f64_data(&[1.0]);
+        b.create_dataset("alpha").with_f64_data(&[1.0]);
+        let mut g = b.create_group("mid");
         g.create_dataset("x").with_f64_data(&[1.0]);
         let finished = g.finish();
         b.add_group(finished);
         let bytes = b.finish().unwrap();
         let file = Arc::new(clawhdf5_rs::File::from_bytes(bytes).unwrap());
-        let _grp = PyGroup::from_read(file, "grp".into());
-    }
-
-    #[test]
-    fn write_group_state() {
-        let state = WriteGroupState {
-            name: "test".into(),
-            datasets: vec![],
-            attrs: Arc::new(Mutex::new(vec![])),
-        };
-        let arc = Arc::new(Mutex::new(state));
-        let _grp = PyGroup::from_write(arc);
+        let root = file.superblock().root_group_address;
+        let top = ReadGroup::new(Arc::clone(&file), String::new(), root);
+        assert_eq!(top.member_names().unwrap(), ["alpha", "mid", "zeta"]);
+        let (path, addr) = top.locate("mid").unwrap();
+        assert_eq!(path, "mid");
+        let mid = ReadGroup::new(Arc::clone(&file), path, addr);
+        assert_eq!(mid.member_names().unwrap(), ["x"]);
+        assert!(top.contains("mid/x"));
+        assert!(mid.contains("/alpha"));
+        assert!(mid.contains("x") && mid.contains("./x"));
+        assert!(!top.contains("nope"));
+        assert!(!mid.contains("alpha"));
     }
 
     #[test]

crates/clawhdf5-py/src/lib.rs

@@ -10,9 +10,12 @@
 //! ```
 
 mod attrs;
+mod convert;
 mod dataset;
 mod file;
 mod group;
+mod node;
+mod select;
 
 use pyo3::prelude::*;
 
@@ -21,6 +24,42 @@ pub(crate) use dataset::PyDataset;
 pub(crate) use file::PyFile;
 pub(crate) use group::PyGroup;
 
+pyo3::create_exception!(
+    clawhdf5,
+    InternalError,
+    pyo3::exceptions::PyRuntimeError,
+    "A bug in clawhdf5 met while reading or writing a file (a Rust panic, \
+     caught). Derived from RuntimeError, so `except Exception` handles it."
+);
+
+/// The text of a caught panic.
+pub(crate) fn panic_text(payload: &(dyn std::any::Any + Send)) -> String {
+    payload
+        .downcast_ref::<&str>()
+        .map(|s| (*s).to_string())
+        .or_else(|| payload.downcast_ref::<String>().cloned())
+        .unwrap_or_else(|| "unknown panic".to_string())
+}
+
+/// Run `f`, turning a panic in the library into [`InternalError`] instead of
+/// PyO3's `PanicException` (a `BaseException`, which `except Exception`
+/// does not catch). Wraps every call into the library.
+pub(crate) fn no_panic<T>(f: impl FnOnce() -> PyResult<T>) -> PyResult<T> {
+    std::panic::catch_unwind(std::panic::AssertUnwindSafe(f)).unwrap_or_else(|p| {
+        Err(InternalError::new_err(format!(
+            "clawhdf5 internal error (please report it): {}",
+            panic_text(&*p)
+        )))
+    })
+}
+
+/// A test hook: panics inside [`no_panic`], so the tests can check that a
+/// library panic reaches Python as an ordinary exception.
+#[pyfunction]
+fn _panic_for_test() -> PyResult<()> {
+    no_panic(|| panic!("deliberate panic for the test suite"))
+}
+
 /// Convert a `clawhdf5_rs::Error` into a `PyErr`.
 ///
 /// Maps different error variants to more specific Python exception types:
@@ -46,6 +85,54 @@ pub(crate) fn to_py_err(e: clawhdf5_rs::Error) -> PyErr {
     }
 }
 
+/// The value of a dataset or attribute with a null dataspace: a type but no
+/// data. Mirrors `h5py.Empty`.
+#[pyclass(name = "Empty", frozen)]
+pub struct PyEmpty {
+    dtype: Py<PyAny>,
+}
+
+impl PyEmpty {
+    pub(crate) fn new(dtype: Py<PyAny>) -> Self {
+        Self { dtype }
+    }
+}
+
+#[pymethods]
+impl PyEmpty {
+    #[new]
+    fn py_new(py: Python<'_>, dtype: &Bound<'_, PyAny>) -> PyResult<Self> {
+        let dtype = py.import("numpy")?.getattr("dtype")?.call1((dtype,))?;
+        Ok(Self::new(dtype.unbind()))
+    }
+
+    #[getter]
+    fn dtype(&self, py: Python<'_>) -> Py<PyAny> {
+        self.dtype.clone_ref(py)
+    }
+
+    #[getter]
+    fn shape(&self, py: Python<'_>) -> Py<PyAny> {
+        py.None()
+    }
+
+    #[getter]
+    fn size(&self, py: Python<'_>) -> Py<PyAny> {
+        py.None()
+    }
+
+    fn __eq__(&self, py: Python<'_>, other: &Bound<'_, PyAny>) -> PyResult<bool> {
+        match other.cast::<PyEmpty>() {
+            Ok(o) => self.dtype.bind(py).eq(o.get().dtype.bind(py)),
+            Err(_) => Ok(false),
+        }
+    }
+
+    fn __repr__(&self, py: Python<'_>) -> PyResult<String> {
+        Ok(format!("Empty(dtype={})", self.dtype.bind(py).repr()?))
+    }
+}
+
 /// The data payload for a dataset being written.
 #[derive(Clone)]
 pub(crate) enum DatasetData {
@@ -219,10 +306,14 @@ pub(crate) fn extract_numpy_data(
 /// The clawhdf5 Python module.
 #[pymodule]
 fn clawhdf5(m: &Bound<'_, PyModule>) -> PyResult<()> {
+    m.add("__version__", env!("CARGO_PKG_VERSION"))?;
     m.add_class::<PyFile>()?;
     m.add_class::<PyDataset>()?;
     m.add_class::<PyGroup>()?;
     m.add_class::<PyAttrs>()?;
+    m.add_class::<PyEmpty>()?;
+    m.add("InternalError", m.py().get_type::<InternalError>())?;
+    m.add_function(wrap_pyfunction!(_panic_for_test, m)?)?;
     Ok(())
 }
 
@@ -232,9 +323,9 @@ mod tests {
 
     #[test]
     fn owned_attr_value_roundtrip() {
-        let val = OwnedAttrValue::F64(3.14);
+        let val = OwnedAttrValue::F64(2.5);
         let attr: clawhdf5_rs::AttrValue = val.into();
-        assert!(matches!(attr, clawhdf5_rs::AttrValue::F64(v) if (v - 3.14).abs() < 1e-10));
+        assert!(matches!(attr, clawhdf5_rs::AttrValue::F64(v) if (v - 2.5).abs() < 1e-10));
     }
 
     #[test]

crates/clawhdf5-py/src/node.rs

@@ -0,0 +1,215 @@
+//! Resolving paths to objects in a file opened for reading.
+
+use std::sync::Arc;
+
+use clawhdf5_format::attribute::AttributeMessage;
+use clawhdf5_format::dataspace::{Dataspace, DataspaceType};
+use clawhdf5_format::message_type::MessageType;
+use clawhdf5_format::object_header::ObjectHeader;
+use pyo3::exceptions::{PyKeyError, PyTypeError, PyValueError};
+use pyo3::prelude::*;
+
+use crate::dataset::PyDataset;
+use crate::group::PyGroup;
+
+/// Join `key` onto the group path `base` the way h5py does: an absolute key
+/// starts from the root, a relative one from `base`. Paths are kept without
+/// a leading `/`; the root is `""`.
+pub(crate) fn join(base: &str, key: &str) -> String {
+    let parts = if key.starts_with('/') {
+        key.split('/').collect::<Vec<_>>()
+    } else {
+        base.split('/').chain(key.split('/')).collect()
+    };
+    parts
+        .into_iter()
+        .filter(|p| !p.is_empty() && *p != ".")
+        .collect::<Vec<_>>()
+        .join("/")
+}
+
+/// The HDF5 name (`/a/b`) of a path.
+pub(crate) fn name(path: &str) -> String {
+    format!("/{path}")
+}
+
+/// The address of the object at `path`, resolved from the root group.
+pub(crate) fn address(file: &clawhdf5_rs::File, path: &str) -> PyResult<u64> {
+    resolve_from(file, file.superblock().root_group_address, path, path)
+}
+
+/// The address of `rel` resolved from the group at `group` (`full` is the
+/// resulting path, for the error message).
+pub(crate) fn resolve_from(
+    file: &clawhdf5_rs::File,
+    group: u64,
+    rel: &str,
+    full: &str,
+) -> PyResult<u64> {
+    if rel.is_empty() {
+        return Ok(group);
+    }
+    crate::no_panic(|| {
+        clawhdf5_format::group_v2::resolve_path_from(file.as_bytes(), file.superblock(), group, rel)
+            .map_err(|e| {
+                PyKeyError::new_err(format!(
+                    "Unable to open object (object '{}' doesn't exist): {e}",
+                    name(full)
+                ))
+            })
+    })
+}
+
+/// The object header at `addr` (the object at `path`).
+pub(crate) fn header_at(file: &clawhdf5_rs::File, addr: u64, path: &str) -> PyResult<ObjectHeader> {
+    crate::no_panic(|| {
+        let sb = file.superblock();
+        let at = usize::try_from(addr)
+            .map_err(|_| PyValueError::new_err(format!("{}: address out of range", name(path))))?;
+        ObjectHeader::parse(file.as_bytes(), at, sb.offset_size, sb.length_size)
+            .map_err(|e| PyValueError::new_err(format!("{}: {e}", name(path))))
+    })
+}
+
+/// What an object header describes.
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub(crate) enum Kind {
+    Dataset,
+    Group,
+    Datatype,
+}
+
+pub(crate) fn kind(hdr: &ObjectHeader) -> Option<Kind> {
+    let has = |t: MessageType| hdr.messages.iter().any(|m| m.msg_type == t);
+    if has(MessageType::DataLayout) {
+        Some(Kind::Dataset)
+    } else if has(MessageType::LinkInfo)
+        || has(MessageType::Link)
+        || has(MessageType::SymbolTable)
+        || has(MessageType::GroupInfo)
+    {
+        Some(Kind::Group)
+    } else if has(MessageType::Datatype) {
+        Some(Kind::Datatype)
+    } else {
+        None
+    }
+}
+
+/// Open the object at `addr` (whose path is `path`) as a `Dataset` or
+/// `Group`. Both keep the address, so later reads resolve nothing.
+pub(crate) fn open(
+    py: Python<'_>,
+    file: &Arc<clawhdf5_rs::File>,
+    path: String,
+    addr: u64,
+) -> PyResult<Py<PyAny>> {
+    let hdr = header_at(file, addr, &path)?;
+    match kind(&hdr) {
+        Some(Kind::Dataset) => Ok(PyDataset::open(py, Arc::clone(file), path, addr, &hdr)?
+            .into_pyobject(py)?
+            .into_any()
+            .unbind()),
+        Some(Kind::Group) => Ok(PyGroup::from_read(Arc::clone(file), path, addr)
+            .into_pyobject(py)?
+            .into_any()
+            .unbind()),
+        Some(Kind::Datatype) => Err(PyTypeError::new_err(format!(
+            "{}: committed (named) datatypes are not supported by clawhdf5",
+            name(&path)
+        ))),
+        None => Err(PyValueError::new_err(format!(
+            "{}: not a dataset, group or datatype",
+            name(&path)
+        ))),
+    }
+}
+
+/// The dataspace message of an object header.
+pub(crate) fn dataspace(file: &clawhdf5_rs::File, hdr: &ObjectHeader) -> PyResult<Dataspace> {
+    crate::no_panic(|| {
+        let sb = file.superblock();
+        let msg = hdr
+            .messages
+            .iter()
+            .find(|m| m.msg_type == MessageType::Dataspace)
+            .ok_or_else(|| PyValueError::new_err("object has no dataspace message"))?;
+        let data = clawhdf5_format::shared_message::message_data(
+            file.as_bytes(),
+            msg,
+            sb.offset_size,
+            sb.length_size,
+        )
+        .map_err(|e| PyValueError::new_err(e.to_string()))?;
+        Dataspace::parse(&data, sb.length_size).map_err(|e| PyValueError::new_err(e.to_string()))
+    })
+}
+
+/// The chunk shape of a chunked dataset (one entry per dataset dimension),
+/// or `None` for other layouts or a layout message that does not parse.
+pub(crate) fn chunk_shape(
+    file: &clawhdf5_rs::File,
+    hdr: &ObjectHeader,
+    rank: usize,
+) -> Option<Vec<u64>> {
+    let sb = file.superblock();
+    let msg = hdr
+        .messages
+        .iter()
+        .find(|m| m.msg_type == MessageType::DataLayout)?;
+    match clawhdf5_format::data_layout::DataLayout::parse(&msg.data, sb.offset_size, sb.length_size)
+        .ok()?
+    {
+        clawhdf5_format::data_layout::DataLayout::Chunked {
+            chunk_dimensions, ..
+        } if chunk_dimensions.len() >= rank => Some(
+            chunk_dimensions[..rank]
+                .iter()
+                .map(|&d| u64::from(d))
+                .collect(),
+        ),
+        _ => None,
+    }
+}
+
+pub(crate) fn is_null(space: &Dataspace) -> bool {
+    space.space_type == DataspaceType::Null
+}
+
+/// The attributes of the object at `addr` (whose path is `path`), sorted by
+/// name (h5py's order). Attributes whose messages cannot be parsed are left
+/// out, as the facade's `attrs()` does.
+pub(crate) fn attributes(
+    file: &clawhdf5_rs::File,
+    addr: u64,
+    path: &str,
+) -> PyResult<Vec<AttributeMessage>> {
+    let hdr = header_at(file, addr, path)?;
+    crate::no_panic(|| {
+        let sb = file.superblock();
+        let (mut attrs, _errors) = clawhdf5_format::attribute::extract_attributes_tolerant(
+            file.as_bytes(),
+            &hdr,
+            sb.offset_size,
+            sb.length_size,
+        )
+        .map_err(|e| PyValueError::new_err(format!("{}: {e}", name(path))))?;
+        attrs.sort_by(|a, b| a.name.as_bytes().cmp(b.name.as_bytes()));
+        Ok(attrs)
+    })
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn join_paths() {
+        assert_eq!(join("", "a"), "a");
+        assert_eq!(join("a", "b/c"), "a/b/c");
+        assert_eq!(join("a/b", "/x"), "x");
+        assert_eq!(join("a", "/"), "");
+        assert_eq!(join("", "/a//b/"), "a/b");
+        assert_eq!(join("a", "./b"), "a/b");
+    }
+}

crates/clawhdf5-py/src/select.rs

@@ -0,0 +1,554 @@
+//! h5py-style indexing (`ds[1, 2:10:3, ...]`) mapped onto hyperslab
+//! selections, so the library reads the selection rather than the whole
+//! dataset (it still decodes everything for large selections; see the
+//! facade's `Dataset::read_selection`).
+//!
+//! The rules and error messages follow h5py's `selections.py`: integers
+//! (negative from the end) drop their axis, slices must have a positive
+//! step, one `Ellipsis` fills the unmentioned axes, a single increasing list
+//! of integers may index one axis, and strings name compound fields.
+//! Everything else (`None`/`np.newaxis`, boolean masks, several index lists)
+//! is refused with the error h5py gives.
+
+use clawhdf5_format::selection::Selection;
+use pyo3::exceptions::{PyIndexError, PyTypeError, PyValueError};
+use pyo3::prelude::*;
+use pyo3::types::{PyEllipsis, PySlice, PyString, PyTuple};
+
+/// The selection along one axis.
+#[derive(Clone, Debug, PartialEq)]
+pub(crate) enum Axis {
+    /// A single index: the axis is dropped from the result.
+    Index(u64),
+    /// `start, start + step, ...`, `count` of them.
+    Slice { start: u64, step: u64, count: u64 },
+    /// Increasing, distinct indices.
+    List(Vec<u64>),
+}
+
+impl Axis {
+    fn len(&self) -> u64 {
+        match self {
+            Axis::Index(_) => 1,
+            Axis::Slice { count, .. } => *count,
+            Axis::List(v) => v.len() as u64,
+        }
+    }
+}
+
+/// A parsed index expression.
+#[derive(Clone, Debug, PartialEq)]
+pub(crate) struct Plan {
+    /// One entry per dataset axis.
+    pub axes: Vec<Axis>,
+    /// Compound field names to keep (empty: all).
+    pub fields: Vec<String>,
+    /// For a scalar dataset: `ds[()]` gives a scalar, `ds[...]` a 0-d array.
+    /// For other datasets: every axis was an integer, so h5py gives a scalar.
+    pub scalar: bool,
+}
+
+impl Plan {
+    /// The shape of the result.
+    pub fn out_shape(&self) -> Vec<usize> {
+        self.axes
+            .iter()
+            .filter(|a| !matches!(a, Axis::Index(_)))
+            .map(|a| a.len() as usize)
+            .collect()
+    }
+
+    /// The shape of the result before the integer-indexed axes are dropped
+    /// (they have length 1 here): the shape of the joined reads.
+    pub fn read_shape(&self) -> Vec<usize> {
+        self.axes.iter().map(|a| a.len() as usize).collect()
+    }
+
+    /// The axis indexed by a list, if any.
+    pub fn list_axis(&self) -> Option<usize> {
+        self.axes.iter().position(|a| matches!(a, Axis::List(_)))
+    }
+
+    /// Whether the selection is empty.
+    pub fn is_empty(&self) -> bool {
+        self.axes.iter().any(|a| a.len() == 0)
+    }
+
+    /// The hyperslab reads that make up this selection, each with the shape
+    /// of its block (index axes kept at length 1). More than one only when an
+    /// axis is indexed by a list; those are joined along `list_axis`
+    /// afterwards (`join_along`), after keeping each read's `pick` rows.
+    ///
+    /// A list is read in groups, each one hyperslab over a stretch of the
+    /// axis, not once per index: every read decodes the chunks it touches
+    /// (and lists the dataset's chunks), so a read per run of indices decoded
+    /// the same chunk again and again. For a chunked dataset (`chunk_len` is
+    /// the chunk's length along the list axis) a group ends only where a
+    /// whole chunk holds no selected index, so no chunk is decoded twice or
+    /// without need. Otherwise a group ends at a gap of more than
+    /// [`MAX_GAP_BYTES`] of unselected data.
+    pub fn reads(
+        &self,
+        dims: &[u64],
+        chunk_len: Option<u64>,
+        elem_size: usize,
+    ) -> (Vec<Read>, Option<usize>) {
+        let list_axis = self.list_axis();
+        let groups: Vec<&[u64]> = match list_axis.map(|i| &self.axes[i]) {
+            Some(Axis::List(idx)) => {
+                let row_bytes = self.row_bytes(elem_size);
+                group_indices(idx, |last, next| match chunk_len {
+                    Some(c) if c > 0 => next / c <= last / c + 1,
+                    _ => (next - last - 1).saturating_mul(row_bytes) <= MAX_GAP_BYTES,
+                })
+            }
+            _ => vec![&[]],
+        };
+        let mut out = Vec::with_capacity(groups.len());
+        for group in groups {
+            let (first, span) = match (group.first(), group.last()) {
+                (Some(&f), Some(&l)) => (f, l - f + 1),
+                _ => (0, 0),
+            };
+            let pick = (span != group.len() as u64)
+                .then(|| group.iter().map(|&i| (i - first) as usize).collect());
+            let mut start = Vec::with_capacity(dims.len());
+            let mut stride = Vec::with_capacity(dims.len());
+            let mut count = Vec::with_capacity(dims.len());
+            for axis in &self.axes {
+                let (s, st, c) = match axis {
+                    Axis::Index(i) => (*i, 1, 1),
+                    Axis::Slice { start, step, count } => (*start, *step, *count),
+                    Axis::List(_) => (first, 1, span),
+                };
+                start.push(s);
+                // A stride only matters between blocks; keep it >= 1.
+                stride.push(if c <= 1 { 1 } else { st });
+                count.push(c);
+            }
+            let shape: Vec<usize> = count.iter().map(|&c| c as usize).collect();
+            let whole = start.iter().all(|&s| s == 0)
+                && stride.iter().all(|&s| s == 1)
+                && count.as_slice() == dims;
+            let sel = if whole {
+                Selection::All
+            } else {
+                let block = vec![1; dims.len()];
+                Selection::Hyperslab {
+                    start,
+                    stride,
+                    count,
+                    block,
+                }
+            };
+            out.push(Read { sel, shape, pick });
+        }
+        (out, list_axis)
+    }
+
+    /// Bytes of one step along the list axis within a read's bounding box.
+    fn row_bytes(&self, elem_size: usize) -> u64 {
+        self.axes
+            .iter()
+            .map(|a| match a {
+                Axis::Slice { step, count, .. } if *count > 0 => (count - 1) * step + 1,
+                _ => 1,
+            })
+            .fold(elem_size as u64, u64::saturating_mul)
+    }
+}
+
+/// Unselected data a read of a non-chunked dataset copies through rather
+/// than start another read.
+pub(crate) const MAX_GAP_BYTES: u64 = 64 * 1024;
+
+/// One hyperslab read of a selection.
+#[derive(Clone, Debug, PartialEq)]
+pub(crate) struct Read {
+    pub sel: Selection,
+    /// The block's shape (index axes at length 1).
+    pub shape: Vec<usize>,
+    /// For a list: the positions along the list axis, within the block, to
+    /// keep (`None`: all of them).
+    pub pick: Option<Vec<usize>>,
+}
+
+/// Split increasing indices into groups; `joins(last, next)` says whether
+/// `next` extends the group whose last index is `last`.
+fn group_indices(idx: &[u64], joins: impl Fn(u64, u64) -> bool) -> Vec<&[u64]> {
+    let mut groups = Vec::new();
+    let mut from = 0;
+    for k in 1..idx.len() {
+        if !joins(idx[k - 1], idx[k]) {
+            groups.push(&idx[from..k]);
+            from = k;
+        }
+    }
+    if from < idx.len() {
+        groups.push(&idx[from..]);
+    }
+    groups
+}
+
+/// Keep the elements at positions `pick` along `axis` of a row-major block.
+pub(crate) fn gather_along(
+    bytes: &[u8],
+    shape: &[usize],
+    axis: usize,
+    pick: &[usize],
+    elem_size: usize,
+) -> Vec<u8> {
+    let outer: usize = shape[..axis].iter().product();
+    let inner: usize = shape[axis + 1..].iter().product::<usize>() * elem_size;
+    let len = shape[axis];
+    let mut out = Vec::with_capacity(outer * pick.len() * inner);
+    for o in 0..outer {
+        for &p in pick {
+            let at = (o * len + p) * inner;
+            out.extend_from_slice(&bytes[at..at + inner]);
+        }
+    }
+    out
+}
+
+/// Join row-major blocks of `elem_size`-byte elements whose shapes differ
+/// only along `axis` into one buffer, in order along that axis. Whole
+/// elements are copied, so compound padding keeps the bytes that were read.
+pub(crate) fn join_along(
+    blocks: &[(Vec<u8>, Vec<usize>)],
+    axis: usize,
+    elem_size: usize,
+) -> Vec<u8> {
+    let Some((_, first)) = blocks.first() else {
+        return Vec::new();
+    };
+    let outer: usize = first[..axis].iter().product();
+    let inner: usize = first[axis + 1..].iter().product::<usize>() * elem_size;
+    let total: usize = blocks.iter().map(|(_, s)| s[axis]).sum();
+    let mut out = vec![0u8; outer * total * inner];
+    let mut at = 0;
+    for (bytes, shape) in blocks {
+        let len = shape[axis] * inner;
+        for o in 0..outer {
+            let dst = (o * total) * inner + at;
+            out[dst..dst + len].copy_from_slice(&bytes[o * len..(o + 1) * len]);
+        }
+        at += len;
+    }
+    out
+}
+
+/// Parse `key` for a dataset of shape `dims`.
+pub(crate) fn parse(key: &Bound<'_, PyAny>, dims: &[u64]) -> PyResult<Plan> {
+    let items: Vec<Bound<'_, PyAny>> = match key.cast::<PyTuple>() {
+        Ok(t) => t.iter().collect(),
+        Err(_) => vec![key.clone()],
+    };
+    let mut fields = Vec::new();
+    let mut args = Vec::new();
+    for item in items {
+        if let Ok(s) = item.cast::<PyString>() {
+            fields.push(s.to_str()?.to_owned());
+        } else {
+            args.push(item);
+        }
+    }
+
+    if args.iter().any(|a| a.is_none()) {
+        return Err(PyTypeError::new_err(
+            "Indexing with None (or np.newaxis) is not supported",
+        ));
+    }
+
+    let rank = dims.len();
+    if rank == 0 {
+        return match args.as_slice() {
+            [] => Ok(Plan {
+                axes: vec![],
+                fields,
+                scalar: true,
+            }),
+            [a] if a.is_instance_of::<PyEllipsis>() => Ok(Plan {
+                axes: vec![],
+                fields,
+                scalar: false,
+            }),
+            _ => Err(PyValueError::new_err(
+                "Illegal slicing argument for scalar dataspace",
+            )),
+        };
+    }
+
+    // Expand the ellipsis (at most one) to full slices.
+    let n_ellipsis = args
+        .iter()
+        .filter(|a| a.is_instance_of::<PyEllipsis>())
+        .count();
+    if n_ellipsis > 1 {
+        return Err(PyValueError::new_err("Only one ellipsis may be used."));
+    }
+    let explicit = args.len() - n_ellipsis;
+    if explicit > rank {
+        return Err(PyValueError::new_err(format!(
+            "{explicit} indexing arguments for {rank} dimensions"
+        )));
+    }
+    let py = key.py();
+    let mut expanded: Vec<Option<Bound<'_, PyAny>>> = Vec::with_capacity(rank);
+    for a in args {
+        if a.is_instance_of::<PyEllipsis>() {
+            for _ in 0..(rank - explicit) {
+                expanded.push(None);
+            }
+        } else {
+            expanded.push(Some(a));
+        }
+    }
+    while expanded.len() < rank {
+        expanded.push(None);
+    }
+
+    let mut axes = Vec::with_capacity(rank);
+    for (arg, &n) in expanded.iter().zip(dims) {
+        axes.push(match arg {
+            None => Axis::Slice {
+                start: 0,
+                step: 1,
+                count: n,
+            },
+            Some(a) => parse_axis(py, a, n)?,
+        });
+    }
+    if axes.iter().filter(|a| matches!(a, Axis::List(_))).count() > 1 {
+        return Err(PyTypeError::new_err(
+            "Only one indexing vector or array is currently allowed for fancy indexing",
+        ));
+    }
+    let scalar = axes.iter().all(|a| matches!(a, Axis::Index(_)));
+    Ok(Plan {
+        axes,
+        fields,
+        scalar,
+    })
+}
+
+fn parse_axis(py: Python<'_>, a: &Bound<'_, PyAny>, n: u64) -> PyResult<Axis> {
+    if a.is_none() {
+        return Err(PyTypeError::new_err(
+            "Indexing with None (or np.newaxis) is not supported",
+        ));
+    }
+    if let Ok(s) = a.cast::<PySlice>() {
+        let n_isize = isize::try_from(n)
+            .map_err(|_| PyValueError::new_err("dimension too large to slice"))?;
+        let ind = s.indices(n_isize)?;
+        if ind.step < 1 {
+            return Err(PyValueError::new_err(format!(
+                "Step must be >= 1 (got {})",
+                ind.step
+            )));
+        }
+        // `slicelength` is the number of elements selected, >= 0.
+        let count = ind.slicelength as u64;
+        let start = if count == 0 { 0 } else { ind.start as u64 };
+        return Ok(Axis::Slice {
+            start,
+            step: ind.step as u64,
+            count,
+        });
+    }
+    let np = py.import("numpy")?;
+    let is_bool =
+        a.is_instance_of::<pyo3::types::PyBool>() || a.is_instance(&np.getattr("bool_")?)?;
+    let is_array_like = a.is_instance(&np.getattr("ndarray")?)?
+        || a.is_instance_of::<pyo3::types::PyList>()
+        || a.is_instance_of::<PyTuple>();
+    // A 0-d integer array (`ds[np.array(1)]`) is an integer index, as in h5py.
+    if a.is_instance(&np.getattr("ndarray")?)? && a.getattr("ndim")?.extract::<usize>()? == 0 {
+        let kind: String = a.getattr("dtype")?.getattr("kind")?.extract()?;
+        if kind == "i" || kind == "u" {
+            let i: i128 = a.call_method0("item")?.extract()?;
+            return Ok(Axis::Index(normalize(i, n)?));
+        }
+    }
+    if !is_bool && !is_array_like && a.hasattr("__index__")? {
+        let i: i128 = a.call_method0("__index__")?.extract()?;
+        return Ok(Axis::Index(normalize(i, n)?));
+    }
+    if is_array_like {
+        let arr = np.call_method1("asarray", (a,))?;
+        let kind: String = arr.getattr("dtype")?.getattr("kind")?.extract()?;
+        if kind == "b" {
+            return Err(PyTypeError::new_err(
+                "Boolean mask indexing is not supported by clawhdf5",
+            ));
+        }
+        let ndim: usize = arr.getattr("ndim")?.extract()?;
+        let size: usize = arr.getattr("size")?.extract()?;
+        if size > 0 && kind != "i" && kind != "u" {
+            return Err(PyTypeError::new_err(
+                "Indexing arrays must have integer dtypes",
+            ));
+        }
+        if ndim > 1 {
+            return Err(PyTypeError::new_err(
+                "Only 1-D integer lists or arrays can be used for fancy indexing",
+            ));
+        }
+        let vals: Vec<i128> = arr.call_method0("tolist")?.extract()?;
+        let mut idx = Vec::with_capacity(vals.len());
+        for v in vals {
+            idx.push(normalize(v, n)?);
+        }
+        if idx.windows(2).any(|w| w[0] >= w[1]) {
+            return Err(PyTypeError::new_err(
+                "Indexing elements must be in increasing order",
+            ));
+        }
+        return Ok(Axis::List(idx));
+    }
+    Err(PyTypeError::new_err(format!(
+        "Illegal index type for clawhdf5 datasets: {}",
+        a.get_type().name()?
+    )))
+}
+
+fn normalize(i: i128, n: u64) -> PyResult<u64> {
+    let n_i = i128::from(n);
+    let j = if i < 0 { i + n_i } else { i };
+    if j < 0 || j >= n_i {
+        let hi = n_i - 1;
+        return Err(PyIndexError::new_err(format!(
+            "Index ({i}) out of range for (0-{hi})"
+        )));
+    }
+    Ok(j as u64)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn indices_group_by_chunk() {
+        let chunked = |c: u64| move |last: u64, next: u64| next / c <= last / c + 1;
+        // Chunks of 10: 3, 5 and 15 are in neighbouring chunks; 42 skips two.
+        let idx = [3, 5, 15, 42, 43, 99];
+        assert_eq!(
+            group_indices(&idx, chunked(10)),
+            vec![&[3, 5, 15][..], &[42, 43], &[99]]
+        );
+        assert_eq!(group_indices(&[], chunked(10)), Vec::<&[u64]>::new());
+    }
+
+    #[test]
+    fn a_list_reads_once_per_group() {
+        let plan = Plan {
+            axes: vec![
+                Axis::List(vec![0, 2, 3, 40]),
+                Axis::Slice {
+                    start: 0,
+                    step: 1,
+                    count: 5,
+                },
+            ],
+            fields: vec![],
+            scalar: false,
+        };
+        // Chunks of 8 rows: rows 0-3 are one read, row 40 another.
+        let (reads, axis) = plan.reads(&[50, 5], Some(8), 4);
+        assert_eq!(axis, Some(0));
+        assert_eq!(reads.len(), 2);
+        assert_eq!(reads[0].shape, vec![4, 5]);
+        assert_eq!(reads[0].pick, Some(vec![0, 2, 3]));
+        assert_eq!(reads[1].shape, vec![1, 5]);
+        assert_eq!(reads[1].pick, None);
+        // Not chunked: a gap under MAX_GAP_BYTES is read through.
+        let (reads, _) = plan.reads(&[50, 5], None, 4);
+        assert_eq!(reads.len(), 1);
+        assert_eq!(reads[0].pick, Some(vec![0, 2, 3, 40]));
+    }
+
+    #[test]
+    fn gather_keeps_picked_rows() {
+        // A 2x3 block of 1-byte elements; keep columns 0 and 2.
+        let block = [1, 2, 3, 4, 5, 6];
+        assert_eq!(
+            gather_along(&block, &[2, 3], 1, &[0, 2], 1),
+            vec![1, 3, 4, 6]
+        );
+    }
+
+    #[test]
+    fn blocks_join_along_the_list_axis() {
+        // Two 2x1 and 2x2 blocks of 1-byte elements, joined along axis 1.
+        let a = (vec![1, 2], vec![2, 1]);
+        let b = (vec![3, 4, 5, 6], vec![2, 2]);
+        assert_eq!(join_along(&[a, b], 1, 1), vec![1, 3, 4, 2, 5, 6]);
+        // Along axis 0 it is concatenation; 2-byte elements stay whole.
+        let a = (vec![1, 2, 3, 4], vec![1, 2]);
+        let b = (vec![5, 6, 7, 8], vec![1, 2]);
+        assert_eq!(join_along(&[a, b], 0, 2), vec![1, 2, 3, 4, 5, 6, 7, 8]);
+    }
+
+    #[test]
+    fn full_selection_reads_everything() {
+        let plan = Plan {
+            axes: vec![
+                Axis::Slice {
+                    start: 0,
+                    step: 1,
+                    count: 4,
+                },
+                Axis::Slice {
+                    start: 0,
+                    step: 1,
+                    count: 3,
+                },
+            ],
+            fields: vec![],
+            scalar: false,
+        };
+        let (reads, list) = plan.reads(&[4, 3], None, 8);
+        assert_eq!(list, None);
+        assert_eq!(
+            reads,
+            vec![Read {
+                sel: Selection::All,
+                shape: vec![4, 3],
+                pick: None
+            }]
+        );
+    }
+
+    #[test]
+    fn index_and_step_map_to_a_hyperslab() {
+        let plan = Plan {
+            axes: vec![
+                Axis::Index(2),
+                Axis::Slice {
+                    start: 1,
+                    step: 3,
+                    count: 2,
+                },
+            ],
+            fields: vec![],
+            scalar: false,
+        };
+        let (reads, _) = plan.reads(&[4, 8], None, 8);
+        assert_eq!(
+            reads,
+            vec![Read {
+                sel: Selection::Hyperslab {
+                    start: vec![2, 1],
+                    stride: vec![1, 3],
+                    count: vec![1, 2],
+                    block: vec![1, 1],
+                },
+                shape: vec![1, 2],
+                pick: None
+            }]
+        );
+        assert_eq!(plan.out_shape(), vec![2]);
+    }
+}

crates/clawhdf5-py/tests/conftest.py

@@ -0,0 +1,18 @@
+"""Shared fixtures for the clawhdf5 Python binding tests."""
+
+import os
+
+import pytest
+
+
+@pytest.fixture(scope="session")
+def h5py():
+    """h5py, or a skip — unless CLAWHDF5_REQUIRE_INTEROP=1, which makes a
+    missing h5py a failure (as for the Rust interop suites)."""
+    try:
+        import h5py as mod
+    except ImportError:
+        if os.environ.get("CLAWHDF5_REQUIRE_INTEROP") == "1":
+            pytest.fail("h5py is required (CLAWHDF5_REQUIRE_INTEROP=1) but not importable")
+        pytest.skip("h5py not installed")
+    return mod

crates/clawhdf5-py/tests/test_read_vs_h5py.py

@@ -0,0 +1,653 @@
+"""Every read through clawhdf5 compared against h5py (libhdf5) on files h5py
+writes: dtypes, shapes, values and the type of what comes back (array,
+numpy scalar, bytes, str, Empty), for every datatype the bindings map and a
+spread of index expressions; plus the errors h5py gives for the same keys."""
+
+import threading
+from concurrent.futures import ThreadPoolExecutor
+
+import numpy as np
+import pytest
+
+import clawhdf5
+
+# ---------------------------------------------------------------------------
+# The generated file
+# ---------------------------------------------------------------------------
+
+NUMERIC = [
+    "<i1", "<i2", "<i4", "<i8", "<u1", "<u2", "<u4", "<u8",
+    ">i2", ">i4", ">i8", ">u2", ">u4", ">u8",
+    "<f2", "<f4", "<f8", ">f2", ">f4", ">f8",
+]
+
+
+def _values(dtype, shape, seed):
+    rng = np.random.default_rng(seed)
+    dt = np.dtype(dtype)
+    n = int(np.prod(shape))
+    if dt.kind == "f":
+        return rng.standard_normal(n).astype(dt).reshape(shape)
+    info = np.iinfo(dt)
+    return rng.integers(info.min, info.max, size=n, dtype=dt.newbyteorder("=")).astype(dt).reshape(shape)
+
+
+def _compound_dtype():
+    return np.dtype(
+        {
+            "names": ["id", "pos", "label", "flag", "vec"],
+            "formats": ["<i4", ">f8", "S6", "u1", ("<f4", (3,))],
+            "offsets": [0, 8, 16, 22, 24],
+            "itemsize": 40,
+        }
+    )
+
+
+def _nested_dtype():
+    inner = np.dtype([("a", "<i2"), ("b", "<f4")])
+    return np.dtype([("x", "<u8"), ("inner", inner), ("c", "<c8")])
+
+
+def _write_fixture(h5py, path):
+    str_dt = h5py.string_dtype()
+    ascii_dt = h5py.string_dtype("ascii")
+    with h5py.File(path, "w") as f:
+        # Numeric types, both byte orders, 1-D contiguous and 2-D chunked+gzip.
+        for i, dt in enumerate(NUMERIC):
+            name = dt.replace("<", "le_").replace(">", "be_")
+            f.create_dataset(f"num/{name}_1d", data=_values(dt, (37,), i))
+            f.create_dataset(
+                f"num/{name}_2d_gzip",
+                data=_values(dt, (13, 11), 100 + i),
+                chunks=(4, 5),
+                compression="gzip",
+            )
+        f.create_dataset("num/f8_3d", data=_values("<f8", (6, 7, 5), 7), chunks=(2, 3, 5))
+        f.create_dataset(
+            "num/i4_3d_shuffle",
+            data=_values("<i4", (5, 9, 4), 8),
+            chunks=(3, 4, 2),
+            shuffle=True,
+            fletcher32=True,
+            compression="gzip",
+        )
+        f.create_dataset("num/scalar_f8", data=np.float64(3.25))
+        f.create_dataset("num/scalar_i2_be", data=np.array(-7, dtype=">i2"))
+        f.create_dataset("num/zero_size", shape=(0, 3), dtype="<f4")
+        f.create_dataset("num/empty", data=h5py.Empty("<f8"))
+        sparse = f.create_dataset("num/sparse_fill", shape=(40,), chunks=(8,), dtype="<i4", fillvalue=-3)
+        sparse[10:14] = [1, 2, 3, 4]
+        f.create_dataset("num/resizable", data=np.arange(12.0), maxshape=(None,), chunks=(5,))
+        # Compact layout (low level: h5py's create_dataset cannot ask for it).
+        dcpl = h5py.h5p.create(h5py.h5p.DATASET_CREATE)
+        dcpl.set_layout(h5py.h5d.COMPACT)
+        space = h5py.h5s.create_simple((9,))
+        dsid = h5py.h5d.create(f.id, b"num/compact_u2", h5py.h5t.py_create(np.dtype("<u2")), space, dcpl=dcpl)
+        dsid.write(h5py.h5s.ALL, h5py.h5s.ALL, np.arange(9, dtype="<u2") * 7)
+
+        # bool, enum, complex.
+        f.create_dataset("misc/bool", data=np.array([True, False, True, True, False]))
+        enum_dt = h5py.enum_dtype({"RED": 0, "GREEN": 1, "BLUE": 42}, basetype="i2")
+        f.create_dataset("misc/enum", data=np.array([0, 42, 1, 1, 0], dtype="<i2"), dtype=enum_dt)
+        uenum_dt = h5py.enum_dtype({"LO": 0, "HI": 2**40}, basetype="u8")
+        f.create_dataset("misc/enum_u8", data=np.array([0, 2**40], dtype="<u8"), dtype=uenum_dt)
+        f.create_dataset("misc/c8", data=(np.arange(6) + 1j * np.arange(6)).astype("<c8"))
+        f.create_dataset("misc/c16_2d", data=(np.arange(12) - 2j).reshape(3, 4).astype("<c16"))
+        f.create_dataset("misc/opaque", data=np.array([b"\x00\x01\x02", b"\xff\xfe\xfd"], dtype="V3"))
+
+        # Strings.
+        f.create_dataset("str/fixed", data=np.array([b"alpha", b"be", b"", b"gamma!"], dtype="S6"))
+        utf8_4 = h5py.string_dtype("utf-8", 4)
+        f.create_dataset(
+            "str/fixed_2d_utf8",
+            data=np.array([["é".encode(), b"b"], [b"c", b"dd"]], dtype=utf8_4),
+        )
+        f.create_dataset("str/vlen", data=["", "one", "twø", "a" * 300, "ünïcödé"], dtype=str_dt)
+        f.create_dataset("str/vlen_ascii", data=[b"x", b"yy", b"zzz"], dtype=ascii_dt)
+        f.create_dataset(
+            "str/vlen_2d_gzip",
+            data=np.array([[f"r{r}c{c}" * (r + c) for c in range(5)] for r in range(6)], dtype=object),
+            dtype=str_dt,
+            chunks=(2, 2),
+            compression="gzip",
+        )
+        f.create_dataset("str/vlen_scalar", data="just one", dtype=str_dt)
+
+        # Variable-length sequences.
+        vl_i = h5py.vlen_dtype(np.dtype("<i4"))
+        seqs = np.empty(4, dtype=object)
+        seqs[:] = [np.arange(3, dtype="<i4"), np.array([], dtype="<i4"), np.arange(10, dtype="<i4") * -1, np.array([7], dtype="<i4")]
+        f.create_dataset("vlen/i4", data=seqs, dtype=vl_i)
+        vl_f = h5py.vlen_dtype(np.dtype(">f8"))
+        fseqs = np.empty(3, dtype=object)
+        fseqs[:] = [np.linspace(0, 1, 5).astype(">f8"), np.array([2.5], dtype=">f8"), np.array([], dtype=">f8")]
+        f.create_dataset("vlen/f8_be", data=fseqs, dtype=vl_f)
+
+        # Compounds.
+        cdt = _compound_dtype()
+        rec = np.zeros(10, dtype=cdt)
+        rec["id"] = np.arange(10) * 3
+        rec["pos"] = np.linspace(-1, 1, 10)
+        rec["label"] = [f"n{i}".encode() for i in range(10)]
+        rec["flag"] = np.arange(10) % 2
+        rec["vec"] = np.arange(30, dtype="<f4").reshape(10, 3)
+        f.create_dataset("cmp/padded", data=rec)
+        f.create_dataset("cmp/padded_chunked", data=rec, chunks=(3,), compression="gzip")
+        ndt = _nested_dtype()
+        nrec = np.zeros((4, 3), dtype=ndt)
+        nrec["x"] = np.arange(12).reshape(4, 3)
+        nrec["inner"]["a"] = -np.arange(12).reshape(4, 3)
+        nrec["inner"]["b"] = np.arange(12).reshape(4, 3) / 4
+        nrec["c"] = np.arange(12).reshape(4, 3) * (1 + 1j)
+        f.create_dataset("cmp/nested_2d", data=nrec)
+        vdt = np.dtype([("n", "<i4"), ("s", h5py.string_dtype())])
+        vrec = np.array([(1, "a"), (2, "bb")], dtype=vdt)
+        f.create_dataset("cmp/with_vlen", data=vrec)
+
+        # A true HDF5 array datatype (h5py's high level would widen the shape).
+        tid = h5py.h5t.array_create(h5py.h5t.py_create(np.dtype("<i4")), (2, 3))
+        space = h5py.h5s.create_simple((4,))
+        dsid = h5py.h5d.create(f.id, b"cmp/array_type", tid, space)
+        dsid.write(h5py.h5s.ALL, h5py.h5s.ALL, np.arange(24, dtype="<i4").reshape(4, 2, 3), mtype=tid)
+
+        # Unsupported: references.
+        f.create_dataset("unsupported/refs", data=[f["num"].ref, f["misc"].ref], dtype=h5py.ref_dtype)
+
+        # Groups and attributes.
+        g = f.create_group("deep/er/est")
+        g.create_dataset("leaf", data=np.arange(5))
+        f.create_group("empty_group")
+        f.attrs["i4"] = np.int32(-5)
+        f.attrs["u8"] = np.uint64(2**63 + 1)
+        f.attrs["f2"] = np.float16(1.5)
+        f.attrs["f8_arr"] = np.array([1.0, 2.5, -3.0])
+        f.attrs["f8_one"] = np.array([4.0])
+        f.attrs["i2_2d_be"] = np.arange(6, dtype=">i2").reshape(2, 3)
+        f.attrs["bool"] = True
+        f.attrs["bool_arr"] = np.array([True, False])
+        f.attrs["vstr"] = "héllo"
+        f.attrs["vstr_arr"] = ["a", "bcd", ""]
+        f.attrs["fstr"] = np.bytes_(b"fixed")
+        f.attrs["fstr_arr"] = np.array([b"x", b"yz"], dtype="S2")
+        f.attrs["empty"] = h5py.Empty("<i4")
+        f.attrs["complex"] = np.complex128(1 - 2j)
+        f.attrs["compound"] = np.array((7, 2.5), dtype=[("a", "<i2"), ("b", "<f8")])
+        f.attrs["enum"] = np.array(42, dtype=enum_dt)
+        f["num/le_f8_1d"].attrs["units"] = "m/s"
+        f["num/le_f8_1d"].attrs["scale"] = np.float32(0.5)
+        g.attrs["depth"] = np.int8(3)
+
+
+@pytest.fixture(scope="module")
+def pair(h5py, tmp_path_factory):
+    path = str(tmp_path_factory.mktemp("h5") / "fixture.h5")
+    _write_fixture(h5py, path)
+    theirs = h5py.File(path, "r")
+    ours = clawhdf5.File(path, "r")
+    yield ours, theirs, path
+    theirs.close()
+    ours.close()
+
+
+def _all_datasets(h5py, f):
+    names = []
+    f.visititems(lambda n, o: names.append(n) if isinstance(o, h5py.Dataset) else None)
+    return sorted(names)
+
+
+# ---------------------------------------------------------------------------
+# Comparison helpers
+# ---------------------------------------------------------------------------
+
+
+def assert_same(ours, theirs, what=""):
+    if type(theirs).__name__ == "Empty":
+        assert isinstance(ours, clawhdf5.Empty), what
+        assert ours.dtype == theirs.dtype, what
+        return
+    assert type(ours) is type(theirs), f"{what}: {type(ours)} vs {type(theirs)}"
+    if isinstance(theirs, np.ndarray):
+        assert ours.shape == theirs.shape, what
+        assert ours.dtype == theirs.dtype, f"{what}: {ours.dtype} vs {theirs.dtype}"
+        if theirs.dtype == object:
+            for a, b in zip(ours.ravel(), theirs.ravel()):
+                assert_same(a, b, what)
+        elif theirs.dtype.kind == "V":
+            # Structured and opaque: every byte, padding included (h5py's
+            # padding is zero; uninitialised memory there would leak).
+            if theirs.dtype.names is not None:
+                np.testing.assert_array_equal(ours, theirs, err_msg=what)
+            assert ours.tobytes() == theirs.tobytes(), f"{what}: bytes differ"
+        else:
+            np.testing.assert_array_equal(ours, theirs, err_msg=what)
+    elif isinstance(theirs, np.generic):
+        assert ours.dtype == theirs.dtype, what
+        if theirs.dtype.names is not None:
+            np.testing.assert_array_equal(np.asarray(ours), np.asarray(theirs), err_msg=what)
+            assert ours.tobytes() == theirs.tobytes(), f"{what}: bytes differ"
+        else:
+            assert ours == theirs or (ours != ours and theirs != theirs), what
+    else:
+        assert ours == theirs, what
+
+
+def keys_for(shape):
+    if shape == ():
+        return [(), Ellipsis]
+    keys = [(), Ellipsis, 0, -1, slice(None), slice(None, None, 2), slice(1, None, 3), slice(0, 0), np.int64(0),
+            np.array(0), np.array(-1, dtype="i1")]
+    n0 = shape[0]
+    if n0 == 0:
+        return [(), Ellipsis, slice(None), slice(None, None, 2), slice(0, 0)]
+    keys += [slice(n0 // 2, None), slice(-3, None), [0, n0 - 1] if n0 > 1 else [0], (Ellipsis,)]
+    if n0 > 5:
+        keys += [[1, 2, 3, 5], [0, 4, 5]]
+    if len(shape) >= 2:
+        n1 = shape[1]
+        keys += [
+            (0, 0),
+            (-1, -1),
+            (slice(1, 3), slice(None, None, 2)),
+            (Ellipsis, 1),
+            (1, Ellipsis),
+            (slice(None), [0, n1 - 1] if n1 > 1 else [0]),
+            (slice(None, None, 2), 1),
+            (slice(0, 2), slice(3, 1)),
+            (np.array(1) if n0 > 1 else np.array(0), slice(None)),
+            (slice(None), np.array(n1 - 1, dtype="u2")),
+        ]
+    if len(shape) >= 3:
+        keys += [(0, slice(None), -1), (slice(1, None, 2), 2, slice(None, None, 3)), (Ellipsis, 0, 0), (0, Ellipsis, 1)]
+    return keys
+
+
+# h5py 3.16 (HDF5 2.0) returns the elements of a variable-length sequence of
+# big-endian floats unswapped (0.25 comes back as 2.6e-319), so these are
+# checked against the values written instead (test_vlen_big_endian).
+H5PY_MISREADS = {"vlen/f8_be"}
+
+ERROR_KEYS_1D = [
+    slice(None, None, -1), 10**6, -(10**6), None, (0, 0, 0, 0, 0), [3, 1], (Ellipsis, Ellipsis), 1.5, "nope",
+    np.array(1.0), np.array(True), np.array(10**6), [0, 0], [], (), Ellipsis, (0,), [-1], np.array([1, 2]),
+]
+
+
+# ---------------------------------------------------------------------------
+# Tests
+# ---------------------------------------------------------------------------
+
+
+def test_every_dataset_matches_h5py(h5py, pair):
+    ours, theirs, _ = pair
+    checked = 0
+    for name in _all_datasets(h5py, theirs):
+        if name.startswith("unsupported/") or name in H5PY_MISREADS or name == "cmp/with_vlen":
+            continue
+        t = theirs[name]
+        o = ours[name]
+        assert o.shape == t.shape, name
+        assert o.dtype == t.dtype, f"{name}: {o.dtype} vs {t.dtype}"
+        assert dict(o.dtype.metadata or {}) == dict(t.dtype.metadata or {}), name
+        assert o.ndim == t.ndim and o.size == t.size, name
+        assert o.maxshape == t.maxshape, name
+        assert o.name == t.name, name
+        if t.shape is None:
+            assert_same(o[()], t[()], name)
+            continue
+        for key in keys_for(t.shape):
+            what = f"{name}[{key!r}]"
+            try:
+                expected = t[key]
+            except Exception as e:  # noqa: BLE001 - h5py refuses: so must we
+                with pytest.raises(type(e)):
+                    o[key]
+                continue
+            assert_same(o[key], expected, what)
+            checked += 1
+    assert checked > 500
+
+
+def test_vlen_big_endian(pair):
+    ours, _, _ = pair
+    ds = ours["vlen/f8_be"]
+    assert ds.dtype.metadata["vlen"] == np.dtype(">f8")
+    got = ds[()]
+    expected = [np.linspace(0, 1, 5), np.array([2.5]), np.array([])]
+    assert got.shape == (3,)
+    for g, e in zip(got, expected):
+        assert g.dtype == np.dtype(">f8")
+        np.testing.assert_array_equal(g, e)
+    np.testing.assert_array_equal(ds[1], [2.5])
+
+
+def test_errors_match_h5py(h5py, pair):
+    ours, theirs, _ = pair
+    for name in ["num/le_i4_1d", "num/be_f8_2d_gzip", "str/vlen", "num/scalar_f8"]:
+        for key in ERROR_KEYS_1D:
+            try:
+                expected = theirs[name][key]
+            except Exception as e:  # noqa: BLE001
+                with pytest.raises(type(e)):
+                    ours[name][key]
+            else:
+                # h5py reads it, so must we (and the same values).
+                assert_same(ours[name][key], expected, f"{name}[{key!r}]")
+
+
+def test_compound_fields_match_h5py(pair):
+    ours, theirs, _ = pair
+    for name in ["cmp/padded", "cmp/padded_chunked", "cmp/nested_2d"]:
+        t, o = theirs[name], ours[name]
+        for field in t.dtype.names:
+            assert_same(o[field], t[field], f"{name}[{field}]")
+            assert_same(o[field, 1:3], t[field, 1:3], f"{name}[{field}, 1:3]")
+        two = list(t.dtype.names[:2])
+        expected = t[tuple(two)]
+        got = o[tuple(two)]
+        assert got.dtype.names == expected.dtype.names
+        for field in two:
+            np.testing.assert_array_equal(got[field], expected[field])
+    with pytest.raises(ValueError):
+        ours["cmp/padded"]["no_such_field"]
+    with pytest.raises(ValueError):
+        ours["num/le_i4_1d"]["id"]
+
+
+def test_numpy_asarray_and_len(pair):
+    ours, theirs, _ = pair
+    assert_same(np.asarray(ours["num/f8_3d"]), np.asarray(theirs["num/f8_3d"]))
+    assert len(ours["num/f8_3d"]) == len(theirs["num/f8_3d"])
+    with pytest.raises(TypeError):
+        len(ours["num/scalar_f8"])
+
+
+def test_attributes_match_h5py(h5py, pair):
+    ours, theirs, _ = pair
+    for path in ["/", "num/le_f8_1d", "deep/er/est"]:
+        t = theirs[path].attrs
+        o = ours[path].attrs
+        assert list(o.keys()) == sorted(t.keys()), path
+        assert len(o) == len(t)
+        for k in t.keys():
+            assert k in o
+            assert_same(o[k], t[k], f"{path}.attrs[{k}]")
+        assert o.get("missing", 5) == 5
+        with pytest.raises(KeyError):
+            o["missing"]
+        assert [k for k, _ in o.items()] == list(o.keys())
+
+
+def test_groups_match_h5py(h5py, pair):
+    ours, theirs, _ = pair
+    assert list(ours.keys()) == list(theirs.keys())
+    assert len(ours) == len(theirs)
+    for path in ["num", "deep", "deep/er", "deep/er/est", "empty_group"]:
+        o, t = ours[path], theirs[path]
+        assert isinstance(o, clawhdf5.Group)
+        assert list(o.keys()) == list(t.keys()), path
+        assert list(o) == list(t), path
+        assert len(o) == len(t), path
+        assert o.name == t.name
+    g = ours["deep/er"]
+    assert isinstance(g["est"], clawhdf5.Group)
+    assert isinstance(g["est/leaf"], clawhdf5.Dataset)
+    assert g["/deep/er/est/leaf"].name == "/deep/er/est/leaf"
+    assert "est/leaf" in g and "/num" in g and "nope" not in g
+    assert ours["/"].name == "/"
+    assert "deep/er/est/leaf" in ours
+    assert ours.get("nope") is None
+    with pytest.raises(KeyError):
+        ours["deep/nope"]
+    names = [k for k, v in ours["deep/er/est"].items()]
+    assert names == ["leaf"]
+    assert isinstance(ours["deep/er/est"].values()[0], clawhdf5.Dataset)
+
+
+def test_unsupported_types_are_errors_not_data(pair):
+    ours, _, _ = pair
+    ds = ours["unsupported/refs"]  # opening works
+    with pytest.raises(TypeError):
+        ds.dtype
+    with pytest.raises(TypeError):
+        ds[()]
+    with pytest.raises(TypeError):
+        ours["cmp/with_vlen"][()]
+
+
+def test_boolean_masks_are_refused(pair):
+    ours, _, _ = pair
+    with pytest.raises(TypeError):
+        ours["num/le_i4_1d"][np.ones(37, dtype=bool)]
+
+
+def test_reads_hand_numpy_the_rust_buffer(pair):
+    """A fixed-size read is a view over the buffer the library filled, not a
+    copy of it."""
+    ours, _, _ = pair
+    arr = ours["num/le_f8_2d_gzip"][2:9, 1:4]
+    assert not arr.flags.owndata
+    assert arr.base is not None
+    assert arr.flags.aligned and arr.flags.c_contiguous
+
+
+def test_only_the_selected_chunks_are_read(h5py, tmp_path):
+    """Damage one chunk: a selection that avoids it still reads, one that
+    touches it fails. Reading everything and slicing afterwards (what the
+    bindings did before) failed both. (The library reads the whole dataset
+    anyway when a selection's bounding box covers more than half of it, so
+    the selections here stay below that.)"""
+    path = str(tmp_path / "damaged.h5")
+    data = np.arange(1000, dtype="<f8")
+    with h5py.File(path, "w") as f:
+        f.create_dataset("d", data=data, chunks=(100,), compression="gzip")
+    with h5py.File(path, "r") as f:
+        info = f["d"].id.get_chunk_info(9)  # the last chunk
+    with open(path, "r+b") as fh:
+        fh.seek(info.byte_offset)
+        fh.write(b"\xff" * info.size)
+    with clawhdf5.File(path, "r") as f:
+        ds = f["d"]
+        np.testing.assert_array_equal(ds[0:400], data[0:400])
+        np.testing.assert_array_equal(ds[805:900], data[805:900])
+        np.testing.assert_array_equal(ds[5:450:7], data[5:450:7])
+        np.testing.assert_array_equal(ds[[3, 450, 899]], data[[3, 450, 899]])
+        with pytest.raises(Exception):
+            ds[950]
+        with pytest.raises(Exception):
+            ds[:]
+
+
+def test_threads_read_the_same_file(pair):
+    """Reads from many threads at once return exactly what h5py returns."""
+    ours, theirs, _ = pair
+    names = ["num/le_f8_2d_gzip", "num/i4_3d_shuffle", "str/vlen_2d_gzip", "cmp/padded_chunked", "num/be_i8_1d"]
+    expected = {n: theirs[n][()] for n in names}
+    errors = []
+    barrier = threading.Barrier(8)
+
+    def work(i):
+        barrier.wait()
+        for r in range(40):
+            n = names[(i + r) % len(names)]
+            got = ours[n][()]
+            try:
+                assert_same(got, expected[n], n)
+            except AssertionError as e:
+                errors.append(e)
+
+    with ThreadPoolExecutor(8) as pool:
+        list(pool.map(work, range(8)))
+    assert not errors, errors[:3]
+
+
+def test_reads_release_the_gil(h5py, tmp_path):
+    """While one thread is inside a long read, another Python thread keeps
+    running. With the GIL held for the read, the other thread would stall
+    for the whole read; the test measures its longest stall."""
+    import sys
+    import time
+
+    path = str(tmp_path / "gil.h5")
+    data = np.arange(2048 * 4096, dtype="<f4").reshape(2048, 4096)
+    with h5py.File(path, "w") as f:
+        f.create_dataset("d", data=data, chunks=(64, 4096), compression="gzip", compression_opts=1)
+    ds = clawhdf5.File(path, "r")["d"]
+    key = (slice(0, 1000), slice(None))  # under half: the uncached selection path
+    t0 = time.perf_counter()
+    ds[key]
+    one_read = time.perf_counter() - t0
+    assert one_read > 0.03, f"a read took only {one_read:.3f} s; too short to measure"
+
+    old = sys.getswitchinterval()
+    sys.setswitchinterval(0.001)
+    stop = threading.Event()
+    gaps = []
+
+    def spin():
+        last = time.perf_counter()
+        worst = 0.0
+        while not stop.is_set():
+            now = time.perf_counter()
+            worst = max(worst, now - last)
+            last = now
+        gaps.append(worst)
+
+    try:
+        t = threading.Thread(target=spin)
+        t.start()
+        time.sleep(0.01)
+        t0 = time.perf_counter()
+        for _ in range(2):
+            ds[key]
+        reading = time.perf_counter() - t0
+        stop.set()
+        t.join()
+    finally:
+        sys.setswitchinterval(old)
+    np.testing.assert_array_equal(ds[key], data[:1000])
+    # Held, the spinner would stall for about one read.
+    assert gaps[0] < one_read / 3, f"spinner stalled {gaps[0]:.3f} s during reads of {one_read:.3f} s ({reading:.3f} s)"
+
+
+def _v4_index_fixture(h5py, path):
+    """One 2-D dataset per v4 chunk index (HDF5 1.10+ layout, libver='latest')."""
+    data = (np.arange(37 * 23, dtype="<i4") * 7 - 1000).reshape(37, 23)
+    early = h5py.h5p.create(h5py.h5p.DATASET_CREATE)
+    early.set_alloc_time(h5py.h5d.ALLOC_TIME_EARLY)
+    with h5py.File(path, "w", libver="latest") as f:
+        f.create_dataset("implicit", data=data, chunks=(5, 4), dcpl=early)
+        f.create_dataset("fixed_array", data=data, chunks=(5, 4), compression="gzip")
+        f.create_dataset("extensible_array", data=data, chunks=(5, 4), maxshape=(None, 23), compression="gzip")
+        f.create_dataset("btree2", data=data, chunks=(5, 4), maxshape=(None, None), compression="gzip")
+        f.create_dataset("single_chunk", data=data, chunks=(37, 23), compression="gzip")
+
+
+V4_KEYS = [
+    slice(0, 3), slice(0, 30), (slice(7, 16), slice(3, 9)), (36, 22), (slice(None, None, 3), slice(1, None, 4)),
+    (slice(1, None, 2), Ellipsis), (Ellipsis, slice(2, 22)), [0, 5, 6, 36], (slice(None), [0, 3, 22]), -1, (),
+]
+
+
+def test_every_v4_chunk_index_matches_h5py(h5py, tmp_path):
+    """Partial reads of each v4 chunk index. The implicit index (early
+    allocation, no filters) used to panic in the library for any selection
+    covering more than half the dataset, e.g. ds[0:30]."""
+    path = str(tmp_path / "v4.h5")
+    _v4_index_fixture(h5py, path)
+    with h5py.File(path, "r") as theirs, clawhdf5.File(path, "r") as ours:
+        for name in theirs:
+            for key in V4_KEYS:
+                assert_same(ours[name][key], theirs[name][key], f"{name}[{key!r}]")
+
+
+def test_a_library_panic_is_an_ordinary_exception():
+    """PyO3 turns a Rust panic into PanicException, a BaseException that
+    `except Exception` does not catch. Every call into the library is
+    guarded, so a panic surfaces as clawhdf5.InternalError instead."""
+    assert issubclass(clawhdf5.InternalError, RuntimeError)
+    with pytest.raises(clawhdf5.InternalError, match="deliberate panic"):
+        clawhdf5._panic_for_test()
+    try:
+        clawhdf5._panic_for_test()
+    except Exception:  # noqa: BLE001 - the point of the test
+        pass
+
+
+def test_compound_padding_bytes_match_h5py(pair):
+    """Every byte of a padded compound, padding included, is h5py's, for
+    index lists with many runs as well as slices. Joining the runs with
+    np.concatenate left the padding uninitialised: process memory ended up
+    in tobytes()."""
+    ours, theirs, _ = pair
+    keys = [[0, 3, 6], [1, 2, 5, 9], [0, 2, 4, 6, 8], slice(None), slice(1, 9, 3), 4, [9]]
+    for name in ["cmp/padded", "cmp/padded_chunked"]:
+        for _ in range(20):  # garbage varies between runs; zeros do not
+            for key in keys:
+                assert ours[name][key].tobytes() == theirs[name][key].tobytes(), f"{name}[{key!r}]"
+    for key in [(slice(None), [0, 2]), ([0, 2, 3], slice(None)), ([1, 3], 1)]:
+        assert ours["cmp/nested_2d"][key].tobytes() == theirs["cmp/nested_2d"][key].tobytes(), key
+
+
+def test_a_long_index_list_decodes_each_chunk_once(h5py, tmp_path):
+    """An index list is read one group of chunks at a time, not one
+    hyperslab per run of indices: 5000 runs over 20 gzip chunks used to
+    decode the chunks 5000 times (8 s, against h5py's 0.014 s)."""
+    import time
+
+    path = str(tmp_path / "long_list.h5")
+    data = np.arange(200000, dtype="<f8")
+    grid = np.arange(400 * 3000, dtype="<i4").reshape(400, 3000)
+    with h5py.File(path, "w") as f:
+        f.create_dataset("d", data=data, chunks=(10000,), compression="gzip")
+        f.create_dataset("grid", data=grid, chunks=(50, 100), compression="gzip")
+        f.create_dataset("flat", data=data)  # contiguous
+    rng = np.random.default_rng(3)
+    cases = [
+        ("d", list(range(0, 200000, 40))),
+        ("d", sorted(rng.choice(200000, 3000, replace=False).tolist())),
+        ("flat", list(range(0, 200000, 40))),
+        ("flat", [0, 7, 199999]),
+        ("grid", (slice(None), list(range(0, 3000, 3)))),
+        ("grid", (sorted(rng.choice(400, 150, replace=False).tolist()), slice(5, 2900, 7))),
+        ("grid", (7, [0, 1, 2, 2000, 2999])),
+    ]
+    with h5py.File(path, "r") as theirs, clawhdf5.File(path, "r") as ours:
+        for name, key in cases:
+            t0 = time.perf_counter()
+            got = ours[name][key]
+            took = time.perf_counter() - t0
+            assert_same(got, theirs[name][key], f"{name}[{len(key)}-key]")
+            assert took < 2.0, f"{name}: {took:.2f} s"
+
+
+@pytest.mark.parametrize("libver", ["earliest", "latest"])
+def test_big_groups_are_not_quadratic(h5py, tmp_path, libver):
+    """A dataset or group remembers where its object is, and a group its
+    links, so reads and walks over a large group do not resolve every path
+    from the root again (it was O(n) per access: O(n^2) to visit a group)."""
+    import time
+
+    path = str(tmp_path / f"big_{libver}.h5")
+    n = 4000
+    with h5py.File(path, "w", libver=libver) as f:
+        g = f.create_group("g")
+        for i in range(n):
+            g.create_dataset(f"d{i:05d}", data=np.int32(i))
+        g.create_group("sub").create_dataset("leaf", data=np.arange(3))
+    with h5py.File(path, "r") as theirs, clawhdf5.File(path, "r") as ours:
+        t0 = time.perf_counter()
+        g = ours["g"]
+        assert list(g.keys()) == list(theirs["g"].keys())
+        total = sum(int(v[()]) for k, v in g.items() if k.startswith("d"))
+        assert total == n * (n - 1) // 2
+        seen = 0
+        for k in g:
+            if k.startswith("d"):
+                seen += int(g[k][()]) == int(k[1:])
+        assert seen == n
+        ds = g["d00007"]
+        assert all(ds[()] == 7 for _ in range(2000))
+        assert list(g["sub"]["leaf"][:]) == [0, 1, 2]
+        assert ours["/g/sub/leaf"][1] == 1 and g["/g/d00003"][()] == 3
+        took = time.perf_counter() - t0
+        assert took < 5.0, f"{took:.2f} s"

crates/clawhdf5-py/tests/test_write_read.py

@@ -1,4 +1,4 @@
-"""Tests for rustyhdf5 Python bindings."""
+"""Tests for clawhdf5 Python bindings."""
 
 import os
 import tempfile
@@ -6,7 +6,7 @@ import tempfile
 import numpy as np
 import pytest
 
-import rustyhdf5
+import clawhdf5
 
 
 @pytest.fixture
@@ -18,7 +18,7 @@ def tmp_h5(tmp_path):
 @pytest.fixture
 def sample_read_file(tmp_h5):
     """Create a sample HDF5 file for reading tests."""
-    with rustyhdf5.File(tmp_h5, "w") as f:
+    with clawhdf5.File(tmp_h5, "w") as f:
         f.create_dataset("temperatures", data=np.array([22.5, 23.1, 21.8]))
         f.create_dataset("counts", data=np.array([10, 20, 30], dtype=np.int32))
         f.attrs["version"] = 1
@@ -29,7 +29,7 @@ def sample_read_file(tmp_h5):
 @pytest.fixture
 def grouped_read_file(tmp_h5):
     """Create an HDF5 file with groups for reading tests."""
-    with rustyhdf5.File(tmp_h5, "w") as f:
+    with clawhdf5.File(tmp_h5, "w") as f:
         f.create_dataset("root_data", data=np.array([0.0, 1.0]))
         grp = f.create_group("sensors")
         grp.create_dataset("temperature", data=np.array([22.5, 23.1, 21.8]))
@@ -46,7 +46,7 @@ def grouped_read_file(tmp_h5):
 
 
 def test_open_and_read_f64(sample_read_file):
-    f = rustyhdf5.File(sample_read_file, "r")
+    f = clawhdf5.File(sample_read_file, "r")
     ds = f["temperatures"]
     data = ds[:]
     np.testing.assert_array_almost_equal(data, [22.5, 23.1, 21.8])
@@ -54,7 +54,7 @@ def test_open_and_read_f64(sample_read_file):
 
 
 def test_open_and_read_i32(sample_read_file):
-    f = rustyhdf5.File(sample_read_file, "r")
+    f = clawhdf5.File(sample_read_file, "r")
     ds = f["counts"]
     data = ds[:]
     np.testing.assert_array_equal(data, [10, 20, 30])
@@ -68,13 +68,13 @@ def test_open_and_read_i32(sample_read_file):
 
 
 def test_dataset_shape(sample_read_file):
-    with rustyhdf5.File(sample_read_file, "r") as f:
+    with clawhdf5.File(sample_read_file, "r") as f:
         ds = f["temperatures"]
         assert ds.shape == (3,)
 
 
 def test_dataset_dtype(sample_read_file):
-    with rustyhdf5.File(sample_read_file, "r") as f:
+    with clawhdf5.File(sample_read_file, "r") as f:
         assert f["temperatures"].dtype == "float64"
         assert f["counts"].dtype == "int32"
 
@@ -85,24 +85,24 @@ def test_dataset_dtype(sample_read_file):
 
 
 def test_read_root_attrs(sample_read_file):
-    with rustyhdf5.File(sample_read_file, "r") as f:
+    with clawhdf5.File(sample_read_file, "r") as f:
         assert f.attrs["version"] == 1
-        assert f.attrs["description"] == "test file"
+        assert f.attrs["description"] == b"test file"  # fixed-length string: numpy.bytes_, as in h5py
 
 
 def test_attrs_len(sample_read_file):
-    with rustyhdf5.File(sample_read_file, "r") as f:
+    with clawhdf5.File(sample_read_file, "r") as f:
         assert len(f.attrs) >= 2
 
 
 def test_attrs_contains(sample_read_file):
-    with rustyhdf5.File(sample_read_file, "r") as f:
+    with clawhdf5.File(sample_read_file, "r") as f:
         assert "version" in f.attrs
         assert "nonexistent" not in f.attrs
 
 
 def test_attrs_keys(sample_read_file):
-    with rustyhdf5.File(sample_read_file, "r") as f:
+    with clawhdf5.File(sample_read_file, "r") as f:
         keys = f.attrs.keys()
         assert "version" in keys
         assert "description" in keys
@@ -114,7 +114,7 @@ def test_attrs_keys(sample_read_file):
 
 
 def test_read_group_keys(grouped_read_file):
-    with rustyhdf5.File(grouped_read_file, "r") as f:
+    with clawhdf5.File(grouped_read_file, "r") as f:
         keys = f.keys()
         assert "sensors" in keys
         assert "metadata" in keys
@@ -122,7 +122,7 @@ def test_read_group_keys(grouped_read_file):
 
 
 def test_read_group_dataset(grouped_read_file):
-    with rustyhdf5.File(grouped_read_file, "r") as f:
+    with clawhdf5.File(grouped_read_file, "r") as f:
         grp = f["sensors"]
         ds = grp["temperature"]
         data = ds[:]
@@ -130,14 +130,14 @@ def test_read_group_dataset(grouped_read_file):
 
 
 def test_read_group_attrs(grouped_read_file):
-    with rustyhdf5.File(grouped_read_file, "r") as f:
+    with clawhdf5.File(grouped_read_file, "r") as f:
         grp = f["sensors"]
-        assert grp.attrs["location"] == "lab"
+        assert grp.attrs["location"] == b"lab"
 
 
 def test_nested_path_access(grouped_read_file):
     """Test f['group/dataset'] path navigation."""
-    with rustyhdf5.File(grouped_read_file, "r") as f:
+    with clawhdf5.File(grouped_read_file, "r") as f:
         ds = f["sensors/temperature"]
         data = ds[:]
         np.testing.assert_array_almost_equal(data, [22.5, 23.1, 21.8])
@@ -149,7 +149,7 @@ def test_nested_path_access(grouped_read_file):
 
 
 def test_context_manager(sample_read_file):
-    with rustyhdf5.File(sample_read_file, "r") as f:
+    with clawhdf5.File(sample_read_file, "r") as f:
         data = f["temperatures"][:]
         np.testing.assert_array_almost_equal(data, [22.5, 23.1, 21.8])
     # File should be closed after with block
@@ -162,30 +162,30 @@ def test_context_manager(sample_read_file):
 
 
 def test_write_simple(tmp_h5):
-    with rustyhdf5.File(tmp_h5, "w") as f:
+    with clawhdf5.File(tmp_h5, "w") as f:
         f.create_dataset("data", data=np.array([1.0, 2.0, 3.0]))
     # Verify by reading back
-    with rustyhdf5.File(tmp_h5, "r") as f:
+    with clawhdf5.File(tmp_h5, "r") as f:
         data = f["data"][:]
         np.testing.assert_array_almost_equal(data, [1.0, 2.0, 3.0])
 
 
 def test_write_with_attrs(tmp_h5):
-    with rustyhdf5.File(tmp_h5, "w") as f:
+    with clawhdf5.File(tmp_h5, "w") as f:
         f.create_dataset("values", data=np.array([10, 20], dtype=np.int32))
         f.attrs["author"] = "test"
         f.attrs["count"] = 42
-    with rustyhdf5.File(tmp_h5, "r") as f:
-        assert f.attrs["author"] == "test"
+    with clawhdf5.File(tmp_h5, "r") as f:
+        assert f.attrs["author"] == b"test"
         assert f.attrs["count"] == 42
 
 
 def test_write_with_group(tmp_h5):
-    with rustyhdf5.File(tmp_h5, "w") as f:
+    with clawhdf5.File(tmp_h5, "w") as f:
         grp = f.create_group("experiment")
         grp.create_dataset("results", data=np.array([3.14, 2.72]))
         grp.attrs["version"] = 1
-    with rustyhdf5.File(tmp_h5, "r") as f:
+    with clawhdf5.File(tmp_h5, "r") as f:
         ds = f["experiment/results"]
         np.testing.assert_array_almost_equal(ds[:], [3.14, 2.72])
         grp = f["experiment"]
@@ -199,9 +199,9 @@ def test_write_with_group(tmp_h5):
 
 def test_roundtrip_float64(tmp_h5):
     original = np.array([1.1, 2.2, 3.3], dtype=np.float64)
-    with rustyhdf5.File(tmp_h5, "w") as f:
+    with clawhdf5.File(tmp_h5, "w") as f:
         f.create_dataset("data", data=original)
-    with rustyhdf5.File(tmp_h5, "r") as f:
+    with clawhdf5.File(tmp_h5, "r") as f:
         result = f["data"][:]
         np.testing.assert_array_almost_equal(result, original)
         assert result.dtype == np.float64
@@ -209,9 +209,9 @@ def test_roundtrip_float64(tmp_h5):
 
 def test_roundtrip_float32(tmp_h5):
     original = np.array([1.5, 2.5, 3.5], dtype=np.float32)
-    with rustyhdf5.File(tmp_h5, "w") as f:
+    with clawhdf5.File(tmp_h5, "w") as f:
         f.create_dataset("data", data=original)
-    with rustyhdf5.File(tmp_h5, "r") as f:
+    with clawhdf5.File(tmp_h5, "r") as f:
         result = f["data"][:]
         np.testing.assert_array_almost_equal(result, original)
         assert result.dtype == np.float32
@@ -219,9 +219,9 @@ def test_roundtrip_float32(tmp_h5):
 
 def test_roundtrip_int32(tmp_h5):
     original = np.array([-10, 0, 10, 100], dtype=np.int32)
-    with rustyhdf5.File(tmp_h5, "w") as f:
+    with clawhdf5.File(tmp_h5, "w") as f:
         f.create_dataset("data", data=original)
-    with rustyhdf5.File(tmp_h5, "r") as f:
+    with clawhdf5.File(tmp_h5, "r") as f:
         result = f["data"][:]
         np.testing.assert_array_equal(result, original)
         assert result.dtype == np.int32
@@ -229,9 +229,9 @@ def test_roundtrip_int32(tmp_h5):
 
 def test_roundtrip_int64(tmp_h5):
     original = np.array([-1, 0, 1, 2**40], dtype=np.int64)
-    with rustyhdf5.File(tmp_h5, "w") as f:
+    with clawhdf5.File(tmp_h5, "w") as f:
         f.create_dataset("data", data=original)
-    with rustyhdf5.File(tmp_h5, "r") as f:
+    with clawhdf5.File(tmp_h5, "r") as f:
         result = f["data"][:]
         np.testing.assert_array_equal(result, original)
         assert result.dtype == np.int64
@@ -239,9 +239,9 @@ def test_roundtrip_int64(tmp_h5):
 
 def test_roundtrip_uint8(tmp_h5):
     original = np.array([0, 127, 255], dtype=np.uint8)
-    with rustyhdf5.File(tmp_h5, "w") as f:
+    with clawhdf5.File(tmp_h5, "w") as f:
         f.create_dataset("data", data=original)
-    with rustyhdf5.File(tmp_h5, "r") as f:
+    with clawhdf5.File(tmp_h5, "r") as f:
         result = f["data"][:]
         np.testing.assert_array_equal(result, original)
         assert result.dtype == np.uint8
@@ -254,7 +254,7 @@ def test_roundtrip_uint8(tmp_h5):
 
 def test_chunked_gzip(tmp_h5):
     original = np.arange(100, dtype=np.float64)
-    with rustyhdf5.File(tmp_h5, "w") as f:
+    with clawhdf5.File(tmp_h5, "w") as f:
         f.create_dataset(
             "compressed",
             data=original,
@@ -262,7 +262,7 @@ def test_chunked_gzip(tmp_h5):
             compression="gzip",
             compression_opts=6,
         )
-    with rustyhdf5.File(tmp_h5, "r") as f:
+    with clawhdf5.File(tmp_h5, "r") as f:
         result = f["compressed"][:]
         np.testing.assert_array_equal(result, original)
 
@@ -276,7 +276,7 @@ def test_h5py_can_read_our_file(tmp_h5):
     """Verify that h5py can read files we create."""
     import h5py
 
-    with rustyhdf5.File(tmp_h5, "w") as f:
+    with clawhdf5.File(tmp_h5, "w") as f:
         f.create_dataset("values", data=np.array([1.0, 2.0, 3.0]))
         f.attrs["meta"] = "hello"
     with h5py.File(tmp_h5, "r") as f:
@@ -292,7 +292,7 @@ def test_we_can_read_h5py_file(tmp_h5):
     with h5py.File(tmp_h5, "w") as f:
         f.create_dataset("data", data=np.array([10.0, 20.0, 30.0]))
         f.attrs["version"] = 2
-    with rustyhdf5.File(tmp_h5, "r") as f:
+    with clawhdf5.File(tmp_h5, "r") as f:
         data = f["data"][:]
         np.testing.assert_array_equal(data, [10.0, 20.0, 30.0])
         assert f.attrs["version"] == 2
@@ -305,9 +305,9 @@ def test_we_can_read_h5py_file(tmp_h5):
 
 def test_2d_array_roundtrip(tmp_h5):
     original = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]], dtype=np.float64)
-    with rustyhdf5.File(tmp_h5, "w") as f:
+    with clawhdf5.File(tmp_h5, "w") as f:
         f.create_dataset("matrix", data=original)
-    with rustyhdf5.File(tmp_h5, "r") as f:
+    with clawhdf5.File(tmp_h5, "r") as f:
         ds = f["matrix"]
         assert ds.shape == (2, 3)
         result = ds[:]
@@ -321,15 +321,15 @@ def test_2d_array_roundtrip(tmp_h5):
 
 def test_open_nonexistent_file():
     with pytest.raises(OSError):
-        rustyhdf5.File("/nonexistent/path.h5", "r")
+        clawhdf5.File("/nonexistent/path.h5", "r")
 
 
 def test_invalid_mode(tmp_h5):
     with pytest.raises(ValueError):
-        rustyhdf5.File(tmp_h5, "x")
+        clawhdf5.File(tmp_h5, "x")
 
 
 def test_key_error_on_missing_dataset(sample_read_file):
-    with rustyhdf5.File(sample_read_file, "r") as f:
+    with clawhdf5.File(sample_read_file, "r") as f:
         with pytest.raises(KeyError):
             f["nonexistent"]

crates/clawhdf5/src/reader.rs

@@ -176,6 +176,22 @@ impl File {
         })
     }
 
+    /// A `Dataset` handle for the object header at `address` (an address
+    /// from a group listing, or one kept from an earlier lookup), without
+    /// resolving a path. Resolving a path walks every group on it, which in
+    /// a large group costs a scan of its links; keep the address instead to
+    /// open the same dataset repeatedly.
+    pub fn dataset_at(&self, address: u64) -> Result<Dataset<'_>, Error> {
+        let hdr = self.parse_header(address)?;
+        if !has_message(&hdr, MessageType::DataLayout) {
+            return Err(Error::NotADataset(format!("object at address {address}")));
+        }
+        Ok(Dataset {
+            file: self,
+            header: hdr,
+        })
+    }
+
     /// Resolve a path and return a `Group` handle.
     ///
     /// The path uses `/` separators (e.g., `"sensors"`).
@@ -606,8 +622,16 @@ impl<'f> Dataset<'f> {
 
     /// Read selected elements as raw bytes.
     ///
-    /// Only the elements matching the [`clawhdf5_format::selection::Selection`] are returned. For chunked
-    /// datasets, only intersecting chunks are decompressed.
+    /// Only the elements matching the [`clawhdf5_format::selection::Selection`] are returned.
+    ///
+    /// What is read to get them: when the selection's bounding box covers at
+    /// most half the dataset, only that box — the chunks overlapping it, or
+    /// the rows of a contiguous dataset. The whole dataset is decoded instead
+    /// when the box covers more than half (a strided selection spanning the
+    /// dataset does), for compact and virtual layouts, for a dataset with no
+    /// storage, and for a chunked dataset with a non-default fill value.
+    /// [`Selection::All`](clawhdf5_format::selection::Selection::All) goes
+    /// through the file's chunk cache; other selections do not.
     pub fn read_selection(
         &self,
         selection: &clawhdf5_format::selection::Selection,

crates/clawhdf5/tests/integration_tests.rs

@@ -986,3 +986,35 @@ fn u64_data_roundtrip() {
         values
     );
 }
+
+// ---------------------------------------------------------------------------
+// Opening a dataset by address
+// ---------------------------------------------------------------------------
+
+#[test]
+fn dataset_at_opens_the_same_dataset_as_its_path() {
+    let mut b = FileBuilder::new();
+    let mut g = b.create_group("grp");
+    g.create_dataset("vals").with_f64_data(&[1.0, 2.5, -3.0]);
+    b.add_group(g.finish());
+    let file = File::from_bytes(b.finish().unwrap()).unwrap();
+
+    let addr =
+        clawhdf5_format::group_v2::resolve_path_any(file.as_bytes(), file.superblock(), "grp/vals")
+            .unwrap();
+    let by_addr = file.dataset_at(addr).unwrap();
+    assert_eq!(by_addr.read_f64().unwrap(), vec![1.0, 2.5, -3.0]);
+    assert_eq!(
+        by_addr.shape().unwrap(),
+        file.dataset("grp/vals").unwrap().shape().unwrap()
+    );
+
+    // The group's own header is not a dataset.
+    let group_addr =
+        clawhdf5_format::group_v2::resolve_path_any(file.as_bytes(), file.superblock(), "grp")
+            .unwrap();
+    assert!(matches!(
+        file.dataset_at(group_addr),
+        Err(clawhdf5::Error::NotADataset(_))
+    ));
+}

crates/clawhdf5/tests/v4_chunk_index_selection.rs

@@ -0,0 +1,241 @@
+//! Partial hyperslab reads of every v4 (`libver='latest'`) chunk index type,
+//! compared element for element with h5py.
+//!
+//! `Dataset::read_selection` takes two routes: `partial_read` materialises
+//! the selection's bounding box when it covers at most half the dataset, and
+//! `data_read::read_raw_data_selection` handles the rest. The second route
+//! once passed the layout's full chunk dimensions (which carry the element
+//! size as an extra, last dimension) to the implicit-index chunk generator,
+//! which then indexed past the dataset's rank and panicked. So every case
+//! below reads both a small window and one covering most of the dataset.
+//!
+//! Each case asserts which chunk index the file actually uses (parsed from
+//! the layout message), so a change in how h5py lays the file out can't turn
+//! this into a test of the wrong index.
+//!
+//! Skipped when python3 with h5py is unavailable, unless
+//! `CLAWHDF5_REQUIRE_INTEROP=1`.
+
+use std::process::Command;
+
+use clawhdf5::File;
+use clawhdf5_format::data_layout::DataLayout;
+use clawhdf5_format::message_type::MessageType;
+use clawhdf5_format::object_header::ObjectHeader;
+use clawhdf5_format::selection::Selection;
+use clawhdf5_format::superblock::Superblock;
+
+fn python() -> String {
+    std::env::var("CLAWHDF5_PYTHON").unwrap_or_else(|_| "python3".to_string())
+}
+
+fn interop_required() -> bool {
+    std::env::var("CLAWHDF5_REQUIRE_INTEROP").is_ok_and(|v| v == "1")
+}
+
+fn python_available() -> bool {
+    Command::new(python())
+        .args(["-c", "import h5py"])
+        .output()
+        .map(|o| o.status.success())
+        .unwrap_or(false)
+}
+
+fn run_python(script: &str) -> String {
+    let output = Command::new(python())
+        .args(["-c", script])
+        .output()
+        .expect("failed to run python");
+    assert!(
+        output.status.success(),
+        "Python script failed:\nSTDOUT: {}\nSTDERR: {}",
+        String::from_utf8_lossy(&output.stdout),
+        String::from_utf8_lossy(&output.stderr)
+    );
+    String::from_utf8_lossy(&output.stdout).trim().to_string()
+}
+
+/// The v4 chunk index type recorded in `name`'s layout message
+/// (1 single chunk, 2 implicit, 3 fixed array, 4 extensible array, 5 B-tree v2).
+fn chunk_index_type(path: &std::path::Path, name: &str) -> u8 {
+    let data = std::fs::read(path).unwrap();
+    let sb = Superblock::parse(&data, 0).unwrap();
+    let addr = clawhdf5_format::group_v2::resolve_path_any(&data, &sb, name).unwrap();
+    let hdr = ObjectHeader::parse(&data, addr as usize, sb.offset_size, sb.length_size).unwrap();
+    let msg = hdr
+        .messages
+        .iter()
+        .find(|m| m.msg_type == MessageType::DataLayout)
+        .expect("layout message");
+    match DataLayout::parse(&msg.data, sb.offset_size, sb.length_size).unwrap() {
+        DataLayout::Chunked {
+            version: 4,
+            chunk_index_type: Some(t),
+            ..
+        } => t,
+        other => panic!("{name}: expected a v4 chunked layout, got {other:?}"),
+    }
+}
+
+struct Case {
+    name: &'static str,
+    /// Python keyword arguments to `create_dataset` besides `data`.
+    kwargs: &'static str,
+    index_type: u8,
+}
+
+const SHAPE: [u64; 2] = [37, 23];
+
+const CASES: &[Case] = &[
+    Case {
+        name: "implicit",
+        // Early allocation, no filters, fixed maximum: the implicit index.
+        kwargs: "chunks=(5, 4), dcpl=early()",
+        index_type: 2,
+    },
+    Case {
+        name: "fixed_array",
+        kwargs: "chunks=(5, 4), compression='gzip'",
+        index_type: 3,
+    },
+    Case {
+        name: "extensible_array",
+        kwargs: "chunks=(5, 4), maxshape=(None, 23), compression='gzip'",
+        index_type: 4,
+    },
+    Case {
+        name: "btree2",
+        kwargs: "chunks=(5, 4), maxshape=(None, None), compression='gzip'",
+        index_type: 5,
+    },
+    Case {
+        name: "single_chunk",
+        kwargs: "chunks=(37, 23), compression='gzip'",
+        index_type: 1,
+    },
+    Case {
+        name: "single_chunk_unfiltered",
+        kwargs: "chunks=(37, 23)",
+        index_type: 1,
+    },
+];
+
+/// `(start, stride, count, block)` per dimension; the first few stay below
+/// half the dataset (bounding-box path), the rest exceed it (full path).
+/// `(start, stride, count, block)` of a 2-D hyperslab.
+type Hyperslab2 = ([u64; 2], [u64; 2], [u64; 2], [u64; 2]);
+
+fn selections() -> Vec<Hyperslab2> {
+    vec![
+        ([0, 0], [1, 1], [3, 23], [1, 1]),  // ds[0:3]
+        ([7, 3], [1, 1], [9, 6], [1, 1]),   // interior window across chunks
+        ([36, 22], [1, 1], [1, 1], [1, 1]), // last element (edge chunk)
+        ([2, 1], [3, 4], [4, 3], [1, 1]),   // strided, small
+        ([0, 0], [1, 1], [30, 23], [1, 1]), // most rows
+        ([1, 0], [2, 1], [18, 23], [1, 1]), // every other row, spanning all
+        ([0, 2], [1, 1], [37, 20], [1, 1]), // columns 2..22 of every row
+        ([3, 1], [5, 3], [7, 7], [2, 2]),   // strided blocks over everything
+    ]
+}
+
+fn py_slice(start: u64, stride: u64, count: u64, block: u64) -> String {
+    // Each case is expressible as a numpy index when block == 1; with a block
+    // the selected indices are listed explicitly.
+    let idx: Vec<String> = (0..count)
+        .flat_map(|c| (0..block).map(move |b| start + c * stride + b))
+        .map(|i| i.to_string())
+        .collect();
+    format!("[{}]", idx.join(","))
+}
+
+#[test]
+fn partial_hyperslabs_of_every_v4_chunk_index_match_h5py() {
+    if !python_available() {
+        assert!(
+            !interop_required(),
+            "CLAWHDF5_REQUIRE_INTEROP=1 but python3 with h5py is not available"
+        );
+        eprintln!("SKIP: python3 with h5py not available");
+        return;
+    }
+    let dir = tempfile::tempdir().unwrap();
+    let path = dir.path().join("v4_index_selection.h5");
+    let path_str = path.display().to_string();
+
+    // h5py writes the file, then reads every selection back and prints the
+    // values, one line per (case, selection).
+    let mut script = format!(
+        "import h5py, numpy as np\n\
+         def early():\n\
+         \x20   p = h5py.h5p.create(h5py.h5p.DATASET_CREATE)\n\
+         \x20   p.set_alloc_time(h5py.h5d.ALLOC_TIME_EARLY)\n\
+         \x20   return p\n\
+         data = (np.arange({n}, dtype='<i4') * 7 - 1000).reshape({r}, {c})\n\
+         with h5py.File(r'{path_str}', 'w', libver='latest') as f:\n",
+        n = SHAPE[0] * SHAPE[1],
+        r = SHAPE[0],
+        c = SHAPE[1],
+    );
+    for case in CASES {
+        script += &format!(
+            "    f.create_dataset('{}', data=data, {})\n",
+            case.name, case.kwargs
+        );
+    }
+    script += &format!("with h5py.File(r'{path_str}', 'r') as f:\n");
+    for case in CASES {
+        for (start, stride, count, block) in selections() {
+            let rows = py_slice(start[0], stride[0], count[0], block[0]);
+            let cols = py_slice(start[1], stride[1], count[1], block[1]);
+            script += &format!(
+                "    print(' '.join(map(str, f['{}'][{rows}][:, {cols}].ravel())))\n",
+                case.name
+            );
+        }
+    }
+    let out = run_python(&script);
+    let mut expected = out.lines();
+
+    let file = File::open(&path).unwrap();
+    for case in CASES {
+        assert_eq!(
+            chunk_index_type(&path, case.name),
+            case.index_type,
+            "{}: h5py did not produce the intended chunk index",
+            case.name
+        );
+        let ds = file.dataset(case.name).unwrap();
+        assert_eq!(ds.shape().unwrap(), SHAPE);
+        for (start, stride, count, block) in selections() {
+            let sel = Selection::Hyperslab {
+                start: start.to_vec(),
+                stride: stride.to_vec(),
+                count: count.to_vec(),
+                block: block.to_vec(),
+            };
+            let want: Vec<i32> = expected
+                .next()
+                .expect("h5py printed too few lines")
+                .split_whitespace()
+                .map(|v| v.parse().unwrap())
+                .collect();
+            let raw = ds
+                .read_selection(&sel)
+                .unwrap_or_else(|e| panic!("{}: read_selection {sel:?} failed: {e}", case.name));
+            let got: Vec<i32> = raw
+                .as_chunks::<4>()
+                .0
+                .iter()
+                .map(|b| i32::from_le_bytes(*b))
+                .collect();
+            assert_eq!(got, want, "{}: selection {sel:?}", case.name);
+            assert_eq!(
+                ds.read_i32_selection(&sel).unwrap(),
+                want,
+                "{}: read_i32_selection {sel:?}",
+                case.name
+            );
+        }
+    }
+    assert!(expected.next().is_none(), "h5py printed extra lines");
+}

docs/QUICKSTART.md

@@ -395,19 +395,22 @@ clawhdf5 --path agent.h5 snapshot backup_2026-03-19.h5
 Read HDF5 files from Python without libhdf5:
 
 ```bash
-pip install clawhdf5  # coming soon — build from source for now
-cd crates/clawhdf5-py && maturin develop
+# Not on PyPI yet: build from source into a virtualenv
+pip install maturin numpy
+cd crates/clawhdf5-py && maturin develop --release
 ```
 
 ```python
 import clawhdf5
 
-# Read
-f = clawhdf5.open("data.h5")
-temps = f.read_f64("temperatures")
-print(temps)  # [22.5, 23.1, 21.8]
+# Read (h5py-style)
+with clawhdf5.File("data.h5", "r") as f:
+    temps = f["temperatures"][:]
+    print(temps)  # [22.5 23.1 21.8]
 ```
 
+See `crates/clawhdf5-py/README.md` for the supported types and indexing.
+
 ---
 
 ## Common Patterns

docs/known-issues.md

@@ -7,6 +7,23 @@ deleting it.
 
 ---
 
+## Selection reads that decode more than the selection
+
+**Status:** open (documented 2026-09-26). `Dataset::read_selection` (and so
+the Python `ds[...]`) materialises only the selection's bounding box when
+that box covers at most half the dataset (`partial_read`). It decodes the
+whole dataset and extracts the selection instead when:
+- the bounding box covers more than half the dataset — which a strided
+  selection across a chunked dataset (`ds[::100]`) always does, although
+  it may touch few chunks;
+- the dataset is compact or virtual, or has no storage;
+- it is chunked with a non-default fill value (the box path does not fill
+  unallocated chunks, so the fill-aware full read is used).
+Values are correct in every case; this is cost only. Selections other than
+`Selection::All` also bypass the file's chunk cache. The bounding-box
+heuristic's other cost is measured under "Concurrent and contiguous read
+performance" below.
+
 ## Concurrent and contiguous read performance (measured 2026-09-26)
 
 **Status:** open for chunked full reads (one cause fixed 2026-09-26); the

scripts/ci-test.sh

@@ -7,7 +7,9 @@
 #
 # Environment:
 #   CLAWHDF5_REQUIRE_INTEROP=1  Fail (instead of skip) when python3 with
-#                               h5py/netCDF4/xarray is missing. CI sets this.
+#                               h5py/netCDF4/xarray is missing, or without
+#                               maturin/pytest for the Python package step.
+#                               CI sets this.
 #                               Unset locally, the interop steps are skipped
 #                               if python3+h5py is not importable.
 #   CLAWHDF5_FUZZ_SECONDS=N     Run each cargo-fuzz target for N seconds
@@ -55,7 +57,8 @@ run_step "cargo fmt --check" cargo fmt --check
 
 # 2. Clippy over every target (lib, bins, tests, benches, examples). Without
 #    --all-targets, test and bench code is never linted. clawhdf5-py is
-#    excluded because it needs PyO3/Python headers.
+#    excluded here: PyO3 needs a Python interpreter to build, so it is
+#    linted in the Python package step (5b) instead.
 run_step "cargo clippy --all-targets" cargo clippy \
     --workspace \
     --exclude clawhdf5-py \
@@ -216,6 +219,36 @@ else
     STEPS+=("SKIP: h5py interop (format, ignored tests)")
 fi
 
+# 5b. The Python package (crates/clawhdf5-py): lint it, build the wheel with
+#     maturin and run its pytest suite, which compares every read with h5py.
+#     The wheel is unpacked under target/ and put on PYTHONPATH, so the
+#     interpreter's environment is left as it was. Needs maturin and pytest
+#     in $PYTHON (CI installs both); skipped without them, and a failure
+#     instead when CLAWHDF5_REQUIRE_INTEROP=1.
+python_package() {
+    local root="$SCRIPT_DIR/.." out
+    out="${CARGO_TARGET_DIR:-$root/target}/py-package"
+    rm -rf "$out" && mkdir -p "$out/wheel" "$out/site" || return 1
+    cargo clippy -p clawhdf5-py --all-targets -- -D warnings || return 1
+    "$PYTHON" -m maturin build \
+        -m "$root/crates/clawhdf5-py/Cargo.toml" \
+        -i "$PYTHON" \
+        --out "$out/wheel" || return 1
+    "$PYTHON" -m pip install --quiet --no-deps --target "$out/site" "$out"/wheel/*.whl || return 1
+    PYTHONPATH="$out/site" "$PYTHON" -m pytest -q -p no:cacheprovider \
+        "$root/crates/clawhdf5-py/tests"
+}
+if "$PYTHON" -m maturin --version >/dev/null 2>&1 && "$PYTHON" -c "import pytest" >/dev/null 2>&1; then
+    run_step "Python package (maturin build + pytest vs h5py)" python_package
+elif [ "${CLAWHDF5_REQUIRE_INTEROP:-0}" = "1" ]; then
+    run_step "Python package (maturin build + pytest vs h5py)" \
+        bash -c "echo \"maturin and pytest are required in $PYTHON (CLAWHDF5_REQUIRE_INTEROP=1)\"; exit 1"
+else
+    echo ""
+    echo "==> [Python package] SKIPPED: needs maturin and pytest in $PYTHON"
+    STEPS+=("SKIP: Python package (maturin build + pytest)")
+fi
+
 # 6. Benches must keep compiling (they are not run).
 run_step "cargo bench --no-run" cargo bench \
     --workspace \