OpenPMD HDF5Reader Plasma subclass

plasmapy/classes/__init__.py

@@ -1,4 +1,5 @@
 from .plasma_base import BasePlasma, GenericPlasma
-from . import sources
 from .plasma_factory import Plasma
+
+from . import sources
 from .species import Species

plasmapy/classes/sources/__init__.py

@@ -1,2 +1,3 @@
 from .plasma3d import Plasma3D
 from .plasmablob import PlasmaBlob
+from .openpmd_hdf5 import HDF5Reader

plasmapy/classes/sources/openpmd_hdf5.py

@@ -0,0 +1,164 @@
+import h5py
+import numpy as np
+import astropy.units as u
+
+from plasmapy.classes import GenericPlasma
+from plasmapy.utils import DataStandardError
+
+import os
+from distutils.version import StrictVersion
+
+
+_OUTDATED_VERSION = "1.1.0"
+_NEWER_VERSION = "2.0.0"
+
+# This is the order what OpenPMD uses to store unit
+# dimensions for a record.
+_UNITS = (u.meter,
+          u.kilogram,
+          u.second,
+          u.ampere,
+          u.Kelvin,
+          u.mol,
+          u.candela)
+
+
+def _fetch_units(openPMD_dims):
+    """
+    Converts a collection of OpenPMD dimensions to astropy.units.
+    """
+
+    units = u.dimensionless_unscaled
+    for factor, unit in zip(openPMD_dims, _UNITS):
+        units *= (unit ** factor)
+    units, *_ = units.compose()
+    return units
+
+
+def _valid_version(openPMD_version,
+                   outdated=_OUTDATED_VERSION, newer=_NEWER_VERSION):
+    """
+    Checks if the passed version is supported or not.
+    """
+
+    parsed_version = StrictVersion(openPMD_version)
+    outdated_version = StrictVersion(outdated)
+    newer_version = StrictVersion(newer)
+    return outdated_version <= parsed_version < newer_version
+
+
+class HDF5Reader(GenericPlasma):
+    def __init__(self, hdf5, **kwargs):
+        """
+        Core class for accessing various attributes on HDF5 files that
+        are based on OpenPMD standards.
+
+        Attributes
+        ----------
+        electric_field : `astropy.units.Quantity`
+            An (x, y, z) array containing electric field data.
+        charge_density : `astropy.units.Quantity`
+            An array containing charge density data.
+
+        Parameters
+        ----------
+        hdf5 : `str`
+            Path to HDF5 file.
+
+        References
+        ----------
+        .. [1] http://openpmd.org/
+        """
+
+        if not os.path.isfile(hdf5):
+            raise FileNotFoundError(f"Could not find file: '{hdf5}'")
+
+        h5 = h5py.File(hdf5)
+        self.h5 = h5
+
+        self._check_valid_openpmd_version()
+
+        self.subname = tuple(self.h5['data'])[0]
+
+    def _check_valid_openpmd_version(self):
+        try:
+            openPMD_version = self.h5.attrs["openPMD"].decode('utf-8')
+            if _valid_version(openPMD_version):
+                return True
+            else:
+                raise DataStandardError(f"We currently only support HDF5 versions"
+                                        f"starting from v{_OUTDATED_VERSION} and "
+                                        f"lower than v{_NEWER_VERSION}. You can "
+                                        f"however convert your HDF5 to a supported "
+                                        f"version. For more information; see "
+                                        f"https://github.com/openPMD/openPMD-updater")
+        except KeyError:
+            raise DataStandardError("Input HDF5 file does not go on with "
+                                    "standards defined by OpenPMD")
+
+    @property
+    def electric_field(self):
+        path = f"data/{self.subname}/fields/E"
+        if path in self.h5:
+            units = _fetch_units(self.h5[path].attrs["unitDimension"])
+            axes = [self.h5[path][axis]
+                    for axis in self.h5[path]]
+            return np.array(axes) * units
+        else:
+            raise AttributeError("No electric field data available "
+                                 "in HDF5 file")
+
+    @property
+    def charge_density(self):
+        path = f"data/{self.subname}/fields/rho"
+        if path in self.h5:
+            units = _fetch_units(self.h5[path].attrs["unitDimension"])
+            return np.array(self.h5[path]) * units
+        else:
+            raise AttributeError("No charge density data available "
+                                 "in HDF5 file")
+
+    @property
+    def magnetic_field(self):
+        path = f"data/{self.subname}/fields/B"
+        if path in self.h5:
+            units = _fetch_units(self.h5[path].attrs["unitDimension"])
+            axes = [self.h5[path][axis]
+                    for axis in self.h5[path]]
+            return np.array(axes) * units
+        else:
+            raise AttributeError("No magnetic field data available "
+                                 "in HDF5 file")
+
+    @property
+    def electric_current(self):
+        path = f"data/{self.subname}/fields/J"
+        if path in self.h5:
+            units = _fetch_units(self.h5[path].attrs["unitDimension"])
+            axes = [self.h5[path][axis]
+                    for axis in self.h5[path]]
+            return np.array(axes) * units
+        else:
+            raise AttributeError("No electric current data available "
+                                 "in HDF5 file")
+
+    @classmethod
+    def is_datasource_for(cls, **kwargs):
+        if "hdf5" not in kwargs:
+            return False
+
+        hdf5 = kwargs.get("hdf5")
+        openPMD = kwargs.get("openPMD")
+
+        isfile = os.path.isfile(hdf5)
+        if not isfile:
+            raise FileNotFoundError(f"Could not find file: '{hdf5}'")
+
+        if "openPMD" not in kwargs:
+            h5 = h5py.File(hdf5)
+            try:
+                openPMD = h5.attrs["openPMD"]
+            except KeyError:
+                openPMD = False
+
+        return openPMD

