590188220

Build Type

github

Committed by GitHub

Commit Message

<a href="https://github.com/GeoStat-Framework/GSTools/commit/<a class=hub.com/GeoStat-Framework/GSTools/commit/9e01292e0e47201e2669340b4a8bdba5eab51d25">9e01292e0<a href="https://github.com/GeoStat-Framework/GSTools/commit/9e01292e0e47201e2669340b4a8bdba5eab51d25">">Merge </a><a class="double-link" href="https://github.com/GeoStat-Framework/GSTools/commit/<a class="double-link" href="https://github.com/GeoStat-Framework/GSTools/commit/3a3fd47e4756f59dd654d7b5a2586e57a84d40fc">3a3fd47e4</a>">3a3fd47e4</a><a href="https://github.com/GeoStat-Framework/GSTools/commit/9e01292e0e47201e2669340b4a8bdba5eab51d25"> into f8e6339a0">f8e6339a0</a>

Pull Request Pull Request #124: Mean, Normalizer, Trend and Fitting

Run Details

522 of 522 new or added lines in 14 files covered. (100.0%)

2461 of 2702 relevant lines covered (91.08%)

13.66 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

58.0

/gstools/field/tools.py

# -*- coding: utf-8 -*-
"""
GStools subpackage providing tools for Fields.

.. currentmodule:: gstools.field.tools

The following classes and functions are provided

.. autosummary::
   to_vtk_helper
   mesh_call
"""
import numpy as np
import meshio

from gstools.normalizer import Normalizer
from gstools.tools.export import to_vtk, vtk_export
from gstools.tools.misc import list_format


__all__ = ["fmt_mean_norm_trend", "to_vtk_helper", "mesh_call"]


def _fmt_func_val(f_cls, func_val):  # pragma: no cover
    if func_val is None:
        return str(None)
    if callable(func_val):
        return "<function>"  # or format(func_val.__name__)
    if np.size(func_val) > 1:
        return list_format(func_val, prec=f_cls.model._prec)
    return "{0:.{p}}".format(float(func_val), p=f_cls.model._prec)


def _fmt_normalizer(f_cls):  # pragma: no cover
    norm = f_cls.normalizer
    return str(None) if norm.__class__ is Normalizer else norm.name


def fmt_mean_norm_trend(f_cls):  # pragma: no cover
    """Format string repr. for mean, normalizer and trend of a field."""
    args = [
        "mean=" + _fmt_func_val(f_cls, f_cls.mean),
        "normalizer=" + _fmt_normalizer(f_cls),
        "trend=" + _fmt_func_val(f_cls, f_cls.trend),
    ]
    return "".join([", " + arg for arg in args if not arg.endswith("None")])


def to_vtk_helper(
    f_cls, filename=None, field_select="field", fieldname="field"
):  # pragma: no cover
    """Create a VTK/PyVista grid of the field or save it as a VTK file.

    This is an internal helper that will handle saving or creating objects

    Parameters
    ----------
    f_cls : :any:`Field`
        Field class in use.
    filename : :class:`str`
        Filename of the file to be saved, including the path. Note that an
        ending (.vtr or .vtu) will be added to the name. If ``None`` is
        passed, a PyVista dataset of the appropriate type will be returned.
    field_select : :class:`str`, optional
        Field that should be stored. Can be:
        "field", "raw_field", "krige_field", "err_field" or "krige_var".
        Default: "field"
    fieldname : :class:`str`, optional
        Name of the field in the VTK file. Default: "field"
    """
    if f_cls.value_type == "vector":
        if hasattr(f_cls, field_select):
            field = getattr(f_cls, field_select)
        else:
            field = None
        if not (f_cls.pos is None or field is None or f_cls.mesh_type is None):
            suf = ["_X", "_Y", "_Z"]
            fields = {}
            for i in range(f_cls.model.dim):
                fields[fieldname + suf[i]] = field[i]
            if filename is None:
                return to_vtk(f_cls.pos, fields, f_cls.mesh_type)
            else:
                return vtk_export(filename, f_cls.pos, fields, f_cls.mesh_type)
    elif f_cls.value_type == "scalar":
        if hasattr(f_cls, field_select):
            field = getattr(f_cls, field_select)
        else:
            field = None
        if not (f_cls.pos is None or field is None or f_cls.mesh_type is None):
            if filename is None:
                return to_vtk(f_cls.pos, {fieldname: field}, f_cls.mesh_type)
            else:
                return vtk_export(
                    filename, f_cls.pos, {fieldname: field}, f_cls.mesh_type
                )
        else:
            print("Field.to_vtk: '{}' not available.".format(field_select))
    else:
        raise ValueError(
            "Unknown field value type: {}".format(f_cls.value_type)
        )


