Merge pull request #112 from GeoStat-Framework/field_refactor

[Field] refactoring and update

docs/requirements_doc.txt

@@ -3,3 +3,4 @@ sphinx-gallery
 matplotlib
 pyvista
 pykrige
+meshzoo

examples/01_random_field/05_mesh_ensemble.py

@@ -0,0 +1,94 @@
+"""
+Generating fields on meshes
+---------------------------
+
+GSTools provides an interface for meshes, to support
+`meshio <https://github.com/nschloe/meshio>`_ and
+`ogs5py <https://github.com/GeoStat-Framework/ogs5py>`_ meshes.
+
+When using `meshio`, the generated fields will be stored imediatly in the mesh
+container.
+
+One has two options to generate a field on a given mesh:
+
+- `points="points"` will generate a field on the mesh points
+- `points="centroids"` will generate a field on the cell centroids
+
+In this example, we will generate a simple mesh with the aid of
+`meshzoo <https://github.com/nschloe/meshzoo>`_.
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.tri as tri
+import meshzoo
+import meshio
+import gstools as gs
+
+# generate a triangulated hexagon with meshzoo
+points, cells = meshzoo.ngon(6, 4)
+mesh = meshio.Mesh(points, {"triangle": cells})
+
+###############################################################################
+# Now we prepare the SRF class as always. We will generate an ensemble of
+# fields on the generated mesh.
+
+# number of fields
+fields_no = 12
+# model setup
+model = gs.Gaussian(dim=2, len_scale=0.5)
+srf = gs.SRF(model, mean=1)
+
+###############################################################################
+# To generate fields on a mesh, we provide a separate method: :any:`SRF.mesh`.
+# First we generate fields on the mesh-centroids controlled by a seed.
+# You can specify the field name by the keyword `name`.
+
+for i in range(fields_no):
+    srf.mesh(mesh, points="centroids", name="c-field-{}".format(i), seed=i)
+
+###############################################################################
+# Now we generate fields on the mesh-points again controlled by a seed.
+
+for i in range(fields_no):
+    srf.mesh(mesh, points="points", name="p-field-{}".format(i), seed=i)
+
+###############################################################################
+# To get an impression we now want to plot the generated fields.
+# Luckily, matplotlib supports triangular meshes.
+
+
+triangulation = tri.Triangulation(points[:, 0], points[:, 1], cells)
+# figure setup
+cols = 4
+rows = int(np.ceil(fields_no / cols))
+
+###############################################################################
+# Cell data can be easily visualized with matplotlibs `tripcolor`.
+# To highlight the cell structure, we use `triplot`.
+
+fig = plt.figure(figsize=[2 * cols, 2 * rows])
+for i, field in enumerate(mesh.cell_data, 1):
+    ax = fig.add_subplot(rows, cols, i)
+    ax.tripcolor(triangulation, mesh.cell_data[field][0])
+    ax.triplot(triangulation, linewidth=0.5, color="k")
+    ax.set_aspect("equal")
+fig.tight_layout()
+
+###############################################################################
+# Point data is plotted via `tricontourf`.
+
+fig = plt.figure(figsize=[2 * cols, 2 * rows])
+for i, field in enumerate(mesh.point_data, 1):
+    ax = fig.add_subplot(rows, cols, i)
+    ax.tricontourf(triangulation, mesh.point_data[field])
+    ax.triplot(triangulation, linewidth=0.5, color="k")
+    ax.set_aspect("equal")
+fig.tight_layout()
+plt.show()
+
+###############################################################################
+# Last but not least, `meshio` can be used for what is does best: Exporting.
+# Tada!
+
+mesh.write("mesh_ensemble.vtk")

examples/01_random_field/06_higher_dimensions.py

@@ -0,0 +1,75 @@
+"""
+Higher Dimensions
+-----------------
+
+GSTools provides experimental support for higher dimensions.
+
+Anisotropy is the same as in lower dimensions:
+
+- in `n` dimensions we need `(n-1)` anisotropy ratios
+
+Rotation on the other hand is a bit more complex.
+With increasing dimensions more and more rotation angles are added in order
+to properply describe the rotated axes of anisotropy.
+
+By design the first rotation angles coincide with the lower ones:
+
+- 2D (rotation in x-y plane) -> 3D: first angle describes xy-plane rotation
+- 3D (tait-bryan angles) -> 4D: first 3 angles coincide with tait-bryan angles
+
+By increasing the dimension from `n` to `(n+1)`, `n` angles are added:
+
+- 2D (1 angle) -> 3D: 3 angles (2 added)
+- 3D (3 angles) -> 4D: 6 angles (3 added)
+
+the following list of rotation-planes are described by the list of
+angles in the model:
+
+1. x-y plane
+2. x-z plane
+3. y-z plane
+4. x-v plane
+5. y-v plane
+6. z-v plane
+7. ...
+
+The rotation direction in these planes have alternating signs
+in order to match tait-bryan in 3D.
+
+Let's have a look at a 4D example, where we naively add a 4th dimension.
+"""
+
+import matplotlib.pyplot as plt
+import gstools as gs
+
+dim = 4
+size = 20
+pos = [range(size)] * dim
+model = gs.Exponential(dim=dim, len_scale=5)
+srf = gs.SRF(model, seed=20170519)
+field = srf.structured(pos)
+
+###############################################################################
+# In order to "prove" correctnes, we can calculate an empirical variogram
+# of the generated field and fit our model to it.
+
+bin_edges = range(size)
+bin_center, vario = gs.vario_estimate(
+    pos, field, bin_edges, sampling_size=2000, mesh_type="structured"
+)
+model.fit_variogram(bin_center, vario)
+print(model)
+
+###############################################################################
+# As you can see, the estimated variance and length scale match our input
+# quite good.
+#
+# Let's have a look at the fit and a x-y cross-section of the 4D field:
+
+f, a = plt.subplots(1, 2, gridspec_kw={"width_ratios": [2, 1]}, figsize=[9, 3])
+model.plot(x_max=size + 1, ax=a[0])
+a[0].scatter(bin_center, vario)
+a[1].imshow(field[:, :, 0, 0].T, origin="lower")
+a[0].set_title("isotropic empirical variogram with fitted model")
+a[1].set_title("x-y cross-section")
+f.show()

examples/03_variogram/00_fit_variogram.py

@@ -30,5 +30,5 @@
 # Plot the fitting result.
 
 ax = fit_model.plot(x_max=40)
-ax.plot(bin_center, gamma)
+ax.scatter(bin_center, gamma)
 print(fit_model)

examples/03_variogram/02_find_best_model.py

@@ -25,22 +25,24 @@
 # Define a set of models to test.
 
 models = {
-    "gaussian": gs.Gaussian,
-    "exponential": gs.Exponential,
-    "matern": gs.Matern,
-    "stable": gs.Stable,
-    "rational": gs.Rational,
-    "linear": gs.Linear,
-    "circular": gs.Circular,
-    "spherical": gs.Spherical,
+    "Gaussian": gs.Gaussian,
+    "Exponential": gs.Exponential,
+    "Matern": gs.Matern,
+    "Stable": gs.Stable,
+    "Rational": gs.Rational,
+    "Linear": gs.Linear,
+    "Circular": gs.Circular,
+    "Spherical": gs.Spherical,
+    "SuperSpherical": gs.SuperSpherical,
+    "JBessel": gs.JBessel,
 }
 scores = {}
 
 ###############################################################################
 # Iterate over all models, fit their variogram and calculate the r2 score.
 
 # plot the estimated variogram
-plt.scatter(bin_center, gamma, label="data")
+plt.scatter(bin_center, gamma, color="k", label="data")
 ax = plt.gca()
 
 # fit all models to the estimated variogram

examples/03_variogram/04_directional_2d.py

@@ -48,8 +48,8 @@
 
 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=[10, 5])
 
-ax1.plot(bin_center, dir_vario[0], label="emp. vario: pi/8")
-ax1.plot(bin_center, dir_vario[1], label="emp. vario: pi*5/8")
+ax1.scatter(bin_center, dir_vario[0], label="emp. vario: pi/8")
+ax1.scatter(bin_center, dir_vario[1], label="emp. vario: pi*5/8")
 ax1.legend(loc="lower right")
 
 model.plot("vario_axis", axis=0, ax=ax1, x_max=40, label="fit on axis 0")

examples/03_variogram/05_directional_3d.py

@@ -86,12 +86,12 @@
 ax2.set_title("Tait-Bryan main axis")
 ax2.legend(loc="lower left")
 
-ax3.plot(bin_center, dir_vario[0], label="0. axis")
-ax3.plot(bin_center, dir_vario[1], label="1. axis")
-ax3.plot(bin_center, dir_vario[2], label="2. axis")
-model.plot("vario_axis", axis=0, ax=ax3, label="fit on axis 0", linestyle=":")
-model.plot("vario_axis", axis=1, ax=ax3, label="fit on axis 1", linestyle=":")
-model.plot("vario_axis", axis=2, ax=ax3, label="fit on axis 2", linestyle=":")
+ax3.scatter(bin_center, dir_vario[0], label="0. axis")
+ax3.scatter(bin_center, dir_vario[1], label="1. axis")
+ax3.scatter(bin_center, dir_vario[2], label="2. axis")
+model.plot("vario_axis", axis=0, ax=ax3, label="fit on axis 0")
+model.plot("vario_axis", axis=1, ax=ax3, label="fit on axis 1")
+model.plot("vario_axis", axis=2, ax=ax3, label="fit on axis 2")
 ax3.set_title("Fitting an anisotropic model")
 ax3.legend()
 

gstools/__init__.py

@@ -104,7 +104,7 @@
    vario_estimate
    vario_estimate_axis
 """
