Merge pull request #109 from GeoStat-Framework/cov_refactor

CovModel: Update and Refactoring

docs/source/conf.py

@@ -201,8 +201,8 @@ def setup(app):
     "pointsize": "10pt",
     "papersize": "a4paper",
     "fncychap": "\\usepackage[Glenn]{fncychap}",
+    # 'inputenc': r'\usepackage[utf8]{inputenc}',
 }
-
 # Grouping the document tree into LaTeX files. List of tuples
 # (source start file, target name, title,
 #  author, documentclass [howto, manual, or own class]).

examples/00_misc/README.rst

@@ -3,5 +3,9 @@ Miscellaneous
 
 A few miscellaneous examples
 
-Gallery
--------
+.. only:: html
+
+   Gallery
+   -------
+
+   Below is a gallery of examples

examples/01_random_field/README.rst

@@ -13,5 +13,9 @@ and its discretised modes are evaluated at random frequencies.
 
 GSTools supports arbitrary and non-isotropic covariance models.
 
-Gallery
--------
+.. only:: html
+
+   Gallery
+   -------
+
+   Below is a gallery of examples

examples/02_cov_model/README.rst

@@ -8,32 +8,45 @@ class, which allows you to easily define arbitrary covariance models by
 yourself. The resulting models provide a bunch of nice features to explore the
 covariance models.
 
-
 A covariance model is used to characterize the
 `semi-variogram <https://en.wikipedia.org/wiki/Variogram#Semivariogram>`_,
 denoted by :math:`\gamma`, of a spatial random field.
 In GSTools, we use the following form for an isotropic and stationary field:
 
 .. math::
    \gamma\left(r\right)=
-   \sigma^2\cdot\left(1-\rho\left(r\right)\right)+n
+   \sigma^2\cdot\left(1-\mathrm{cor}\left(s\cdot\frac{r}{\ell}\right)\right)+n
 
 Where:
 
-  - :math:`\rho(r)` is the so called
-    `correlation <https://en.wikipedia.org/wiki/Autocovariance#Normalization>`_
-    function depending on the distance :math:`r`
+  - :math:`r` is the lag distance
+  - :math:`\ell` is the main correlation length
+  - :math:`s` is a scaling factor for unit conversion or normalization
   - :math:`\sigma^2` is the variance
   - :math:`n` is the nugget (subscale variance)
+  - :math:`\mathrm{cor}(h)` is the normalized correlation function depending on
+    the non-dimensional distance :math:`h=s\cdot\frac{r}{\ell}`
+
+Depending on the normalized correlation function, all covariance models in
+GSTools are providing the following functions:
+
+  - :math:`\rho(r)=\mathrm{cor}\left(s\cdot\frac{r}{\ell}\right)`
+    is the so called
+    `correlation <https://en.wikipedia.org/wiki/Autocovariance#Normalization>`_
+    function
+  - :math:`C(r)=\sigma^2\cdot\rho(r)` is the so called
+    `covariance <https://en.wikipedia.org/wiki/Covariance_function>`_
+    function, which gives the name for our GSTools class
 
 .. note::
 
    We are not limited to isotropic models. GSTools supports anisotropy ratios
    for length scales in orthogonal transversal directions like:
 
-   - :math:`x` (main direction)
-   - :math:`y` (1. transversal direction)
-   - :math:`z` (2. transversal direction)
+   - :math:`x_0` (main direction)
+   - :math:`x_1` (1. transversal direction)
+   - :math:`x_2` (2. transversal direction)
+   - ...
 
    These main directions can also be rotated.
    Just have a look at the corresponding examples.
@@ -54,14 +67,21 @@ The following standard covariance models are provided by GSTools
     Linear
     Circular
     Spherical
-    Intersection
+    HyperSpherical
+    SuperSpherical
+    JBessel
 
 As a special feature, we also provide truncated power law (TPL) covariance models
 
 .. autosummary::
     TPLGaussian
     TPLExponential
     TPLStable
+    TPLSimple
+
+.. only:: html
+
+   Gallery
+   -------
 
-Gallery
--------
+   Below is a gallery of examples

examples/03_variogram/04_directional_2d.py

@@ -1,9 +1,11 @@
 """
-Directional variogram estimation in 2D
---------------------------------------
+Directional variogram estimation and fitting in 2D
+--------------------------------------------------
 
 In this example, we demonstrate how to estimate a directional variogram by
 setting the direction angles in 2D.
+
+Afterwards we will fit a model to this estimated variogram and show the result.
 """
 import numpy as np
 import gstools as gs
@@ -20,29 +22,40 @@
 
 ###############################################################################
 # 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.
+# tolerance of 11.25 degree and a bandwith of 8.
 
