Empirical Variogram: Update and Refactoring

README.md

@@ -154,7 +154,7 @@ srf = gs.SRF(model, mean=0, seed=19970221)
 field = srf((x, y))
 # estimate the variogram of the field with 40 bins
 bins = np.arange(40)
-bin_center, gamma = gs.vario_estimate_unstructured((x, y), field, bins)
+bin_center, gamma = gs.vario_estimate((x, y), field, bins)
 # fit the variogram with a stable model. (no nugget fitted)
 fit_model = gs.Stable(dim=2)
 fit_model.fit_variogram(bin_center, gamma, nugget=False)

docs/source/index.rst

@@ -182,7 +182,7 @@ model again.
     field = srf((x, y))
     # estimate the variogram of the field with 40 bins
     bins = np.arange(40)
-    bin_center, gamma = gs.vario_estimate_unstructured((x, y), field, bins)
+    bin_center, gamma = gs.vario_estimate((x, y), field, bins)
     # fit the variogram with a stable model. (no nugget fitted)
     fit_model = gs.Stable(dim=2)
     fit_model.fit_variogram(bin_center, gamma, nugget=False)

examples/00_misc/02_check_rand_meth_sampling.py

@@ -58,7 +58,7 @@ def plot_rand_meth_samples(generator):
     ax.set_xlim([0, np.max(x)])
     ax.set_title("Radius samples shown {}/{}".format(sample_in, len(rad)))
     ax.legend()
-    fig.show()
+    plt.show()
 
 
 model = gs.Stable(dim=3, alpha=1.5)

examples/03_variogram/00_fit_variogram.py

@@ -18,7 +18,7 @@
 # Estimate the variogram of the field with 40 bins.
 
 bins = np.arange(40)
-bin_center, gamma = gs.vario_estimate_unstructured((x, y), field, bins)
+bin_center, gamma = gs.vario_estimate((x, y), field, bins)
 
 ###############################################################################
 # Fit the variogram with a stable model (no nugget fitted).

examples/03_variogram/01_variogram_estimation.py

@@ -146,7 +146,7 @@ def generate_transmissivity():
 
 
 bins = np.linspace(0, 10, 50)