def mesh_call(
    f_cls, mesh, points="centroids", direction="all", name="field", **kwargs
):
    """Generate a field on a given meshio, ogs5py or pyvista mesh.

    Parameters
    ----------
    mesh : meshio.Mesh or ogs5py.MSH or PyVista mesh
        The given meshio, ogs5py, or PyVista mesh
    points : :class:`str`, optional
        The points to evaluate the field at.
        Either the "centroids" of the mesh cells
        (calculated as mean of the cell vertices) or the "points"
        of the given mesh.
        Default: "centroids"
    direction : :class:`str` or :class:`list`, optional
        Here you can state which direction should be choosen for
        lower dimension. For example, if you got a 2D mesh in xz direction,
        you have to pass "xz". By default, all directions are used.
        One can also pass a list of indices.
        Default: "all"
    name : :class:`str` or :class:`list` of :class:`str`, optional
        Name(s) to store the field(s) in the given mesh as point_data or
        cell_data. If to few names are given, digits will be appended.
        Default: "field"
    **kwargs
        Keyword arguments forwareded to `Field.__call__`.

    Notes
    -----
    This will store the field in the given mesh under the given name,
    if a meshio or PyVista mesh was given.

    See: https://github.com/nschloe/meshio

    See: https://github.com/GeoStat-Framework/ogs5py

    See: https://github.com/pyvista/pyvista
    """
    has_pyvista = False
    has_ogs5py = False

    try:
        import pyvista as pv

        has_pyvista = True
    except ImportError:
        pass
    try:
        import ogs5py as ogs

        has_ogs5py = True
    except ImportError:
        pass

    if isinstance(direction, str) and direction == "all":
        select = list(range(f_cls.dim))
    elif isinstance(direction, str):
        select = _get_select(direction)[: f_cls.dim]
    else:
        select = direction[: f_cls.dim]
    if len(select) < f_cls.dim:
        raise ValueError(
            "Field.mesh: need at least {} direction(s), got '{}'".format(
                f_cls.dim, direction
            )
        )
    # convert pyvista mesh
    if has_pyvista and pv.is_pyvista_dataset(mesh):
        if points == "centroids":
            pnts = mesh.cell_centers().points.T[select]
        else:
            pnts = mesh.points.T[select]
        out = f_cls.unstructured(pos=pnts, **kwargs)
        # Deal with the output
        fields = [out] if isinstance(out, np.ndarray) else out
        if f_cls.value_type == "vector":
            fields = [f.T for f in fields]
        for f_name, field in zip(_names(name, len(fields)), fields):
            mesh[f_name] = field
    # convert ogs5py mesh
    elif has_ogs5py and isinstance(mesh, ogs.MSH):
        if points == "centroids":
            pnts = mesh.centroids_flat.T[select]
        else:
            pnts = mesh.NODES.T[select]
        out = f_cls.unstructured(pos=pnts, **kwargs)
    # convert meshio mesh
    elif isinstance(mesh, meshio.Mesh):
        if points == "centroids":
            # define unique order of cells
            offset = []
            length = []
            mesh_dim = mesh.points.shape[1]
            if mesh_dim < f_cls.dim:
                raise ValueError("Field.mesh: mesh dimension too low!")
            pnts = np.empty((0, mesh_dim), dtype=np.double)
            for cell in mesh.cells:
                pnt = np.mean(mesh.points[cell[1]], axis=1)
                offset.append(pnts.shape[0])
                length.append(pnt.shape[0])
                pnts = np.vstack((pnts, pnt))
            # generate pos for __call__
            pnts = pnts.T[select]
            out = f_cls.unstructured(pos=pnts, **kwargs)
            fields = [out] if isinstance(out, np.ndarray) else out
            if f_cls.value_type == "vector":
                fields = [f.T for f in fields]
            f_lists = []
            for field in fields:
                f_list = []
                for of, le in zip(offset, length):
                    f_list.append(field[of : of + le])
                f_lists.append(f_list)
            for f_name, f_list in zip(_names(name, len(f_lists)), f_lists):
                mesh.cell_data[f_name] = f_list
        else:
            out = f_cls.unstructured(pos=mesh.points.T[select], **kwargs)
            fields = [out] if isinstance(out, np.ndarray) else out
            if f_cls.value_type == "vector":
                fields = [f.T for f in fields]
            for f_name, field in zip(_names(name, len(fields)), fields):
                mesh.point_data[f_name] = field
    else:
        raise ValueError("Field.mesh: Unknown mesh format!")
    return out


