Interacting with HDF5 files
Author: Philipp Schreiner
Created: Oct. 16, 2024
Last updated: Jul. 2, 2026 by Philipp Schreiner
Now that we created our first DataHandler object, it is as good a time as ever to talk about it in more detail. cait uses HDF5 files to store data and intermediate results for later use. This can be voltage traces of triggered events, timestamps of testpulses, calculcated recoil energies and many more. The DataHandler class makes interacting with these files very convenient and safe. You should never be in need of opening the HDF5 file yourself – a practice that we heavily discourage as it can lead to corrupted files!
In the following we explore the core features of the DataHandler class.
import os
import numpy as np
import cait as ai
import cait.versatile as vai
Mock data from previous tutorial
First, we have to have a DataHandler before we can play with it. We will use the one created in the Triggering stream data tutorial.
fdirh5 = "tutorial_output"
hdf5_name = "my_first_trigger"
os.makedirs(fdirh5, exist_ok=True)
record_length = 2**14
trigger_config = {
"trigger_channels": ["Ch0"],
"passive_channels": ["Ch1"],
"testpulse_channels": ["TP0", "TP1"],
"controlpulses_above": [9., 9.],
"f_noise": 1000,
"copy_events": True,
}
n_channels = len(trigger_config["trigger_channels"]) + len(trigger_config["passive_channels"])
stream = vai.MockStream(seed=137, rate_Hz=2)
If you haven’t run the previous notebook, the following cell will bring you up to where we left off:
dh = ai.DataHandler(record_length=record_length,
nmbr_channels=n_channels,
sample_frequency=stream.sample_frequency)
dh.set_filepath(path_h5=fdirh5, fname=hdf5_name, appendix=False)
dh.init_empty()
dh.trigger_zscore(stream, **trigger_config)
for group in ["events", "testpulses", "noise", "controlpulses"]:
dh.cmp(group)
If you already have the file ready, you can instantiate the DataHandler like so:
# Create DataHandler instance
dh = ai.DataHandler(record_length=record_length,
nmbr_channels=n_channels,
sample_frequency=stream.sample_frequency)
# Set the path to the desired HDF5 file
dh.set_filepath(path_h5="tutorial_output", fname="my_first_trigger", appendix=False)
DataHandler Instance created.
Inspecting File
A quick way to get a summary of a DataHandler’s properties is to print it.
print(dh)
DataHandler linked to HDF5 file 'tutorial_output/my_first_trigger.h5'
HDF5 file size on disk: 882.7 MiB
Groups in file: ['_ext_stored_data', 'controlpulses', 'event_building-z-score', 'events', 'noise', 'testpulses', 'triggers-z-score'].
The HDF5 file contains references to externally stored events for the following groups: ['events', 'testpulses', 'controlpulses', 'noise']
Time between first and last testpulse: 1.00 h
First testpulse on/at: 2015-03-14 08:26:58+00:00 (UTC)
Last testpulse on/at: 2015-03-14 09:26:48+00:00 (UTC)
Depending on the available groups/datasets in the HDF5 file, this summary is more or less informative (e.g. the testpulse stats can of course only be printed if available in the file).
The filename/path/directory of the HDF5 file can be accessed via the following methods which support both relative and absolute paths. We highly recommend to use those methods when referencing the HDF5 file and to not directly use properties of the DataHandler.
dh.get_filepath() # tutorial_output/my_first_trigger.h5
dh.get_filename() # my_first_trigger
dh.get_filedirectory() # tutorial_output
dh.get_filedirectory(absolute=True) # /Users/.../tutorial_output
We can get a more detailed view into the contents of the HDF5 file with the content-method which takes the group of interest as an argument (you can also use wildcards like "events*" to list all groups with names starting with "events"). If no group is specified, the datasets of all groups are listed.
This method can be used to inspect the datasets’ shapes and data types.
If you want to get an explanation on how to interpret the output, set print_info=True.
dh.content("events")
events
baseline_difference (2, 5318) float64
baseline_offset (2, 5318) float64
baseline_slope (2, 5318) float64
decay_time (2, 5318) float64
decay_time_CAT (2, 5318) float64
event (2, 5318, 16384) float32
hours (5318,) float64
integral (2, 5318) float64
max_deriv (2, 5318) float64
max_deriv_index (2, 5318) int64
maximum (2, 5318) float64
min_deriv (2, 5318) float64
min_deriv_index (2, 5318) int64
onset (2, 5318) float64
onset_CAT (2, 5318) float64
peak_position (2, 5318) float64
pulse_height (2, 5318) float64
rise_time (2, 5318) float64
rise_time_CAT (2, 5318) float64
rms (2, 5318) float64
time_mus (5318,) int32
time_s (5318,) int32
|timestamps (5318,)
trigger_flag (2, 5318) bool
trigger_phs (2, 5318) float32
trigger_timestamps (2, 5318) int64
variance (2, 5318) float64
Access Data in the HDF5
The groups and datasets in the HDF5 file underlying the DataHandler are accessed through the get method. To read out the full content of a dataset you specify the group and the dataset name.
Alternatively, you can also slice the dataset already in the call to get. E.g. you can access only the first channel of the data, or provide a boolean flag (or a list of indices). The get method supports up to 3 indices corresponding to the (up to) 3 dimensions of the data, and they are understood in increasing order. Therefore, to index only the first and third dimension for example, you have to provide None for the second dimension.
Notice that due to a h5py limitation, only one dimension at a time can be indexed with boolean flags/index lists.
# All pulse heights of all channel; shape: (2, 5318)
dh.get("events", "pulse_height")
# Voltage trace of the 37th event in first channel; shape: (1, 1, 16384)
dh.get("events", "event", 0, 37)
# Voltage traces of all events in first channel with pulse height < 1; shape: (1, 5318, 16384)
flag = dh.get("events", "pulse_height", 0) < 1
dh.get("events", "event", 0, flag)
# Pulse heights for all channels where pulse height of first channel <1; shape: (2, 5318)
dh.get("events", "pulse_height", None, flag)
Note that the following two calls to get yield the same values yet the second one is more memory-efficient because the slicing is done on the HDF5 file. The first option loads the entire dataset into memory and performs the slicing afterwards.
dh.get("events", "event")[0, flag] # slicing of numpy.ndarray
dh.get("events", "event", 0, flag) # slicing of HDF5 file
For seasoned cait-users, the following short-cut notation could be a huge time save. It is equivalent to using the dh.get() method but it allows for more concise code:
dh["events/pulse_height"] # equivalent to dh.get("events", "pulse_height")
dh["events/event", 0, 37] # equivalent to dh.get("events", "event", 0, 37)
dh["events/event", 0, flag] # equivalent to dh.get("events", "event", 0, flag)
dh["events/pulse_height", :, flag] # equivalent to dh.get("events", "pulse_height", None, flag)
I.e. you can slice the DataHandler object with keys of the form "<group>/<dataset>" and possibly additional indices. Notice, that here you can actually use the colon (:) operator, which in the get method had to be replaced by None.
Furthermore, this notation allows for TAB-auto-completion in iPython/notebook contexts (type dh[ and press TAB).
Something that you will often want to do after accessing data is to plot it. For that, vai.Histogram exists (look into the cait.versatile plotting class documentation for more plotting types):
vai.Histogram(dh["events/pulse_height", 1], bins=(0, 0.02, 200), xlabel="Pulse height second channel (V)", backend="plotly");
Accessing event data
Please, never (!) access and loop over events manually, e.g. things like for i in range(1000): do_something(dh["events/event", 0, i]). Since this is a reoccuring task, a special mechanism exists in cait, which is the get_event_iterator method. It lets you load events conveniently and most importantly efficiently!
# Create an event iterator to group "events"
events_it = dh.get_event_iterator("events")
# You can conveniently have a look at the events and click through them:
vai.Preview(events_it, backend="plotly");
If you want to perform some calculation on those events, e.g. apply some custom function to each event, use vai.apply:
# Function can be arbitrarily complex. Here just a variance calculation
# for demonstration.
def my_custom_function(event):
return np.var(event, axis=-1)
# Calculate the output for all events and both channels.
# The output has shape (5318, 2), events are stacked as rows.
all_vars = vai.apply(
my_custom_function,
events_it,
n_processes=1, # Only needed for automated test (you can remove this argument)
)
Write Data to the HDF5 File
To write data to the HDF5 file, you use the set method for which there are two main use cases:
writing an array into the file as is (i.e. we keep the shape and dtype the same)
creating a multi-channel array and writing only to one channel (e.g. you calculate some property for the phonon and light channel separately but want to write them to the same HDF5 dataset)
The two use cases are illustrated below.
Notice that you can write multiple datasets at once as long as they are in the same group. Dataset names and their content are parsed as keyword arguments.
If you want the datasets to have a specific dtype, you can hand it to set as well.
data1, data2 = np.random.rand(1, 37, 100), np.random.rand(1, 37 ,100)
# Include 'data1' and 'data2' as datasets 'new_ds1' and 'new_ds2' in group 'new_group'
# ('new_ds1' and 'new_ds2' do not yet exist)
dh.set(group="new_group", new_ds1=data1, new_ds2=data2)
# Include 'data1' and 'data2' as datasets 'ds1_mc' and 'ds2_mc' in group 'new_group2'
# ('data1' and 'data2' are 1-dimensional but we want to create 2-dimensional
# datasets (for different channels e.g.) and write the data into the 0-th channel. This also
# works for writing single channels to already existing multi-channel datasets.)
dh.set(group="new_group2", n_channels=2, channel=0, ds1_mc=data1[0], ds2_mc=data2[0])
If the dataset(s) already exist in the given group, you have to explicitly allow set to change them:
# Include 'data1' and 'data2' as datasets 'ds1' and 'ds2' in group 'new_group'
# (either or both of 'ds1' and 'ds2' already exist and have correct shape/dtype for new
# data)
dh.set(group="new_group", ds1=data1, ds2=data2, change_existing=True)
# Include 'data1' and 'data2' as datasets 'ds1' and 'ds2' in group 'new_group'
# (either or both of 'ds1' and 'ds2' already exist and have incorrect shape/dtype for new
# data, but we want to force the new dtype/shape)
dh.set(group="new_group", ds1=data1, ds2=data2, overwrite_existing=True)
By calling content again, we can look at the newly created datasets.
Renaming and deleting datasets/groups, repackaging
You can rename and delete datasets/groups using the methods rename and drop as follows. Just like with set, you can do this for multiple datasets simultaneously as long as they are within the same group. The syntax once more uses keyword arguments.
Note
Due to the possibility of changing both the group name and dataset names with a single method, renaming groups called group is obviously not possible.
# Rename groups 'new_group' and 'new_group2' to 'new_group_renamed' and 'new_group2_renamed'
dh.rename(new_group="new_group_renamed", new_group2="new_group2_renamed")
# Rename datasets 'ds1' and 'ds2' in group 'new_group_renamed' to 'new_ds1' and 'new_ds2'
dh.rename(group="new_group_renamed", ds1="ds1_renamed", ds2="ds2_renamed")
Groups/datasets are deleted as follows:
# Drop group 'new_group_renamed'
dh.drop("new_group_renamed")
# Drop dataset 'ds1' from group 'new_group_renamed'
dh.drop("new_group2_renamed", "ds1_mc")
Important
Deleting groups/datasets does not decrease the HDF5 file’s size due to the HDF5 file’s data structure. To decrease the size, repackaging has to be performed for which we provide the method repackage. Alternatively, this method can automatically be called when deleting a dataset/group using the repackage=True argument of the drop method.
dh.repackage()
Successfully repackaged 'tutorial_output/my_first_trigger.h5'. Memory saved: 557.2 KiB
Combining multiple HDF5 files
The HDF5 library allows to link separate files into a single file without copying the data. The resulting file looks identical to a regular HDF5 file. To avoid confusion in such cases, a number of quality-of-life features are included in DataHandler.
The need for virtual datasets arises for example when analysing more than one file at a time. Below we show how files can be combined and how to interact with them.
At first, we have to have multiple files. For that, we show an elegant way to load multiple files in dictionaries. In a real analysis, the vai.MockStream will be replaced by an actual Stream object as described in the Triggering stream data tutorial.
fnames = [
"my_file_to_merge_0",
"my_file_to_merge_1",
]
# Create a dictionary whose keys are the file names
# and whose values are stream objects
streams = {
fname: vai.MockStream(rate_Hz=2, duration_h=0.5)
for fname in fnames
}
# Create a DataHandler object for each stream, also in
# a dictionary
dhs = {
fname: ai.DataHandler(
record_length=record_length,
nmbr_channels=n_channels,
sample_frequency=stream.sample_frequency,
)
for fname, stream in streams.items()
}
# Initialize empty HDF5 files
for k, v in dhs.items():
v.set_filepath(path_h5=fdirh5, fname=k, appendix=False)
v.init_empty()
# Trigger all files and calculate main parameters
for stream, dh in zip(streams.values(), dhs.values()):
dh.trigger_zscore(stream, **trigger_config)
for group in ["events", "testpulses", "noise", "controlpulses"]:
dh.cmp(group)
Now we have two files, 'my_file_to_merge_0.h5' and 'my_file_to_merge_1.h5', and two corresponding DataHandler instances (collected in the dictionary). It would be much more convenient to have one single DataHandler which combines all data for the combined analysis. This is achieved with the combine_h5 method:
ai.data.combine_h5(
fname="my_first_merged_file", # name of the merged file
files=fnames, # files to be merged
src_dir=fdirh5, # location of the source files
out_dir=fdirh5, # location where to put the merged file
groups_combine=["testpulses", "noise"], # which of the DataHandler groups to combine
)
Successfully combined files ['my_file_to_merge_0', 'my_file_to_merge_1'] into 'tutorial_output/my_first_merged_file.h5' (56.0 KiB).
Calculating extended hours for all groups with datasets event, hours, time_s, time_mus:
DataHandler Instance created.
Successfully written hours with shape (1000,) and dtype 'float32' to group noise.
Successfully written hours with shape (576,) and dtype 'float32' to group testpulses.
As you can see, the resulting file size is only a few KiB due to no data being copied. We can now create a DataHandler to the combined file and inspect it with our usual methods:
dh_combined = ai.DataHandler(
record_length=record_length,
nmbr_channels=n_channels,
sample_frequency=stream.sample_frequency,
)
dh_combined.set_filepath(path_h5=fdirh5, fname="my_first_merged_file", appendix=False)
DataHandler Instance created.
print(dh_combined)
DataHandler linked to HDF5 file 'tutorial_output/my_first_merged_file.h5'
HDF5 file size on disk: 64.1 KiB
Groups in file: ['_ext_stored_data', 'noise', 'testpulses'].
The HDF5 file contains virtual datasets linked to the following files: {'tutorial_output/my_file_to_merge_1.h5', 'tutorial_output/my_file_to_merge_0.h5'}
All of the external sources are currently available.
The HDF5 file contains references to externally stored events for the following groups: ['testpulses', 'noise']
Time between first and last testpulse: 0.50 h
First testpulse on/at: 2015-03-14 08:26:58+00:00 (UTC)
Last testpulse on/at: 2015-03-14 08:56:48+00:00 (UTC)
We see that print(dh_combined) tells us which files are linked to the HDF5 file of our DataHandler object. In principle it could happend that some of these source files get deleted. In this case, you could still create the DataHandler but subsequent calculations can have unexpected behaviour. This is why print(dh_combined) also tells you whether all sources are available or not.
Using content now also changed slightly as it shows you an indicator (v) next to datasets which are virtual:
dh_combined.content("testpulses")
testpulses
baseline_difference (v) (2, 576) float64
baseline_offset (v) (2, 576) float64
baseline_slope (v) (2, 576) float64
decay_time (v) (2, 576) float64
decay_time_CAT (v) (2, 576) float64
event (v) (2, 576, 16384) float32
hours (576,) float32
integral (v) (2, 576) float64
max_deriv (v) (2, 576) float64
max_deriv_index (v) (2, 576) int64
maximum (v) (2, 576) float64
min_deriv (v) (2, 576) float64
min_deriv_index (v) (2, 576) int64
onset (v) (2, 576) float64
onset_CAT (v) (2, 576) float64
peak_position (v) (2, 576) float64
pulse_height (v) (2, 576) float64
rise_time (v) (2, 576) float64
rise_time_CAT (v) (2, 576) float64
rms (v) (2, 576) float64
testpulseamplitude (v) (2, 576) float32
time_mus (v) (576,) int32
time_s (v) (576,) int32
|timestamps (576,)
variance (v) (2, 576) float64
In principle, this is all the difference there is. You can do analysis with dh_combined in the same way you would do with dh. The only thing worth mentioning is that writing to virtual datasets is special. Depending on your use case you could either want to change the data in all source files (which could lead to confusion and unexpected results), or you could want to detach the virtual dataset and write to a regular dataset of that name. In most cases, the latter will probably be the desired behaviour. As a safety feature, the set method checks whether you are attempting to write to a virtual dataset or not and you have to make this decision upon calling it by setting write_to_virtual=False or True.
To actually merge the files, use cait.data.merge_h5 instead of cait.data.combine_h5. This will create a new file containing all the data, i.e. you could delete the individual files afterwards. We still recommend to use cait.data.combine_h5.