-
+# Hooray!
 from gstools import field, variogram, random, covmodel, tools, krige, transform
 from gstools.field import SRF
 from gstools.tools import (

gstools/covmodel/plot.py

@@ -24,7 +24,7 @@
 import numpy as np
 
 import gstools
-from gstools.field.tools import reshape_axis_from_struct_to_unstruct
+
 
 __all__ = [
     "plot_variogram",
@@ -70,13 +70,11 @@ def plot_vario_spatial(
     field._value_type = "scalar"
     if x_max is None:
         x_max = 3 * model.len_scale
-    field.mesh_type = "structured"
     x_s = np.linspace(-x_max, x_max) + x_min
-    pos = [x_s] * model.dim
-    x, y, z, shape = reshape_axis_from_struct_to_unstruct(model.dim, *pos)
-    vario = model.vario_spatial([x, y, z][: model.dim]).reshape(shape)
-    field.pos = pos
-    field.field = vario
+    iso_pos, shape = field.pre_pos([x_s] * model.dim, "structured")
+    field.field = model.vario_spatial(model.anisometrize(iso_pos)).reshape(
+        shape
+    )
     return field.plot(fig=fig, ax=ax)
 
 
@@ -88,13 +86,9 @@ def plot_cov_spatial(
     field._value_type = "scalar"
     if x_max is None:
         x_max = 3 * model.len_scale
-    field.mesh_type = "structured"
     x_s = np.linspace(-x_max, x_max) + x_min
-    pos = [x_s] * model.dim
-    x, y, z, shape = reshape_axis_from_struct_to_unstruct(model.dim, *pos)
-    vario = model.cov_spatial([x, y, z][: model.dim]).reshape(shape)
-    field.pos = pos
-    field.field = vario
+    iso_pos, shape = field.pre_pos([x_s] * model.dim, "structured")
+    field.field = model.cov_spatial(model.anisometrize(iso_pos)).reshape(shape)
     return field.plot(fig=fig, ax=ax)
 
 
@@ -106,13 +100,9 @@ def plot_cor_spatial(
     field._value_type = "scalar"
     if x_max is None:
         x_max = 3 * model.len_scale
-    field.mesh_type = "structured"
     x_s = np.linspace(-x_max, x_max) + x_min
-    pos = [x_s] * model.dim
-    x, y, z, shape = reshape_axis_from_struct_to_unstruct(model.dim, *pos)
-    vario = model.cor_spatial([x, y, z][: model.dim]).reshape(shape)
-    field.pos = pos
-    field.field = vario
+    iso_pos, shape = field.pre_pos([x_s] * model.dim, "structured")
+    field.field = model.cor_spatial(model.anisometrize(iso_pos)).reshape(shape)
     return field.plot(fig=fig, ax=ax)