Add new data encapsulation

gstools/__init__.py

@@ -33,6 +33,7 @@
 
 .. autosummary::
    SRF
+   Mesh
 
 Covariance Base-Class
 ^^^^^^^^^^^^^^^^^^^^^
@@ -100,6 +101,7 @@
 """
 
 from gstools import field, variogram, random, covmodel, tools, krige, transform
+from gstools.field.base import Mesh
 from gstools.field import SRF
 from gstools.tools import (
     vtk_export,

gstools/field/__init__.py

@@ -17,10 +17,12 @@
 
 .. autosummary::
    SRF
+   Mesh
 
 ----
 """
 
+from gstools.field.base import Mesh
 from gstools.field.srf import SRF
 
-__all__ = ["SRF"]
+__all__ = ["SRF", "Mesh"]

gstools/field/base.py

@@ -8,6 +8,7 @@
 
 .. autosummary::
    Field
+   Mesh
 """
 # pylint: disable=C0103
 
@@ -26,33 +27,264 @@
 )
 from gstools.tools.geometric import pos2xyz, xyz2pos
 
-__all__ = ["Field"]
+__all__ = ["Field", "Mesh"]
 
 
-class Field:
-    """A field base class for random and kriging fields ect.
+def value_type(mesh_type, shape):
+    """Determine the value type ("scalar" or "vector")"""
+    r = "scalar"
+    if mesh_type == "unstructured":
+        # TODO is this the right place for doing these checks?!
+        if len(shape) == 2 and 2 <= shape[0] <= 3:
+            r = "vector"
+    else:
+        # for very small (2-3 points) meshes, this could break
+        # a 100% solution would require dim. information.
+        if len(shape) == shape[0] + 1:
+            r = "vector"
+    return r
+
+
+class Mesh:
+    """A base class encapsulating field data.
+
+    It holds a position array, which define the spatial locations of the
+    field values.
+    It can hold multiple fields in the :any:`self.point_data` dict. This assumes
+    that each field is defined on the same positions.
+    The mesh type must also be specified.
+
+    Parameters
+    ----------
+    pos : :class:`numpy.ndarray`, optional
+        positions of the field values
+    name : :any:`str`, optional
+        key of the field values
+    values : :any:`list`, optional
+        a list of the values of the fields
+    mesh_type : :any:`str`, optional
+        the type of mesh on which the field is defined on, can be
+        * unstructured
+        * structured
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from gstools import Mesh
+    >>> pos = np.linspace(0., 100., 40)
+    >>> z = np.random.random(40)
+    >>> z2 = np.random.random(40)
+    >>> mesh = Mesh((pos,))
+    >>> mesh.add_field(z, "test_field1")
+    >>> mesh.add_field(z2, "test_field2")
+    >>> mesh.set_default_field("test_field2")
+    >>> print(mesh.field)
+
+    """
+
+    def __init__(
+        self, pos=None, name="field", values=None, mesh_type="unstructured",
+    ):
+        # mesh_type needs a special setter, therefore, `set_field_data` is not
+        # used here
+        self.mesh_type = mesh_type
+
+        # the pos/ points of the mesh
+        self._pos = pos
+
+        # data stored at each pos/ point, the "fields"
+        if values is not None:
+            self.point_data = {name: values}
+        else:
+            self.point_data = {}
+
+        # data valid for the global field
+        self.field_data = {}
+
+        self.set_field_data("default_field", name)
+
+        self.field_data["mesh_type"] = mesh_type
+
+    def set_field_data(self, name, value):
+        """Add an attribute to this instance and add it the `field_data`
+
+        This helper method is used to create attributes for easy access
+        while using this class, but it also adds an entry to the dictionary
+        `field_data`, which is used for exporting the data.
+        """
+        setattr(self, name, value)
+        self.field_data[name] = value
+
+    def add_field(
+        self, values, name="field", is_default_field=False,
+    ):
+        """Add a field (point data) to the mesh
+
+        .. note::
+            If no field has existed before calling this method,
+            the given field will be set to the default field.
+
+        .. warning::
+            If point data with the same `name` already exists, it will be
+            overwritten.
+
+        Parameters
+        ----------
+        values : :class:`numpy.ndarray`
+            the point data, has to be the same shape as the mesh
+        name : :class:`str`, optional
+            the name of the point data
+        is_default_field : :class:`bool`, optional
+            is this the default field of the mesh?
+
+        """
+        values = np.array(values)
+        self._check_point_data(values)
+        self.point_data[name] = values
+        # set the default_field to the first field added
+        if len(self.point_data) == 1 or is_default_field:
+            self.set_field_data("default_field", name)
+
+    def __getitem__(self, key):
+        """:any:`numpy.ndarray`: The values of the field."""
+        return self.point_data[key]
+
+    def __setitem__(self, key, value):
+        self.point_data[key] = value
+
+    @property
+    def pos(self):
+        """:any:`numpy.ndarray`: The pos. on which the field is defined."""
+        return self._pos
+
+    @pos.setter
+    def pos(self, value):
+        """
+        Warning: setting new positions deletes all previously stored fields.
+        """
+        self.point_data = {self.default_field: None}
+        self._pos = value
+
+    @property
+    def field(self):
+        """:class:`numpy.ndarray`: The point data of the default field."""
+        return self.point_data[self.default_field]
+
+    @field.setter
+    def field(self, values):
+        self._check_point_data(values)
+        self.point_data[self.default_field] = values
+
+    @property
+    def value_type(self, field="field"):
+        """:any:`str`: The value type of the default field."""
+        if field in self.point_data:
+            r = value_type(self.mesh_type, self.point_data[field].shape)
+        else:
+            r = None
+        return r
+
+    @property
+    def mesh_type(self):
+        """:any:`str`: The mesh type of the fields."""
+        return self._mesh_type
+
+    @mesh_type.setter
+    def mesh_type(self, value):
+        """
+        Warning: setting a new mesh type deletes all previously stored fields.
+        """
+        self._check_mesh_type(value)
+        self.point_data = {}
+        self._mesh_type = value
+
+    def _check_mesh_type(self, mesh_type):
+        if mesh_type != "unstructured" and mesh_type != "structured":
+            raise ValueError("Unknown 'mesh_type': {}".format(mesh_type))
+
+    def _check_point_data(self, values):
+        """Compare field shape to pos shape.
+
+        Parameters
+        ----------
+        values : :class:`numpy.ndarray`
+            the values of the field to be checked
+        """
+        err = True
+        if self.mesh_type == "unstructured":
+            # scalar
+            if len(values.shape) == 1:
+                if values.shape[0] == len(self.pos[0]):
+                    err = False
+            # vector
+            elif len(values.shape) == 2:
+                if (
+                    values.shape[1] == len(self.pos[0])
+                    and 2 <= values.shape[0] <= 3
+                ):
+                    err = False
+            if err:
+                raise ValueError(
+                    "Wrong field shape: {0} does not match mesh shape ({1},)".format(
+                        values.shape, len(self.pos[0])
+                    )
+                )
+        else:
+            # scalar
+            if len(values.shape) == len(self.pos):
+                if all(
+                    [
+                        values.shape[i] == len(self.pos[i])
+                        for i in range(len(self.pos))
+                    ]
+                ):
+                    err = False
+            # vector
+            elif len(values.shape) == len(self.pos) + 1:
+                if all(
+                    [
+                        values.shape[i + 1] == len(self.pos[i])
+                        for i in range(len(self.pos))
+                    ]
+                ) and values.shape[0] == len(self.pos):
+                    err = False
+            if err:
+                raise ValueError(
+                    "Wrong field shape: {0} does not match mesh shape [0/2/3]{1}".format(
+                        list(values.shape),
+                        [len(self.pos[i]) for i in range(len(self.pos))],
+                    )
+                )
+
+
+class Field(Mesh):
+    """A field base class for random and kriging fields, etc.
 
     Parameters
     ----------
     model : :any:`CovModel`
         Covariance Model related to the field.
-    mean : :class:`float`, optional
-        Mean value of the field.
     """
 