def _names(name, cnt):
    name = [name] if isinstance(name, str) else list(name)[:cnt]
    if len(name) < cnt:
        name += [name[-1] + str(i + 1) for i in range(cnt - len(name))]
    return name


def _get_select(direction):
    select = []
    if not (0 < len(direction) < 4):
        raise ValueError(
            "Field.mesh: need 1 to 3 direction(s), got '{}'".format(direction)
        )
    for axis in direction:
        if axis == "x":
            if 0 in select:
                raise ValueError(
                    "Field.mesh: got duplicate directions {}".format(direction)
                )
            select.append(0)
        elif axis == "y":
            if 1 in select:
                raise ValueError(
                    "Field.mesh: got duplicate directions {}".format(direction)
                )
            select.append(1)
        elif axis == "z":
            if 2 in select:
                raise ValueError(
                    "Field.mesh: got duplicate directions {}".format(direction)
                )
            select.append(2)
        else:
            raise ValueError(
                "Field.mesh: got unknown direction {}".format(axis)
            )
    return select

1	# -- coding: utf-8 --
2	"""	15✔
3	GStools subpackage providing tools for Fields.
4
5	.. currentmodule:: gstools.field.tools
6
7	The following classes and functions are provided
8
9	.. autosummary::
10	to_vtk_helper
11	mesh_call
12	"""
13	import numpy as np	15✔
14	import meshio	15✔
15
16	from gstools.normalizer import Normalizer	15✔
17	from gstools.tools.export import to_vtk, vtk_export	15✔
18	from gstools.tools.misc import list_format	15✔
19
20
21	__all__ = ["fmt_mean_norm_trend", "to_vtk_helper", "mesh_call"]	15✔
22
23
24	def _fmt_func_val(f_cls, func_val): # pragma: no cover
25	if func_val is None:
26	return str(None)
27	if callable(func_val):
28	return "<function>" # or format(func_val.__name__)
29	if np.size(func_val) > 1:
30	return list_format(func_val, prec=f_cls.model._prec)
31	return "{0:.{p}}".format(float(func_val), p=f_cls.model._prec)
32
33
34	def _fmt_normalizer(f_cls): # pragma: no cover
35	norm = f_cls.normalizer
36	return str(None) if norm.__class__ is Normalizer else norm.name
37
38
39	def fmt_mean_norm_trend(f_cls): # pragma: no cover
40	"""Format string repr. for mean, normalizer and trend of a field."""
41	args = [
42	"mean=" + _fmt_func_val(f_cls, f_cls.mean),
43	"normalizer=" + _fmt_normalizer(f_cls),
44	"trend=" + _fmt_func_val(f_cls, f_cls.trend),
45	]
46	return "".join([", " + arg for arg in args if not arg.endswith("None")])
47
48
49	def to_vtk_helper(
50	f_cls, filename=None, field_select="field", fieldname="field"
51	): # pragma: no cover
52	"""Create a VTK/PyVista grid of the field or save it as a VTK file.
53
54	This is an internal helper that will handle saving or creating objects
55
56	Parameters
57	----------
58	f_cls : :any:`Field`
59	Field class in use.
60	filename : :class:`str`
61	Filename of the file to be saved, including the path. Note that an
62	ending (.vtr or .vtu) will be added to the name. If ``None`` is
63	passed, a PyVista dataset of the appropriate type will be returned.
64	field_select : :class:`str`, optional
65	Field that should be stored. Can be:
66	"field", "raw_field", "krige_field", "err_field" or "krige_var".
67	Default: "field"
68	fieldname : :class:`str`, optional
69	Name of the field in the VTK file. Default: "field"
70	"""
71	if f_cls.value_type == "vector":
72	if hasattr(f_cls, field_select):
73	field = getattr(f_cls, field_select)
74	else:
75	field = None
76	if not (f_cls.pos is None or field is None or f_cls.mesh_type is None):
77	suf = ["_X", "_Y", "_Z"]
78	fields = {}
79	for i in range(f_cls.model.dim):
80	fields[fieldname + suf[i]] = field[i]
81	if filename is None:
82	return to_vtk(f_cls.pos, fields, f_cls.mesh_type)
83	else:
84	return vtk_export(filename, f_cls.pos, fields, f_cls.mesh_type)
85	elif f_cls.value_type == "scalar":
86	if hasattr(f_cls, field_select):
87	field = getattr(f_cls, field_select)
88	else:
89	field = None
90	if not (f_cls.pos is None or field is None or f_cls.mesh_type is None):
91	if filename is None:
92	return to_vtk(f_cls.pos, {fieldname: field}, f_cls.mesh_type)
93	else:
94	return vtk_export(
95	filename, f_cls.pos, {fieldname: field}, f_cls.mesh_type
96	)
97	else:
98	print("Field.to_vtk: '{}' not available.".format(field_select))
99	else:
100	raise ValueError(
101	"Unknown field value type: {}".format(f_cls.value_type)
102	)
103
104
105	def mesh_call(	15✔
106	f_cls, mesh, points="centroids", direction="all", name="field", **kwargs
107	):
108	"""Generate a field on a given meshio, ogs5py or pyvista mesh.
109
110	Parameters
111	----------
112	mesh : meshio.Mesh or ogs5py.MSH or PyVista mesh
113	The given meshio, ogs5py, or PyVista mesh
114	points : :class:`str`, optional
115	The points to evaluate the field at.
116	Either the "centroids" of the mesh cells
117	(calculated as mean of the cell vertices) or the "points"
118	of the given mesh.
119	Default: "centroids"
120	direction : :class:`str` or :class:`list`, optional
121	Here you can state which direction should be choosen for
122	lower dimension. For example, if you got a 2D mesh in xz direction,
123	you have to pass "xz". By default, all directions are used.
124	One can also pass a list of indices.
125	Default: "all"
126	name : :class:`str` or :class:`list` of :class:`str`, optional
127	Name(s) to store the field(s) in the given mesh as point_data or
128	cell_data. If to few names are given, digits will be appended.
129	Default: "field"
130	**kwargs
131	Keyword arguments forwareded to `Field.__call__`.
132
133	Notes
134	-----
135	This will store the field in the given mesh under the given name,
136	if a meshio or PyVista mesh was given.
137
138	See: https://github.com/nschloe/meshio
139
140	See: https://github.com/GeoStat-Framework/ogs5py
141
142	See: https://github.com/pyvista/pyvista
143	"""
144	has_pyvista = False	15✔
145	has_ogs5py = False	15✔
146
147	try:	15✔
148	import pyvista as pv	15✔
149
150	has_pyvista = True	×
151	except ImportError:	15✔
152	pass	15✔
153	try:	15✔
154	import ogs5py as ogs	15✔
155
156	has_ogs5py = True	×
157	except ImportError:	15✔
158	pass	15✔
159
160	if isinstance(direction, str) and direction == "all":	15✔
161	select = list(range(f_cls.dim))	15✔
162	elif isinstance(direction, str):	×
163	select = _get_select(direction)[: f_cls.dim]	×
164	else:
165	select = direction[: f_cls.dim]	×
166	if len(select) < f_cls.dim:	15✔
167	raise ValueError(	×
168	"Field.mesh: need at least {} direction(s), got '{}'".format(
169	f_cls.dim, direction
170	)
171	)
172	# convert pyvista mesh
173	if has_pyvista and pv.is_pyvista_dataset(mesh):	15✔
174	if points == "centroids":	×
175	pnts = mesh.cell_centers().points.T[select]	×
176	else:
177	pnts = mesh.points.T[select]	×
178	out = f_cls.unstructured(pos=pnts, **kwargs)	×
179	# Deal with the output
180	fields = [out] if isinstance(out, np.ndarray) else out	×
181	if f_cls.value_type == "vector":	×
182	fields = [f.T for f in fields]	×
183	for f_name, field in zip(_names(name, len(fields)), fields):	×
184	mesh[f_name] = field	×
185	# convert ogs5py mesh
186	elif has_ogs5py and isinstance(mesh, ogs.MSH):	15✔
187	if points == "centroids":	×
188	pnts = mesh.centroids_flat.T[select]	×
189	else:
190	pnts = mesh.NODES.T[select]	×
191	out = f_cls.unstructured(pos=pnts, **kwargs)	×
192	# convert meshio mesh
193	elif isinstance(mesh, meshio.Mesh):	15✔
194	if points == "centroids":	15✔
195	# define unique order of cells
196	offset = []	15✔
197	length = []	15✔
198	mesh_dim = mesh.points.shape[1]	15✔
199	if mesh_dim < f_cls.dim:	15✔
200	raise ValueError("Field.mesh: mesh dimension too low!")	×
201	pnts = np.empty((0, mesh_dim), dtype=np.double)	15✔
202	for cell in mesh.cells:	15✔
203	pnt = np.mean(mesh.points[cell[1]], axis=1)	15✔
204	offset.append(pnts.shape[0])	15✔
205	length.append(pnt.shape[0])	15✔
206	pnts = np.vstack((pnts, pnt))	15✔
207	# generate pos for __call__
208	pnts = pnts.T[select]	15✔
209	out = f_cls.unstructured(pos=pnts, **kwargs)	15✔
210	fields = [out] if isinstance(out, np.ndarray) else out	15✔
211	if f_cls.value_type == "vector":	15✔
212	fields = [f.T for f in fields]	×
213	f_lists = []	15✔
214	for field in fields:	15✔
215	f_list = []	15✔
216	for of, le in zip(offset, length):	15✔
217	f_list.append(field[of : of + le])	15✔
218	f_lists.append(f_list)	15✔
219	for f_name, f_list in zip(_names(name, len(f_lists)), f_lists):	15✔
220	mesh.cell_data[f_name] = f_list	15✔
221	else:
222	out = f_cls.unstructured(pos=mesh.points.T[select], **kwargs)	15✔
223	fields = [out] if isinstance(out, np.ndarray) else out	15✔
224	if f_cls.value_type == "vector":	15✔
225	fields = [f.T for f in fields]	×
226	for f_name, field in zip(_names(name, len(fields)), fields):	15✔
227	mesh.point_data[f_name] = field	15✔
228	else:
229	raise ValueError("Field.mesh: Unknown mesh format!")	×
230	return out	15✔
231
232
233	def _names(name, cnt):	15✔
234	name = [name] if isinstance(name, str) else list(name)[:cnt]	15✔
235	if len(name) < cnt:	15✔
236	name += [name[-1] + str(i + 1) for i in range(cnt - len(name))]	×
237	return name	15✔
238
239
240	def _get_select(direction):	15✔
241	select = []	×
242	if not (0 < len(direction) < 4):	×
243	raise ValueError(	×
244	"Field.mesh: need 1 to 3 direction(s), got '{}'".format(direction)
245	)
246	for axis in direction:	×
247	if axis == "x":	×
248	if 0 in select:	×
249	raise ValueError(	×
250	"Field.mesh: got duplicate directions {}".format(direction)
251	)
252	select.append(0)	×
253	elif axis == "y":	×
254	if 1 in select:	×
255	raise ValueError(	×
256	"Field.mesh: got duplicate directions {}".format(direction)
257	)
258	select.append(1)	×
259	elif axis == "z":	×
260	if 2 in select:	×
261	raise ValueError(	×
262	"Field.mesh: got duplicate directions {}".format(direction)
263	)
264	select.append(2)	×
265	else:
266	raise ValueError(	×
267	"Field.mesh: got unknown direction {}".format(axis)
268	)
269	return select	×

GeoStat-Framework / GSTools / 590188220

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous