The slides should be fine if the zoom is 200%. You will need at least this level, as the screen in MD 1-21 is relatively small and difficult to see from the back.

# Data IO (input/output)¶

## Introduction¶

ESRF data (used to) come in (too many) different formats:

- Specfile, EDF, HDF5
- and specific detector formats: MarCCD, Pilatus CBF, Dectris Eiger, …

HDF5 is now the standard ESRF data format, so we will only focus on it today.

Methods for accessing other file formats are described in the io_spec_edf.ipynb notebook.

# HDF5¶

## What is HDF5?¶

HDF5 (for Hierarchical Data Format version 5) is a file format to structure and store complex and large volumes of data.

## Why HDF5?¶

- high-performance (binary)
- portable file format (Windows, Linux, MacOS)
- self-describing extensible types, rich metadata
- support data compression
- free (& open source)
- widely used in scientific computing
- adopted by many institutes (NASA, LIGO, ...)
- adopted by most of the synchrotrons (ESRF, SOLEIL, Desy ...)

## HDF5 data structure overview¶

HDF5 organizes data in a hierarchical structure, similar to a file system. Data can be almost anything: images, tables, graphs, and documents.

The HDF5 container consists of the following main components:

- File: the root of the container, which serves as the top-level node.
- Group: a hierarchical grouping structure containing other groups or datasets that organize and categorize data. Equivalent of a folder.
- Dataset: a multidimensional array of data elements, representing a collection of values stored in memory. Equivalent of a data file (.txt, .dat, .csv ...)

## HDF5 hierarchy¶

## An HDF5 example¶

Here is an example of the file generated by pyFAI.

#### Additionally, HDF5 features include:¶

- Links: references to other files or datasets within the container.
- Attributes: metadata associated with a dataset or group.
- Datatypes: definitions for the structure and storage format of data (integer, string...)
- Virtual datasets: interfaces that provide access to data without storing it in memory.

## Useful tools for inspecting HDF5 files¶

### HDFGroup tools for viewing and editing HDF5¶

Command line and desktop application: `h5ls`, `h5dump`, `hdfview`

    >>> h5ls -r my_first_one.h5
        /                        Group
        /data1                   Dataset {100, 100}
        /group1                  Group
        /group1/data2            Dataset {100, 100}

### `h5glance`¶

Jupyter notebook and command line tool for browsing HDF5 files.

In [ ]:

    # From the Jupyter notebook.
    from h5glance import H5Glance

    H5Glance("data/water.h5")

In [ ]:

    # From the command line.
    !h5glance data/water.h5

### `jupyterlab-h5web`¶

JupyterLab extension to explore and visualize HDF5 file contents.

In [ ]:

    from jupyterlab_h5web import H5Web

    H5Web("data/water.h5")

### `silx view`¶

Desktop application file browser/viewer

    pip install silx
    silx view my_file.h5

In [ ]:

    !silx view data/water.h5

# h5py¶

## What is h5py?¶

h5py is a Python binding for accessing HDF5 files.

- It allows the reading/writing of HDF5 files using a simple Pythonic interface.

- For a python user, `h5py` is the only package to learn and use

### What is h5py?¶

Originally written by Andrew Collette.

## How to install h5py?¶

    pip install h5py

In [ ]:

    import h5py

    print("h5py:", h5py.version.version)
    print("hdf5:", h5py.version.hdf5_version)

## Opening and creating HDF5 files¶

### `h5py.File`¶

- First open the file with h5py.File:

      h5py.File('myfile.hdf5', mode)

- Opening modes:

Mode | Meaning
--- | ---
r | Read-only, file must exist; default with h5py v3
r+ | Read/write, file must exist
w | Create file, truncate if exists
w- or x | Create file, fail if exists
a | Read/write if exists, create otherwise; default with h5py v2

### Opening an existing file¶

In [ ]:

    h5file = h5py.File("data/water.h5", mode="r")
    ...
    h5file.close()

### Creating a new file¶

In [ ]:

    h5file = h5py.File("data/new_data.h5", mode="w")
    ...
    h5file.close()

### Using a context manager¶

- Context managers guarantee that resources are released. In our case, it ensures that the HDF5 file is closed.
- Usually used from the `with` statement.

To safely access a HDF5 file, do:

In [ ]:

    # Instead of
    h5file = h5py.File("data/water.h5", mode="r")
    h5file.close()

    # Use
    with h5py.File("data/water.h5", mode="r") as h5file:
        print(h5file)

### File structure¶

#### `h5py.Group`¶

Documentation: Group

##### Creating a group¶

Documentation: Group.create_group:

In [ ]:

    # List 'entry_0000' group children.
    with h5py.File("data/new_data.h5", mode="w") as h5file:
        my_h5group = h5file.create_group("my_h5group")
        my_h5subgroup = my_h5group.create_group("my_h5subgroup")