-bin_center, gamma = gs.vario_estimate_unstructured(
+bin_center, gamma = gs.vario_estimate(
     (x_u, y_u),
     herten_log_trans.reshape(-1),
     bins,
@@ -209,8 +209,8 @@ def generate_transmissivity():
 # With this much smaller data set, we can immediately estimate the variogram in
 # the x- and y-axis
 
-gamma_x = gs.vario_estimate_structured(herten_trans_skip, direction="x")
-gamma_y = gs.vario_estimate_structured(herten_trans_skip, direction="y")
+gamma_x = gs.vario_estimate_axis(herten_trans_skip, direction="x")
+gamma_y = gs.vario_estimate_axis(herten_trans_skip, direction="y")
 
 ###############################################################################
 # With these two estimated variograms, we can start fitting :any:`Exponential`

examples/03_variogram/02_find_best_model.py

@@ -19,7 +19,7 @@
 # Estimate the variogram of the field with 40 bins and plot the result.
 
 bins = np.arange(40)
-bin_center, gamma = gs.vario_estimate_unstructured((x, y), field, bins)
+bin_center, gamma = gs.vario_estimate((x, y), field, bins)
 
 ###############################################################################
 # Define a set of models to test.

examples/03_variogram/03_multi_vario.py

@@ -0,0 +1,43 @@
+"""
+Multi-field variogram estimation
+--------------------------------
+
+In this example, we demonstrate how to estimate a variogram from multiple
+fields on the same point-set that should have the same statistical properties.
+"""
+import numpy as np
+import gstools as gs
+import matplotlib.pyplot as plt
+
+
+x = np.random.RandomState(19970221).rand(1000) * 100.0
+y = np.random.RandomState(20011012).rand(1000) * 100.0
+model = gs.Exponential(dim=2, var=2, len_scale=8)
+srf = gs.SRF(model, mean=0)
+
+###############################################################################
+# Generate two synthetic fields with an exponential model.
+
+field1 = srf((x, y), seed=19970221)
+field2 = srf((x, y), seed=20011012)
+fields = [field1, field2]
+
+###############################################################################
+# Now we estimate the variograms for both fields individually and then again
+# simultaneously with only one call.
+
+bins = np.arange(40)
+bin_center, gamma1 = gs.vario_estimate((x, y), field1, bins)
+bin_center, gamma2 = gs.vario_estimate((x, y), field2, bins)
+bin_center, gamma = gs.vario_estimate((x, y), fields, bins)
+
+###############################################################################
+# Now we demonstrate that the mean variogram from both fields coincides
+# with the joined estimated one.
+
+plt.plot(bin_center, gamma1, label="field 1")
+plt.plot(bin_center, gamma2, label="field 2")
+plt.plot(bin_center, gamma, label="joined fields")
+plt.plot(bin_center, 0.5 * (gamma1 + gamma2), ":", label="field 1+2 mean")
+plt.legend()
+plt.show()

examples/03_variogram/04_directional_2d.py

@@ -0,0 +1,53 @@
+"""
+Directional variogram estimation in 2D
+--------------------------------------
+
+In this example, we demonstrate how to estimate a directional variogram by
+setting the direction angles in 2D.
+"""
+import numpy as np
+import gstools as gs
+from matplotlib import pyplot as plt
+
+###############################################################################
+# Generating synthetic field with anisotropy and a rotation of 22.5 degree.
+
+angle = np.pi / 8
+model = gs.Exponential(dim=2, len_scale=[10, 5], angles=angle)
+x = y = range(100)
+srf = gs.SRF(model, seed=123456)
+field = srf((x, y), mesh_type="structured")
+
+###############################################################################
+# Now we are going to estimate a directional variogram with an angular
+# tolerance of 11.25 degree and a bandwith of 1.
+# We provide the rotation angle of the covariance model and the orthogonal
+# direction by adding 90 degree.
+
+bins = range(0, 40, 3)
+bin_c, vario, cnt = gs.vario_estimate(
+    *((x, y), field, bins),
+    angles=(angle, angle + np.pi / 2),  # main dir. and transversal dir.
+    angles_tol=np.pi / 16,
+    bandwidth=1.0,
+    mesh_type="structured",
+    return_counts=True,
+)
+
+###############################################################################
+# Plotting.
+
+fig, (ax1, ax2) = plt.subplots(1, 2, figsize=[10, 5])
+
+ax1.plot(bin_c, vario[0], label="emp. vario: pi/8")
+ax1.plot(bin_c, vario[1], label="emp. vario: pi*5/8")
+ax1.legend(loc="lower right")
+
+srf.plot(ax=ax2)
+ax2.set_aspect("equal")
+
+plt.show()
+
+###############################################################################
+# Without fitting a model, we see that the correlation length in the main
+# direction is greater than the transversal one.

examples/03_variogram/05_directional_3d.py

@@ -0,0 +1,77 @@
+"""
+Directional variogram estimation in 3D
+--------------------------------------
+
+In this example, we demonstrate how to estimate a directional variogram by
+setting the estimation directions in 3D.
+"""
+import numpy as np
+import gstools as gs
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+
+###############################################################################
+# Generating synthetic field with anisotropy and rotation by Tait-Bryan angles.
+
+dim = 3
+# rotation around z, y, x
+angles = [np.pi / 2, np.pi / 4, np.pi / 8]
+model = gs.Gaussian(dim=3, len_scale=[16, 8, 4], angles=angles)
+x = y = z = range(50)
+srf = gs.SRF(model, seed=1001)
+field = srf.structured((x, y, z))
+
+###############################################################################
+# Here we generate the rotated coordinate system to get an impression what
+# the rotation angles do.
+
+x1, x2, x3 = (1, 0, 0), (0, 1, 0), (0, 0, 1)
+ret = np.array(gs.field.tools.rotate_mesh(dim, angles, x1, x2, x3))
+dir0 = ret[:, 0]  # main direction
+dir1 = ret[:, 1]  # first lateral direction
+dir2 = ret[:, 2]  # second lateral direction
+
+###############################################################################
+# Now we estimate the variogram along the main axis. When the main axis is
+# unknown, one would need to sample multiple directions and look for the one
+# with the longest correlation length (flattest gradient).
+# Then check the transversal directions and so on.
+
+bins = range(0, 40, 3)
+bin_c, vario = gs.vario_estimate(
+    *([x, y, z], field, bins),
+    direction=(dir0, dir1, dir2),
+    bandwidth=10,
+    sampling_size=2000,
+    sampling_seed=1001,
+    mesh_type="structured"
+)
+
+###############################################################################
+# Plotting.
+
+fig = plt.figure(figsize=[15, 5])
+ax1 = fig.add_subplot(131)
+ax2 = fig.add_subplot(132, projection=Axes3D.name)
+ax3 = fig.add_subplot(133)
+
+srf.plot(ax=ax1)
+ax1.set_aspect("equal")
+
+ax2.plot([0, dir0[0]], [0, dir0[1]], [0, dir0[2]], label="1.")
+ax2.plot([0, dir1[0]], [0, dir1[1]], [0, dir1[2]], label="2.")
+ax2.plot([0, dir2[0]], [0, dir2[1]], [0, dir2[2]], label="3.")
+ax2.set_xlim(-1, 1)
+ax2.set_ylim(-1, 1)
+ax2.set_zlim(-1, 1)
+ax2.set_xlabel("X")
+ax2.set_ylabel("Y")
+ax2.set_zlabel("Z")
+ax2.set_title("Tait-Bryan main axis")
+ax2.legend(loc="lower left")
+
+ax3.plot(bin_c, vario[0], label="1. axis")
+ax3.plot(bin_c, vario[1], label="2. axis")
+ax3.plot(bin_c, vario[2], label="3. axis")
+ax3.legend()
+plt.show()

gstools/__init__.py

@@ -95,8 +95,8 @@
 .. currentmodule:: gstools.variogram
 
 .. autosummary::
-   vario_estimate_structured
-   vario_estimate_unstructured
+   vario_estimate
+   vario_estimate_axis
 """
 
 from gstools import field, variogram, random, covmodel, tools, krige, transform
@@ -110,6 +110,8 @@
     to_vtk_unstructured,
 )
 from gstools.variogram import (
+    vario_estimate,
+    vario_estimate_axis,
     vario_estimate_structured,
     vario_estimate_unstructured,
 )
@@ -154,7 +156,12 @@
     "TPLStable",
 ]
 
-__all__ += ["vario_estimate_structured", "vario_estimate_unstructured"]
+__all__ += [
+    "vario_estimate",
+    "vario_estimate_axis",
+    "vario_estimate_structured",
+    "vario_estimate_unstructured",
+]
 
 __all__ += [
     "SRF",

gstools/covmodel/base.py

@@ -1219,9 +1219,7 @@ def name(self):
     @property
     def do_rotation(self):
         """:any:`bool`: State if a rotation is performed."""
-        return (
-            not np.all(np.isclose(self.angles, 0.0))
-        )
+        return not np.all(np.isclose(self.angles, 0.0))
 
     @property
     def is_isotropic(self):

gstools/field/plot.py

@@ -15,7 +15,6 @@
 from scipy import interpolate as inter
 import matplotlib.pyplot as plt
 from matplotlib.widgets import Slider, RadioButtons
-from mpl_toolkits.mplot3d import Axes3D
 from gstools.tools import pos2xyz
 from gstools.covmodel.plot import _get_fig_ax
 
@@ -91,7 +90,7 @@ def _plot_2d(pos, field, mesh_type, fig=None, ax=None):  # pragma: no cover
 def _plot_3d(pos, field, mesh_type, fig=None, ax=None):  # pragma: no cover
     """Plot 3D field."""
     dir1, dir2 = np.mgrid[0:1:51j, 0:1:51j]
-    levels = np.linspace(field.min(), field.max(), 256, endpoint=True)
+    levels = np.linspace(field.min(), field.max(), 100, endpoint=True)
 
     x_min = pos[0].min()
     x_max = pos[0].max()
@@ -110,7 +109,7 @@ def _plot_3d(pos, field, mesh_type, fig=None, ax=None):  # pragma: no cover
         "y": [y_min, y_max, y_range, y_step],
         "z": [z_min, z_max, z_range, z_step],
     }
-    fig, ax = _get_fig_ax(fig, ax, Axes3D.name)
+    fig, ax = _get_fig_ax(fig, ax)
     title = "Field 3D " + mesh_type + ": " + str(field.shape)
     fig.subplots_adjust(left=0.2, right=0.8, bottom=0.25)
     sax = plt.axes([0.15, 0.1, 0.65, 0.03])
@@ -122,7 +121,7 @@ def _plot_3d(pos, field, mesh_type, fig=None, ax=None):  # pragma: no cover
         valinit=z_min + z_range / 2.0,
         valstep=z_step,
     )
-    rax = plt.axes([0.05, 0.5, 0.1, 0.15])
+    rax = plt.axes([0.05, 0.7, 0.1, 0.15])
     radio = RadioButtons(rax, ("x slice", "y slice", "z slice"), active=2)
     z_dir_tmp = "z"
     # create container
@@ -155,13 +154,11 @@ def get_plane(z_val_in, z_dir):
             plane = inter.griddata(
                 pos, field, (x_io, y_io, z_io), method="linear"
             )
-        if z_dir == "z":
-            z_io = plane
+        if z_dir == "x":
+            return y_io, z_io, plane
         elif z_dir == "y":
-            y_io = plane
-        else:
-            x_io = plane
-        return x_io, y_io, z_io
+            return x_io, z_io, plane
+        return x_io, y_io, plane
 
     def update(__):
         """Widget update."""
@@ -188,17 +185,20 @@ def update(__):
             vmin=field.min(),
             vmax=field.max(),
             levels=levels,
-            zdir=z_dir_in,
-            offset=z_val,
         )
-        cont.cmap.set_under("k", alpha=0.0)
-        cont.cmap.set_bad("k", alpha=0.0)
-        ax.set_xlabel("X")
-        ax.set_ylabel("Y")
-        ax.set_zlabel("Z")
+        # cont.cmap.set_under("k", alpha=0.0)
+        # cont.cmap.set_bad("k", alpha=0.0)
+        if z_dir_in == "x":
+            ax.set_xlabel("Y")
+            ax.set_ylabel("Z")
+        elif z_dir_in == "y":
+            ax.set_xlabel("X")
+            ax.set_ylabel("Z")
+        else:
+            ax.set_xlabel("X")
+            ax.set_ylabel("Y")
         ax.set_xlim([x_min, x_max])
         ax.set_ylim([y_min, y_max])
-        ax.set_zlim([z_min, z_max])
         ax.set_title(title)
         fig.canvas.draw_idle()
         return cont

gstools/krige/base.py

@@ -634,7 +634,9 @@ def __repr__(self):
         """Return String representation."""
         return (
             "{0}(model={1}, cond_no={2}".format(
-                self.name, self.model.name, self.cond_no,
+                self.name,
+                self.model.name,
+                self.cond_no,
             )
             + ")"
         )