plasmapy/classes/sources/tests/test_openpmd_hdf5.py

@@ -0,0 +1,124 @@
+from plasmapy.classes.sources import openpmd_hdf5
+from plasmapy.utils import DataStandardError
+from plasmapy.data.test import rootdir
+
+from astropy import units as u
+from typing import Union, Tuple, List
+import os
+import pytest
+
+
+class TestOpenPMD2D:
+    """Test 2D HDF5 dataset based on OpenPMD."""
+    # Downloaded from
+    # https://github.com/openPMD/openPMD-example-datasets/blob/draft/example-2d.tar.gz
+    # per the Creative Commons Zero v1.0 Universal license
+    h5 = openpmd_hdf5.HDF5Reader(hdf5=os.path.join(rootdir, "data00000255.h5"))
+
+    def test_has_electric_field_with_units(self):
+        assert self.h5.electric_field.to(u.V / u.m)
+
+    def test_correct_shape_electric_field(self):
+        assert self.h5.electric_field.shape == (3, 51, 201)
+
+    def test_has_charge_density_with_units(self):
+        assert self.h5.charge_density.to(u.C / u.m**3)
+
+    def test_correct_shape_charge_density(self):
+        assert self.h5.charge_density.shape == (51, 201)
+
+    def test_has_magnetic_field(self):
+        with pytest.raises(AttributeError):
+            self.h5.magnetic_field
+
+    def test_has_electric_current(self):
+        with pytest.raises(AttributeError):
+            self.h5.electric_current
+
+
+class TestOpenPMD3D:
+    """Test 3D HDF5 dataset based on OpenPMD."""
+    # Downloaded from
+    # https://github.com/openPMD/openPMD-example-datasets/blob/draft/example-3d.tar.gz
+    # per the Creative Commons Zero v1.0 Universal license
+    h5 = openpmd_hdf5.HDF5Reader(hdf5=os.path.join(rootdir, "data00000100.h5"))
+
+    def test_has_electric_field_with_units(self):
+        assert self.h5.electric_field.to(u.V / u.m)
+
+    def test_correct_shape_electric_field(self):
+        assert self.h5.electric_field.shape == (3, 26, 26, 201)
+
+    def test_has_charge_density_with_units(self):
+        assert self.h5.charge_density.to(u.C / u.m**3)
+
+    def test_correct_shape_charge_density(self):
+        assert self.h5.charge_density.shape == (26, 26, 201)
+
+    def test_has_magnetic_field(self):
+        with pytest.raises(AttributeError):
+            self.h5.magnetic_field
+
+    def test_has_electric_current(self):
+        with pytest.raises(AttributeError):
+            self.h5.electric_current
+
+
+class TestOpenPMDThetaMode:
+    """Test thetaMode HDF5 dataset based on OpenPMD."""
+    # Downloaded from
+    # https://github.com/openPMD/openPMD-example-datasets/blob/draft/example-thetaMode.tar.gz
+    # per the Creative Commons Zero v1.0 Universal license
+    h5 = openpmd_hdf5.HDF5Reader(hdf5=os.path.join(rootdir, "data00000200.h5"))
+
+    def test_has_electric_field_with_units(self):
+        assert self.h5.electric_field.to(u.V / u.m)
+
+    def test_correct_shape_electric_field(self):
+        assert self.h5.electric_field.shape == (3, 3, 51, 201)
+
+    def test_has_charge_density_with_units(self):
+        assert self.h5.charge_density.to(u.C / u.m**3)
+
+    def test_correct_shape_charge_density(self):
+        assert self.h5.charge_density.shape == (3, 51, 201)
+
+    def test_has_magnetic_field_with_units(self):
+        assert self.h5.magnetic_field.to(u.T)
+
+    def test_correct_shape_magnetic_field(self):
+        assert self.h5.magnetic_field.shape == (3, 3, 51, 201)
+
+    def test_has_electric_current_with_units(self):
+        assert self.h5.electric_current.to(u.A * u.kg / u.m**3)
+
+    def test_correct_shape_electric_current(self):
+        assert self.h5.electric_current.shape == (3, 3, 51, 201)
+
+
+units_test_table = [
+    ((1., 1., 0., -1., 0., 0., 2.),
+     u.m * u.kg / u.amp * u.cd ** 2),
+    ((1, 0, 1, 2, 0, 0, 0),
+     u.m * u.s * u.amp ** 2),
+    ([-3.,  0.,  1.,  1.,  0.,  0.,  0.],
+     u.coulomb / u.m**3),
+    ([2, 1, -3, -2, 0, 0, 0],
+     u.ohm)
+]
+
+
+@pytest.mark.parametrize("openPMD_dims, expected", units_test_table)
+def test_fetch_units(openPMD_dims, expected: Union[Tuple, List]):
+    units = openpmd_hdf5._fetch_units(openPMD_dims)
+    assert units == expected
+
+
+def test_unavailable_hdf5():
+    with pytest.raises(FileNotFoundError):
+        openpmd_hdf5.HDF5Reader(hdf5="this_file_does_not_exist.h5")
+
+
+def test_non_openpmd_hdf5():
+    with pytest.raises(DataStandardError):
+        openpmd_hdf5.HDF5Reader(hdf5=os.path.join(rootdir, "blank.h5"))

