Merge pull request #141 from GeoStat-Framework/plotter_update

Rework plotting in nD

README.md

@@ -188,7 +188,7 @@ bin_center, gamma = gs.vario_estimate((x, y), field)
 fit_model = gs.Stable(dim=2)
 fit_model.fit_variogram(bin_center, gamma, nugget=False)
 # output
-ax = fit_model.plot(x_max=bin_center[-1])
+ax = fit_model.plot(x_max=max(bin_center))
 ax.scatter(bin_center, gamma)
 print(fit_model)
 ```

docs/source/index.rst

@@ -220,7 +220,7 @@ model again.
     fit_model = gs.Stable(dim=2)
     fit_model.fit_variogram(bin_center, gamma, nugget=False)
     # output
-    ax = fit_model.plot(x_max=bin_center[-1])
+    ax = fit_model.plot(x_max=max(bin_center))
     ax.scatter(bin_center, gamma)
     print(fit_model)
 

examples/01_random_field/07_higher_dimensions.py

@@ -53,9 +53,8 @@
 # In order to "prove" correctness, 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"
+    pos, field, sampling_size=2000, mesh_type="structured"
 )
 model.fit_variogram(bin_center, vario)
 print(model)
@@ -67,9 +66,16 @@
 # 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])
+model.plot(x_max=max(bin_center), 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()
+
+###############################################################################
+# GSTools also provides plotting routines for higher dimensions.
+# Fields are shown by 2D cross-sections, where other dimensions can be
+# controlled via sliders.
+
+srf.plot()

examples/02_cov_model/02_aniso_rotation.py

@@ -52,3 +52,4 @@
 # - in 3D: given by yaw, pitch, and roll (known as
 #   `Tait–Bryan <https://en.wikipedia.org/wiki/Euler_angles#Tait-Bryan_angles>`_
 #   angles)
+# - in nD: See the random field example about higher dimensions

examples/03_variogram/03_directional_2d.py

@@ -16,7 +16,7 @@
 
 angle = np.pi / 8
 model = gs.Exponential(dim=2, len_scale=[10, 5], angles=angle)
-x = y = range(100)
+x = y = range(101)
 srf = gs.SRF(model, seed=123456)
 field = srf((x, y), mesh_type="structured")
 
@@ -56,9 +56,7 @@
 model.plot("vario_axis", axis=1, ax=ax1, x_max=40, label="fit on axis 1")
 ax1.set_title("Fitting an anisotropic model")
 
-srf.plot(ax=ax2)
-ax2.set_aspect("equal")
-
+srf.plot(ax=ax2, show_colorbar=False)
 plt.show()
 
 ###############################################################################

examples/03_variogram/04_directional_3d.py

@@ -17,11 +17,12 @@
 
 dim = 3
 # rotation around z, y, x
-angles = [np.pi / 2, np.pi / 4, np.pi / 8]
+angles = [np.deg2rad(90), np.deg2rad(45), np.deg2rad(22.5)]
 model = gs.Gaussian(dim=3, len_scale=[16, 8, 4], angles=angles)
 x = y = z = range(50)
+pos = (x, y, z)
 srf = gs.SRF(model, seed=1001)
-field = srf.structured((x, y, z))
+field = srf.structured(pos)
 
 ###############################################################################
 # Here we generate the axes of the rotated coordinate system
@@ -37,9 +38,9 @@
 # with the longest correlation length (flattest gradient).
 # Then check the transversal directions and so on.
 