In [ ]:

    # From the Jupyter notebook.
    from h5glance import H5Glance

    H5Glance("data/new_data.h5")

##### Browsing groups¶

The file content can be accessed using a dictionary-like API, using methods like Group.keys(), Group.values(), or Group.items().

In [ ]:

    # Print the available group names at the root level.
    with h5py.File("data/water.h5", mode="r") as h5file:
        print(list(h5file.keys()))

In [ ]:

    from pprint import pprint

    # List the children of the 'entry_0000' group.
    with h5py.File("data/water.h5", mode="r") as h5file:
        h5group = h5file["entry_0000"]
        pprint(dict(h5group.items()))

#### `h5py.Dataset`¶

Documentation: Dataset

##### Reading a dataset¶

This is probably the most common task when dealing with HDF5 files.

Note that datasets have properties similar to those of the NumPy arrays.

In [ ]:

    with h5py.File("data/water.h5", mode="r") as h5file:
        h5dataset = h5file["/entry_0000/4_azimuthal_integration/results/I"]
        print("h5dataset:", h5dataset.shape, h5dataset.dtype, h5dataset.size)

#### You can copy the data to a NumPy array using one of the following ways.¶

In [ ]:

    import numpy as np

    with h5py.File("data/water.h5", mode="r") as h5file:
        h5dataset = h5file["/entry_0000/4_azimuthal_integration/results/I"]
        data = h5dataset[()] # Recommended. Copy the dataset to a NumPy array.
        data = h5dataset[...] # Same as above.
        data = np.array(h5dataset) # Another way.
        data = h5dataset[:] # This also works, but avoid it as it fails in some cases.
        print("data type:", type(data), "; shape", data.shape, "; min.:", data.min())

Note that once the file is closed, the `h5dataset` object will be unavailable.

In [ ]:

    with h5py.File("data/water.h5", mode="r") as h5file:
        h5dataset = h5file["/entry_0000/4_azimuthal_integration/results/I"]
        data = h5dataset[()]  # Copy the whole dataset to a NumPy array.
    print(h5dataset) # Not available anymore, but the copied data is.
    print(data[:5]) # Print the first five elements.

Not very convenient for interactive browsing... this is why silx view, h5web, h5glance ... exists.

##### Writing a dataset¶

h5py will create all missing groups to solve the dataset location. So you usually won't have to call create_group.

In [ ]:

    with h5py.File("data/new_data.h5", mode="a") as h5file:
        h5file["my_h5group/to/my_h5dataset"] = np.random.rand(100, 100)

Alternative: using Group.create_dataset

In [ ]:

    with h5py.File("data/new_data.h5", mode="w") as h5file:
        h5file.create_dataset("my_data", data=np.arange(100))

In [ ]:

    H5Glance("data/new_data.h5")

#### Attributes¶

Attributes are the way to store metadata for a group or a dataset.

Group and dataset have a small `'<obj>.attrs'` attached.

Warning attributes must be of a limited size

##### Writing an attribute¶

In [ ]:

    with h5py.File("data/new_data.h5", "w") as h5file:
        dataset = h5file.create_dataset("my_dataset", data=np.random.rand(10, 10))
        dataset.attrs["description"] = "This is a random dataset"

##### Reading an attribute¶

In [ ]:

    with h5py.File("data/new_data.h5", "r") as h5file:
        print(h5file["my_dataset"].attrs["description"])

## Exercise: Flat field correction¶

Flat-field correction is a technique used to improve quality in digital imaging.

The goal is to normalize images and remove artifacts caused by variations in the pixel-to-pixel sensitivity of the detector and/or by distortions in the optical path. (see https://en.wikipedia.org/wiki/Flat-field_correction)

$$ normalized = \frac{raw - dark}{flat - dark} $$

- `normalized`: Image after flat field correction
- `raw`: Raw image. It is acquired with the sample.
- `flat`: Flat field image. It is the response given out by the detector for a uniform input signal. This image is acquired without the sample.
- `dark`: Also named `background` or `dark current`. It is the response given out by the detector when there is no signal. This image is acquired without the beam.

Here is a function implementing the flat field correction:

Note: make sure you execute the cell for defining this function

In [ ]:

    import numpy as np

    def flatfield_correction(raw, flat, dark):
        """
        Apply a flat-field correction to a raw data using a flat and a dark.
        """
        # Make sure that the computation is done using float
        # to avoid type overflow or loss of precision
        raw = raw.astype(np.float32)
        flat = flat.astype(np.float32)
        dark = dark.astype(np.float32)
        # Do the computation
        return (raw - dark) / (flat - dark)

If you like to plot an image, you can use `matplotlib`'s `imshow` function.

The `%matplotlib` "magic" command should be called once first.

In [ ]:

    %matplotlib inline

    from matplotlib import pyplot as plt

In [ ]:

    plt.imshow(np.random.random((20, 60)))

### Exercise 1¶

1. Browse the file `data/ID16B_diatomee.h5`.
2. Get a single raw dataset, a flat field dataset, and a dark image dataset from this file.
3. Apply the flat field correction.
4. Save the result into a new HDF5 file.

If you are stuck, the solution is provided in the file solutions/exercise1.py

In [ ]:

    from jupyterlab_h5web import H5Web

    H5Web("data/ID16B_diatomee.h5")

In [ ]:

    import h5py

    with h5py.File("data/ID16B_diatomee.h5", mode="r") as h5s:
        pass
        # Step 1: Read the data.

        # raw_data_path = ...
        # raw_data = ...

        # flat_path = ...
        # flat = ...

        # dark_path = ...
        # dark = ...

    # Step 2: Compute the result.

    # normalized = flatfield_correction(raw_data, flat, dark)

    # Step 3: Save the result.

    # ...

### Exercise 2¶

1. Apply the flat field correction to all raw data available (use the same flat and dark for all the images).
2. Save each result into different datasets of the same HDF5 file.

If you are stuck, the solution is provided in the file solutions/exercise2.py

In [ ]:

### Exercise 3¶

From the previous exercise, we can see that the flat field correction was not very good for the last images.

Another flat field was acquired at the end of the acquisition.

We could use this information to compute a flat field closer to the image we want to normalize. This can be done with a linear interpolation of the flat images by using the name of the image as the interpolation factor (which varies between 0 and 500 in this case).

1. For each raw data, compute the corresponding flat field using linear interpolation (between `flatfield/0000` and `flatfield/0500`).
2. Save each result into different datasets in a single HDF5 file.

If you are stuck, the solution is provided in the file solutions/exercise3.py

In [ ]:

## Advanced usage¶

### Dataset compression¶

Install hdf5plugin and `import hdf5plugin`.

HDF5 provides dataset compression support. With `h5py` GZIP and LZF compression are available (see compression-filters). Yet, there are many third-party compression filters for HDF5 available.

hdf5plugin allows usage of some of those compression filters with `h5py` (Blosc, Blosc2, BitShuffle, BZip2, FciDecomp, LZ4, SZ, SZ3, Zfp, ZStd).

In [ ]:

    import h5py
    import hdf5plugin  # Allows to read dataset stored with supported compressions.

To write compressed datasets, see:

- Group.create_dataset `chunks`, `compression` and `compression_opts` parameters.
- "Chunked Storage" documentation
- hdf5plugin documentation

### Soft and external links¶

A HDF5 file can contain links to Group/Dataset:

- within the same file: see h5py.SoftLink
- in another file: see h5py.ExternalLink

Links can be dangling if the destination does not exist.

### External dataset¶

A HDF5 file can contain datasets that are stored in external binary files: See Group.create_dataset `external` parameter.

### Virtual Dataset (aka. VDS)¶

Virtual dataset allows to map multiple datasets into a single one. Once created it behaves as other datasets.

See https://docs.h5py.org/en/stable/vds.html

### HDF5 file locking: A word of caution¶

Do NOT open an HDF5 file that is currently being written (without caution).

By default, HDF5 locks the file even for reading, and other processes cannot open it for writing. This can be an issue, e.g., during acquisition.

Warning: Do not open the same file twice for writing without file locking, or the file will be corrupted.

Workarounds:

- Helper to handle HDF5 file locking: `silx.io.h5py_utils.File`.
- HDF5 file locking can be disabled by setting the `HDF5_USE_FILE_LOCKING` environment variable to `FALSE`.
- With recent version of `h5py` (>= v3.5.0): `h5py.File`'s `locking` argument.

### Chunked storage¶

By default, HDF5 datasets will be contiguous (like C). If you have specific usages or want to improve speed, you may want to define chunks. See h5py chunked storage (for advanced usage).

### Practical tools¶

- Conversion:

  - `silx convert`: To convert EDF, or spec files to HDF5

- Reading/writing HDF5 helpers:

  - `silx.io.dictdump`: `h5todict`, `dicttoh5`
  - `silx.io.utils.h5py_read_dataset`

### A word about NeXus¶

Nexus is a data format for neutron, x-ray, and muon science.

It aims to be a common data format for scientists for greater collaboration.

If you intend to store some data to be shared, it can give you a "standard way" for storing it.

The main advantage is to ensure compatibility between your data files and existing software packages (if they respect the NeXus convention) or from your software to different datasets:

- an example on how to store tomography raw data
- an example to store tomoraphy application (3D reconstruction).

## Conclusion¶

h5py provides access to HDF5 file content from Python through:

- `h5py.File` opens a HDF5 file:
  - Do not forget the mode: `'r'`, `'a'`, `'w'`.
  - Use a `with` statement, or do not forget to `close` the file.
- `h5py.Group` provides a key-value mapping `dict`-like access to the HDF5 structure.
- `h5py.Dataset` gives access to data as `numpy.ndarray`.