plasmapy/classes/tests/test_plasma_factory.py

@@ -1,7 +1,9 @@
 import astropy.units as u
 import numpy as np
+import os
 
 import plasmapy.classes
+from plasmapy.data.test import rootdir
 
 
 class TestPlasma(object):
@@ -30,3 +32,8 @@ def test_PlasmaBlob(self):
                                        particle=particle)
         assert isinstance(blob, plasmapy.classes.sources.PlasmaBlob)
         assert isinstance(blob, plasmapy.classes.BasePlasma)
+
+    def test_HDF5Reader(self):
+        h5 = plasmapy.classes.Plasma(hdf5=os.path.join(rootdir, "data00000255.h5"))
+        assert isinstance(h5, plasmapy.classes.sources.HDF5Reader)
+        assert isinstance(h5, plasmapy.classes.BasePlasma)

plasmapy/data/setup_package.py

@@ -0,0 +1,2 @@
+def get_package_data():
+    return {'plasmapy.data.test': ['*', '*/*']}

plasmapy/data/test/__init__.py

@@ -0,0 +1,7 @@
+import os
+import glob
+
+import plasmapy
+
+rootdir = os.path.join(os.path.dirname(plasmapy.__file__), "data", "test")
+file_list = glob.glob(os.path.join(rootdir, '*.[!p]*'))

plasmapy/utils/__init__.py

@@ -13,6 +13,7 @@
                          InvalidIsotopeError,
                          InvalidElementError,
                          InvalidParticleError,
+                         DataStandardError,
                          PlasmaPyWarning,
                          PhysicsWarning,
                          CouplingWarning,

plasmapy/utils/exceptions.py

@@ -78,6 +78,11 @@ class InvalidParticleError(AtomicError):
     pass
 
 
+class DataStandardError(PlasmaPyError):
+    """An exception for when HDF5 is not defined by OpenPMD standard."""
+    pass
+
+
 # ----------
 # Warnings:
 # ----------

requirements/automated-documentation-tests.txt

@@ -2,6 +2,7 @@ astropy
 coveralls
 mpmath
 roman
+h5py
 lmfit
 matplotlib
 numpy

requirements/environment.txt

@@ -7,6 +7,7 @@ cython
 coveralls
 mpmath
 roman
+h5py
 lmfit
 matplotlib
 sphinx

requirements/environment.yml

@@ -15,6 +15,7 @@ dependencies:
     - coveralls
     - mpmath
     - roman
+    - h5py
     - lmfit
     - matplotlib
     - sphinx

requirements/requirements.txt

@@ -3,6 +3,7 @@ scipy (>= 0.19)
 astropy (>= 2.0)
 mpmath (>= 1.0.0)
 roman (>= 1.4)
+h5py (>= 2.8)
 lmfit (>= 0.9.7)
 matplotlib (>= 2.0)
 cython (>= 0.27.2)

setup.cfg

@@ -9,7 +9,7 @@ edit_on_github = True
 github_project = PlasmaPy/PlasmaPy
 # install_requires should be formatted as a comma-separated list, e.g.:
 # install_requires = astropy, scipy, matplotlib
-install_requires = numpy, scipy, astropy, cython, lmfit, matplotlib, colorama, roman, mpmath
+install_requires = numpy, scipy, astropy, cython, lmfit, matplotlib, colorama, roman, mpmath, h5py
 # version should be PEP386 compatible (http://www.python.org/dev/peps/pep-0386)
 version = 0.1.1.dev