-    def __init__(self, model, mean=0.0):
+    def __init__(
+        self,
+        model,
+        pos=None,
+        name="field",
+        values=None,
+        mesh_type="unstructured",
+        mean=0.0,
+    ):
         # initialize attributes
-        self.pos = None
-        self.mesh_type = None
-        self.field = None
+        super().__init__(
+            pos=pos, name=name, values=values, mesh_type=mesh_type
+        )
         # initialize private attributes
-        self._mean = None
         self._model = None
-        self.mean = mean
         self.model = model
-        self._value_type = None
+        self.mean = mean
 
-    def __call__(*args, **kwargs):
+    def __call__(self, *args, **kwargs):
         """Generate the field."""
         pass
 
@@ -74,7 +306,7 @@ def unstructured(self, *args, **kwargs):
 
     def mesh(
         self, mesh, points="centroids", direction="xyz", name="field", **kwargs
-    ):  # pragma: no cover
+    ):
         """Generate a field on a given meshio or ogs5py mesh.
 
         Parameters
@@ -382,18 +614,13 @@ def model(self, model):
                 "Field: 'model' is not an instance of 'gstools.CovModel'"
             )
 
-    @property
-    def value_type(self):
-        """:class:`str`: Type of the field values (scalar, vector)."""
-        return self._value_type
-
     def __str__(self):
         """Return String representation."""
         return self.__repr__()
 
     def __repr__(self):
         """Return String representation."""
-        return "Field(model={0}, mean={1})".format(self.model, self.mean)
+        return "Field(model={0})".format(self.model)
 
 
 if __name__ == "__main__":  # pragma: no cover

gstools/field/srf.py

@@ -82,7 +82,7 @@ def __init__(
         generator="RandMeth",
         **generator_kwargs
     ):
-        super().__init__(model, mean)
+        super().__init__(model, mean=mean)
         # initialize private attributes
         self._generator = None
         self._upscaling = None

gstools/krige/base.py

@@ -69,7 +69,7 @@ def __init__(
         drift_functions=None,
         trend_function=None,
     ):
-        super().__init__(model, mean)
+        super().__init__(model, mean=mean)
         self.krige_var = None
         # initialize private attributes
         self._unbiased = True

gstools/transform/field.py

@@ -335,12 +335,12 @@ def _zinnharvey(field, conn="high", mean=None, var=None):
         mean = np.mean(field)
     if var is None:
         var = np.var(field)
-    field = np.abs((field - mean) / np.sqrt(var))
+    field = np.abs((field - mean) / var)
     field = 2 * erf(field / np.sqrt(2)) - 1
     field = np.sqrt(2) * erfinv(field)
     if conn == "high":
         field = -field
-    return field * np.sqrt(var) + mean
+    return field * var + mean
 
 
 def _normal_force_moments(field, mean=0, var=1):
@@ -508,5 +508,5 @@ def _uniform_to_uquad(field, a=0, b=1):
     y_raw = 3 * field / al + ga
     out = np.zeros_like(y_raw)
     out[y_raw > 0] = y_raw[y_raw > 0] ** (1 / 3)
-    out[y_raw < 0] = -((-y_raw[y_raw < 0]) ** (1 / 3))
+    out[y_raw < 0] = -(-y_raw[y_raw < 0]) ** (1 / 3)
     return out + be

tests/test_condition.py

@@ -46,6 +46,7 @@ def test_simple(self):
             )
             srf = SRF(model, self.mean, seed=19970221)
             srf.set_condition(self.cond_pos[0], self.cond_val, "simple")
+            # TODO should this be possible?!?
             field_1 = srf.unstructured(self.pos[0])
             field_2 = srf.structured(self.pos[0])
             for i, val in enumerate(self.cond_val):