-bins = range(0, 40, 3)
-bin_c, vario, cnt = gs.vario_estimate(
+bins = range(0, 40, 2)
+bin_center, dir_vario, counts = gs.vario_estimate(
     *((x, y), field, bins),
-    angles=(angle, angle + np.pi / 2),  # main dir. and transversal dir.
+    direction=gs.rotated_main_axes(dim=2, angles=angle),
     angles_tol=np.pi / 16,
-    bandwidth=1.0,
+    bandwidth=8,
     mesh_type="structured",
     return_counts=True,
 )
 
+###############################################################################
+# Afterwards we can use the estimated variogram to fit a model to it:
+
+print("Original:")
+print(model)
+model.fit_variogram(bin_center, dir_vario)
+print("Fitted:")
+print(model)
+
 ###############################################################################
 # 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.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.legend(loc="lower right")
 
+model.plot("vario_axis", axis=0, ax=ax1, x_max=40, label="fit on axis 0")
+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")
 

examples/03_variogram/05_directional_3d.py

@@ -1,9 +1,11 @@
 """
-Directional variogram estimation in 3D
---------------------------------------
+Directional variogram estimation and fitting in 3D
+--------------------------------------------------
 
 In this example, we demonstrate how to estimate a directional variogram by
 setting the estimation directions in 3D.
+
+Afterwards we will fit a model to this estimated variogram and show the result.
 """
 import numpy as np
 import gstools as gs
@@ -22,31 +24,45 @@
 field = srf.structured((x, y, z))
 
 ###############################################################################
-# Here we generate the rotated coordinate system to get an impression what
-# the rotation angles do.
+# Here we generate the axes of 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
+# All 3 axes of the rotated coordinate-system
+main_axes = gs.rotated_main_axes(dim, angles)
+axis1, axis2, axis3 = main_axes
 
 ###############################################################################
-# Now we estimate the variogram along the main axis. When the main axis is
+# Now we estimate the variogram along the main axes. When the main axes are
 # 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(
+bin_center, dir_vario, counts = gs.vario_estimate(
     *([x, y, z], field, bins),
-    direction=(dir0, dir1, dir2),
+    direction=main_axes,
     bandwidth=10,
     sampling_size=2000,
     sampling_seed=1001,
-    mesh_type="structured"
+    mesh_type="structured",
+    return_counts=True,
 )
 
+###############################################################################
+# 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)
+print("Fitted:")
+print(model)
+
 ###############################################################################
 # Plotting.
 
@@ -58,9 +74,9 @@
 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.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)
@@ -70,8 +86,13 @@
 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.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.set_title("Fitting an anisotropic model")
 ax3.legend()
+
 plt.show()

examples/03_variogram/README.rst

@@ -10,5 +10,9 @@ The same `(semi-)variogram <https://en.wikipedia.org/wiki/Variogram#Semivariogra
 :ref:`tutorial_02_cov` is being used
 by this subpackage.
 
-Gallery
--------
+.. only:: html
+
+   Gallery
+   -------
+
+   Below is a gallery of examples

examples/04_vector_field/README.rst

@@ -32,5 +32,9 @@ with the projector
 By calculating :math:`\nabla \cdot \mathbf U = 0`, it can be verified, that
 the resulting field is indeed incompressible.
 
-Gallery
--------
+.. only:: html
+
+   Gallery
+   -------
+
+   Below is a gallery of examples

examples/05_kriging/README.rst

@@ -42,7 +42,7 @@ the features you want:
   one can set `exact` to `False` and provide either individual measurement errors
   for each point or set the nugget as a constant measurement error everywhere.
 * `pseudo_inv`: Sometimes the inversion of the kriging matrix can be numerically unstable.
-  This occures for examples in cases of redundant input values. In this case we provide a switch to
+  This occurs for examples in cases of redundant input values. In this case we provide a switch to
   use the pseudo-inverse of the matrix. Then redundant conditional values will automatically
   be averaged.
 
@@ -87,5 +87,9 @@ submodule :any:`gstools.krige`.
     ExtDrift
     Detrended
 
-Gallery
--------
+.. only:: html
+
+   Gallery
+   -------
+
+   Below is a gallery of examples

examples/06_conditioned_fields/README.rst

@@ -30,5 +30,9 @@ or:
 
     srf.set_condition(cond_pos, cond_val, "ordinary")
 
-Gallery
--------
+.. only:: html
+
+   Gallery
+   -------
+
+   Below is a gallery of examples

examples/07_transformations/README.rst

@@ -33,5 +33,9 @@ Simply import the transform submodule and apply a transformation to the srf clas
     ...
     tf.normal_to_lognormal(srf)
 
-Gallery
--------
+.. only:: html
+
+   Gallery
+   -------
+
+   Below is a gallery of examples