Implement better openpmd metadata

visualpic/data_reading/field_readers.py

@@ -55,19 +55,24 @@ def _readjust_metadata(self, field_metadata, slice_dir_i, slice_dir_j,
         geom = field_metadata['field']['geometry']
         if geom in ['cylindrical', 'thetaMode'] and theta is None:
             r_md = field_metadata['axis']['r']
+            z_md = field_metadata['axis']['z']
             # Check if resolution should be reduced
             if max_resolution_3d is not None:
-                z = field_metadata['axis']['z']['array']
+                z = z_md['array']
                 r = r_md['array']
                 nr = len(r) - 1
                 nz = len(z) - 1
                 max_res_lon, max_res_transv = max_resolution_3d
                 if nz > max_res_lon:
                     excess_z = int(np.round(nz/max_res_lon))
-                    field_metadata['axis']['z']['array'] = z[::excess_z]
+                    z_md['array'] = z[::excess_z]
+                    z_md['min'] = z_md['array'][0]
+                    z_md['max'] = z_md['array'][-1]
                 if nr > max_res_transv:
                     excess_r = int(np.round(nr/max_res_transv))
                     r_md['array'] = r[::excess_r]
+                    r_md['min'] = z_md['array'][0]
+                    r_md['max'] = z_md['array'][-1]
                 # if nr > max_res_transv:
                 #    r = zoom(r, max_res_transv/nr, order=1)
                 #    r_md['array'] = r
@@ -476,48 +481,11 @@ def _read_field_theta(
         return fld
 
     def _read_field_metadata(self, file_path, iteration, field_path):
+        # Get name of field and component.
         field, *comp = field_path.split('/')
         if len(comp) > 0:
             comp = comp[0]
-        md = {}
-        t, params = self._opmd_reader.read_openPMD_params(iteration)
-        md['time'] = {}
-        md['time']['value'] = t
-        md['time']['units'] = 's'
-        field_geometry = params['fields_metadata'][field]['geometry']
-        axis_labels = params['fields_metadata'][field]['axis_labels']
-        field_units = self._determine_field_units(field)
-        md['field'] = {}
-        md['field']['units'] = field_units
-        md['field']['geometry'] = field_geometry
-        md['field']['axis_labels'] = axis_labels
-        ax_el, ax_lims = self._opmd_reader.get_grid_parameters(
-            iteration, [field], params['fields_metadata'])
-        axes = ax_el.keys()
-        md['axis'] = {}
-        for axis in axes:
-            md['axis'][axis] = {}
-            md['axis'][axis]['units'] = 'm'
-            ax_min = ax_lims[axis][0]
-            ax_max = ax_lims[axis][1]
-            ax_els = ax_el[axis]
-            if field_geometry in ['cylindrical', 'thetaMode'] and axis == 'r':
-                ax_min = -ax_max
-                ax_els += ax_el[axis]
-            # FIXME this does not differentiate between
-            # node-centered / cell-centered fields. FieldMetaInformation
-            # does it properly
-            md['axis'][axis]['array'] = np.linspace(ax_min, ax_max, ax_els)
-        return md
-
-    def _determine_field_units(self, field):
-        if field == 'E':
-            return 'V/m'
-        elif field == 'B':
-            return 'T'
-        elif field == 'rho':
-            return 'C/m^3'
-        elif field == 'J':
-            return 'A'
         else:
-            return ''
+            comp = None
+        # Read metadata.
+        return self._opmd_reader.read_field_metadata(iteration, field, comp)

visualpic/data_reading/folder_scanners.py

@@ -12,7 +12,7 @@
 
 import numpy as np
 from h5py import File as H5F
-from openpmd_viewer.openpmd_timeseries.data_reader import DataReader
+from .openpmd_data_reader import OpenPMDDataReader
 
 import visualpic.data_reading.field_readers as fr
 import visualpic.data_reading.particle_readers as pr
@@ -77,7 +77,7 @@ def __init__(self, opmd_backend='h5py'):
             Possible values are 'h5py' or 'openpmd-api'.
 
         """
-        self.opmd_reader = DataReader(opmd_backend)
+        self.opmd_reader = OpenPMDDataReader(opmd_backend)
         self.field_reader = fr.OpenPMDFieldReader(self.opmd_reader)
         self.particle_reader = pr.OpenPMDParticleReader(self.opmd_reader)
         self.unit_converter = uc.OpenPMDUnitConverter()

visualpic/data_reading/openpmd_data_reader.py

@@ -0,0 +1,254 @@
+"""
+This module contains the `OpenPMDDataReader`, a class derived from the
+openPMD-viewer `DataReader` which has been expanded with a method for returning
+only the field metadata.
+
+Some of the methods below include code from openPMD-viewer (see at
+https://github.com/openPMD/openPMD-viewer).
+"""
+
+import h5py
+from openpmd_viewer.openpmd_timeseries.data_reader import DataReader
+from openpmd_viewer.openpmd_timeseries.data_reader.h5py_reader import (
+    field_reader as fr)
+from openpmd_viewer.openpmd_timeseries import FieldMetaInformation
+
+
+class OpenPMDDataReader(DataReader):
+    """
+    Class perfroming the access to openPMD data files.
+
+    This class is derived from the original openPMD `DataReader` to extend its
+    functionality and provide a method which returns only the field metadata.
+    """
+
+    def __init__(self, backend):
+        """ Initialize class. """
+        super().__init__(backend)
+
+    def read_field_metadata(self, iteration, field_name, component_name):
+        """
+        Read the field metadata.
+
+        Parameters:
+        -----------
+        iteration : int
+            The iteration at which the parameters should be extracted.
+        field_name : str
+            Name of the field (e.g., `'E'`, `'B'`, `'rho'`, etc.).
+        component_name : str
+            Name of the field component (e.g., `'r'`, `'x'`, `'t'`, etc.)
+
+        Returns:
+        --------
+        A dictionary with the time, field and axis metadata.
+        """
+        # Initialize metadata dictionary.
+        md = {}
+
+        # Read basic openPMD parameters.
+        t, params = self.read_openPMD_params(iteration)
+
+        # Time metadata.
+        md['time'] = {}
+        md['time']['value'] = t
+        md['time']['units'] = 's'
+
+        # Field metadata.
+        md['field'] = {}
+        field_geometry = params['fields_metadata'][field_name]['geometry']
+        axis_labels = params['fields_metadata'][field_name]['axis_labels']
+        field_units = determine_field_units(field_name)
+        md['field']['geometry'] = field_geometry
+        md['field']['axis_labels'] = axis_labels
+        md['field']['units'] = field_units
+
+        # Axis metadata.
+        md['axis'] = {}
+        info = self.get_field_meta_info(
+            iteration, field_name, component_name, axis_labels, field_geometry)
+        for axis in axis_labels:
+            md['axis'][axis] = {}
+            md['axis'][axis]['units'] = 'm'
+            md['axis'][axis]['array'] = getattr(info, axis)
+            md['axis'][axis]['spacing'] = getattr(info, 'd'+axis)
+            md['axis'][axis]['min'] = getattr(info, axis+'min')
+            md['axis'][axis]['max'] = getattr(info, axis+'max')
+
+        return md
+
+    def get_field_meta_info(self, iteration, field, comp, axis_labels,
+                            geometry):
+        """ Get the `FieldMetaInformation` of the field. """
+        if self.backend == 'h5py':
+            filename = self.iteration_to_file[iteration]
+            if geometry in ['thetaMode']:
+                return read_circ_field_metadata_h5py(
+                    filename, iteration, field, comp)
+            elif geometry in ["1dcartesian", "2dcartesian", "3dcartesian"]:
+                return read_cartesian_field_metadata_h5py(
+                    filename, iteration, field, comp, axis_labels)
+        elif self.backend == 'openpmd-api':
+            if geometry in ['thetaMode']:
+                return read_circ_field_metadata_io(
+                    self.series, iteration, field, comp)
+            elif geometry in ["1dcartesian", "2dcartesian", "3dcartesian"]:
+                return read_cartesian_field_metadata_io(
+                    self.series, iteration, field, comp)
+
+
+def read_cartesian_field_metadata_h5py(filename, iteration, field_name,
+                                       component_name, axis_labels):
+    """
+    Get `FieldMetaInformation` of cartesian field using `h5py` backend.
+
+    Parameters:
+    -----------
+    filename : str
+        The absolute path to the HDF5 file.
+    iteration : int
+        The iteration at which to obtain the data.
+    field_name : string
+       Which field to extract.
+    component_name : string, optional
+       Which component of the field to extract.
+    axis_labels: list of strings
+       The name of the dimensions of the array (e.g. ['x', 'y', 'z']).
+    """
+    # Open the HDF5 file
+    dfile = h5py.File(filename, 'r')
+    # Extract the dataset and and corresponding group
+    if component_name is None:
+        field_path = field_name
+    else:
+        field_path = fr.join_infile_path(field_name, component_name)
+    group, dset = fr.find_dataset(dfile, iteration, field_path)
+
+    # Extract the metainformation
+    shape = list(fr.get_shape(dset))
+    axes = {i: axis_labels[i] for i in range(len(axis_labels))}
+    info = FieldMetaInformation(
+        axes, shape, group.attrs['gridSpacing'],
+        group.attrs['gridGlobalOffset'], group.attrs['gridUnitSI'],
+        dset.attrs['position'])
+    return info
+
+
+def read_cartesian_field_metadata_io(series, iteration, field_name,
+                                     component_name, axis_labels):
+    """
+    Get `FieldMetaInformation` of cartesian field using `io` backend.
+
+    Parameters:
+    -----------
+    series : openpmd_api.Series
+        An open, readable openPMD-api series object.
+    iteration : int
+        The iteration at which to obtain the data.
+    field_name : string
+       Which field to extract.
+    component_name : string, optional
+       Which component of the field to extract.
+    axis_labels: list of strings
+       The name of the dimensions of the array (e.g. ['x', 'y', 'z']).
+    """
+    it = series.iterations[iteration]
+
+    # Extract the dataset and and corresponding group
+    field = it.meshes[field_name]
+    if field.scalar:
+        component = next(field.items())[1]
+    else:
+        component = field[component_name]
+
+    # Extract the metainformation
+    shape = component.shape
+    axes = {i: axis_labels[i] for i in range(len(axis_labels))}
+    info = FieldMetaInformation(
+        axes, shape,
+        field.grid_spacing, field.grid_global_offset,
+        field.grid_unit_SI, component.position)
+    return info
+
+
+def read_circ_field_metadata_h5py(filename, iteration, field_name,
+                                  component_name):
+    """
+    Get `FieldMetaInformation` of thetaMode field using `h5py` backend.
+
+    Parameters:
+    -----------
+    filename : str
+        The absolute path to the HDF5 file.
+    iteration : int
+        The iteration at which to obtain the data.
+    field_name : string
+       Which field to extract.
+    component_name : string, optional
+       Which component of the field to extract.
+    """
+    # Open the HDF5 file
+    dfile = h5py.File(filename, 'r')
+    # Extract the dataset and and corresponding group
+    if component_name is None:
+        field_path = field_name
+    else:
+        field_path = fr.join_infile_path(field_name, component_name)
+    group, dset = fr.find_dataset(dfile, iteration, field_path)
+
+    # Extract the metainformation
+    Nm, Nr, Nz = fr.get_shape(dset)
+    info = FieldMetaInformation(
+        {0: 'r', 1: 'z'}, (Nr, Nz),
+        group.attrs['gridSpacing'], group.attrs['gridGlobalOffset'],
+        group.attrs['gridUnitSI'], dset.attrs['position'], thetaMode=True)
+
+    return info
+
+
+def read_circ_field_metadata_io(series, iteration, field_name, component_name):
+    """
+    Get `FieldMetaInformation` of thetaMode field using `io` backend.
+
+    Parameters:
+    -----------
+    series : openpmd_api.Series
+        An open, readable openPMD-api series object.
+    iteration : int
+        The iteration at which to obtain the data.
+    field_name : string
+       Which field to extract.
+    component_name : string, optional
+       Which component of the field to extract.
+    """
+    it = series.iterations[iteration]
+
+    # Extract the dataset and and corresponding group
+    field = it.meshes[field_name]
+    if field.scalar:
+        component = next(field.items())[1]
+    else:
+        component = field[component_name]
+
+    # Extract the metainformation
+    Nm, Nr, Nz = component.shape
+    info = FieldMetaInformation(
+        {0: 'r', 1: 'z'}, (Nr, Nz),
+        field.grid_spacing, field.grid_global_offset,
+        field.grid_unit_SI, component.position, thetaMode=True)
+    return info
+
+
+def determine_field_units(field_name):
+    """ Return the corresponding units of the field. """
+    # TODO: Make more robust implementation using unit_dimension attributes.
+    if field_name == 'E':
+        return 'V/m'
+    elif field_name == 'B':
+        return 'T'
+    elif field_name == 'rho':
+        return 'C/m^3'
+    elif field_name == 'J':
+        return 'A'
+    else:
+        return ''