-bins = range(0, 40, 3)
 bin_center, dir_vario, counts = gs.vario_estimate(
-    *([x, y, z], field, bins),
+    pos,
+    field,
     direction=main_axes,
     bandwidth=10,
     sampling_size=2000,
@@ -51,48 +52,45 @@
 ###############################################################################
 # Afterwards we can use the estimated variogram to fit a model to it.
 # Note, that the rotation angles need to be set beforehand.
-#
-# We can use the `counts` of data pairs per bin as weights for the fitting
-# routines to give more attention to areas where more data was available.
-# In order to not introduce to much offset at the origin, we disable
-# fitting the nugget.
 
 print("Original:")
 print(model)
-model.fit_variogram(bin_center, dir_vario, weights=counts, nugget=False)
+model.fit_variogram(bin_center, dir_vario)
 print("Fitted:")
 print(model)
 
 ###############################################################################
-# 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, axis1[0]], [0, axis1[1]], [0, axis1[2]], label="0.")
-ax2.plot([0, axis2[0]], [0, axis2[1]], [0, axis2[2]], label="1.")
-ax2.plot([0, axis3[0]], [0, axis3[1]], [0, axis3[2]], label="2.")
-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.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()
+# Plotting main axes and the fitted directional variogram.
+
+fig = plt.figure(figsize=[10, 5])
+ax1 = fig.add_subplot(121, projection=Axes3D.name)
+ax2 = fig.add_subplot(122)
+
+ax1.plot([0, axis1[0]], [0, axis1[1]], [0, axis1[2]], label="0.")
+ax1.plot([0, axis2[0]], [0, axis2[1]], [0, axis2[2]], label="1.")
+ax1.plot([0, axis3[0]], [0, axis3[1]], [0, axis3[2]], label="2.")
+ax1.set_xlim(-1, 1)
+ax1.set_ylim(-1, 1)
+ax1.set_zlim(-1, 1)
+ax1.set_xlabel("X")
+ax1.set_ylabel("Y")
+ax1.set_zlabel("Z")
+ax1.set_title("Tait-Bryan main axis")
+ax1.legend(loc="lower left")
+
+x_max = max(bin_center)
+ax2.scatter(bin_center, dir_vario[0], label="0. axis")
+ax2.scatter(bin_center, dir_vario[1], label="1. axis")
+ax2.scatter(bin_center, dir_vario[2], label="2. axis")
+model.plot("vario_axis", axis=0, ax=ax2, x_max=x_max, label="fit on axis 0")
+model.plot("vario_axis", axis=1, ax=ax2, x_max=x_max, label="fit on axis 1")
+model.plot("vario_axis", axis=2, ax=ax2, x_max=x_max, label="fit on axis 2")
+ax2.set_title("Fitting an anisotropic model")
+ax2.legend()
 
 plt.show()
+
+###############################################################################
+# Also, let's have a look at the field.
+
+srf.plot()

examples/03_variogram/05_auto_fit_variogram.py

@@ -19,7 +19,7 @@
 
 bin_center, gamma = gs.vario_estimate((x, y), field)
 print("estimated bin number:", len(bin_center))
-print("maximal bin distance:", bin_center[-1])
+print("maximal bin distance:", max(bin_center))
 
 ###############################################################################
 # Fit the variogram with a stable model (no nugget fitted).
@@ -31,5 +31,5 @@
 ###############################################################################
 # Plot the fitting result.
 
-ax = fit_model.plot(x_max=bin_center[-1])
+ax = fit_model.plot(x_max=max(bin_center))
 ax.scatter(bin_center, gamma)

gstools/field/base.py

@@ -268,7 +268,9 @@ def vtk_export(
             fieldname=fieldname,
         )
 
-    def plot(self, field="field", fig=None, ax=None):  # pragma: no cover
+    def plot(
+        self, field="field", fig=None, ax=None, **kwargs
+    ):  # pragma: no cover
         """
         Plot the spatial random field.
 
@@ -283,6 +285,8 @@ def plot(self, field="field", fig=None, ax=None):  # pragma: no cover
         ax : :class:`Axes` or :any:`None`
             Axes to plot on. If `None`, a new one will be added to the figure.
             Default: `None`
+        **kwargs
+            Forwarded to the plotting routine.
         """
         # just import if needed; matplotlib is not required by setup
         from gstools.field.plot import plot_field, plot_vec_field
@@ -294,11 +298,11 @@ def plot(self, field="field", fig=None, ax=None):  # pragma: no cover
             )
 
         elif self.value_type == "scalar":
-            r = plot_field(self, field, fig, ax)
+            r = plot_field(self, field, fig, ax, **kwargs)
 
         elif self.value_type == "vector":
             if self.model.dim == 2:
-                r = plot_vec_field(self, field, fig, ax)
+                r = plot_vec_field(self, field, fig, ax, **kwargs)
             else:
                 raise NotImplementedError(
                     "Streamflow plotting only supported for 2d case."