Krige, Field and Covmodel update

CHANGELOG.md

@@ -7,9 +7,20 @@ All notable changes to **GSTools** will be documented in this file.
 
 ### Enhancements
 - different variogram estimator functions can now be used #51
+- the TPLGaussian and TPLExponential now have analytical spectra #67
+- added property ``is_isotropic`` to CovModel #67
+- reworked the whole krige sub-module to provide multiple kriging methods #67
+  - Simple
+  - Ordinary
+  - Universal
+  - External Drift Kriging
+  - Detrended Kriging
 
 ### Changes
 - Python versions 2.7 and 3.4 are no longer supported #40 #43
+- CovModel: in 3D the input of anisotropy is now treated slightly different: #67
+  - single given anisotropy value [e] is converted to [1, e] (it was [e, e] before)
+  - two given length-scales [l_1, l_2] are converted to [l_1, l_2, l_2] (it was [l_1, l_2, l_1] before)
 
 ### Bugfixes
 - a race condition in the structured variogram estimation has been fixed #51

docs/source/conf.py

@@ -178,7 +178,9 @@ def setup(app):
 
 # Output file base name for HTML help builder.
 htmlhelp_basename = "GeoStatToolsdoc"
-
+# logos for the page
+html_logo = "pics/gstools_150.png"
+html_favicon = "pics/gstools.ico"
 
 # -- Options for LaTeX output ---------------------------------------------
 # latex_engine = 'lualatex'

docs/source/krige.rst

@@ -2,10 +2,6 @@ gstools.krige
 =============
 
 .. automodule:: gstools.krige
-   :members:
-   :undoc-members:
-   :inherited-members:
-   :show-inheritance:
 
 .. raw:: latex
 

examples/02_cov_model/01_basic_methods.py

@@ -8,24 +8,24 @@
 One of the following functions defines the main characterization of the
 variogram:
 
-- :any:`CovModel.variogram` : The variogram of the model given by
+- ``CovModel.variogram`` : The variogram of the model given by
 
   .. math::
       \gamma\left(r\right)=
       \sigma^2\cdot\left(1-\rho\left(r\right)\right)+n
 
-- :any:`CovModel.covariance` : The (auto-)covariance of the model given by
+- ``CovModel.covariance`` : The (auto-)covariance of the model given by
 
   .. math::
       C\left(r\right)= \sigma^2\cdot\rho\left(r\right)
 
-- :any:`CovModel.correlation` : The (auto-)correlation
+- ``CovModel.correlation`` : The (auto-)correlation
   (or normalized covariance) of the model given by
 
   .. math::
       \rho\left(r\right)
 
-- :any:`CovModel.cor` : The normalized correlation taking a
+- ``CovModel.cor`` : The normalized correlation taking a
   normalized range given by:
 
   .. math::

examples/05_kriging/04_extdrift_kriging.py

@@ -0,0 +1,24 @@
+"""
+External Drift Kriging
+----------------------
+"""
+import numpy as np
+from gstools import SRF, Gaussian, krige
+
+# synthetic condtions with a drift
+drift_model = Gaussian(dim=1, len_scale=4)
+drift = SRF(drift_model, seed=1010)
+cond_pos = [0.3, 1.9, 1.1, 3.3, 4.7]
+ext_drift = drift(cond_pos)
+cond_val = ext_drift * 2 + 1
+# resulting grid
+gridx = np.linspace(0.0, 15.0, 151)
+grid_drift = drift(gridx)
+# kriging
+model = Gaussian(dim=1, var=2, len_scale=4)
+krig = krige.ExtDrift(model, cond_pos, cond_val, ext_drift)
+krig(gridx, ext_drift=grid_drift)
+ax = krig.plot()
+ax.scatter(cond_pos, cond_val, color="k", zorder=10, label="Conditions")
+ax.plot(gridx, grid_drift, label="drift")
+ax.legend()

examples/05_kriging/05_universal_kriging.py

@@ -0,0 +1,24 @@
+"""
+Universal Kriging
+-----------------
+"""
+import numpy as np
+from gstools import SRF, Gaussian, krige
+
+# synthetic condtions with a drift
+drift_model = Gaussian(dim=1, var=0.1, len_scale=2)
+drift = SRF(drift_model, seed=101)
+cond_pos = np.linspace(0.1, 8, 10)
+cond_val = drift(cond_pos) + cond_pos * 0.1 + 1
+# resulting grid
+gridx = np.linspace(0.0, 15.0, 151)
+drift_field = drift(gridx) + gridx * 0.1 + 1
+# kriging
+model = Gaussian(dim=1, var=0.1, len_scale=2)
+krig = krige.Universal(model, cond_pos, cond_val, "linear")
+krig(gridx)
+ax = krig.plot()
+ax.scatter(cond_pos, cond_val, color="k", zorder=10, label="Conditions")
+ax.plot(gridx, gridx * 0.1 + 1, ":", label="linear drift")
+ax.plot(gridx, drift_field, "--", label="original field")
+ax.legend()

examples/05_kriging/06_detrended_kriging.py

@@ -0,0 +1,30 @@
+"""
+Detrended Kriging
+-----------------
+"""
+import numpy as np
+from gstools import SRF, Gaussian, krige
+
+
+def trend(x):
+    """Example for a simple linear trend."""
+    return x * 0.1 + 1
+
+
+# synthetic condtions with trend/drift
+drift_model = Gaussian(dim=1, var=0.1, len_scale=2)
+drift = SRF(drift_model, seed=101)
+cond_pos = np.linspace(0.1, 8, 10)
+cond_val = drift(cond_pos) + trend(cond_pos)
+# resulting grid
+gridx = np.linspace(0.0, 15.0, 151)
+drift_field = drift(gridx) + trend(gridx)
+# kriging
+model = Gaussian(dim=1, var=0.1, len_scale=2)
+krig_trend = krige.Detrended(model, cond_pos, cond_val, trend)
+krig_trend(gridx)
+ax = krig_trend.plot()
+ax.scatter(cond_pos, cond_val, color="k", zorder=10, label="Conditions")
+ax.plot(gridx, trend(gridx), ":", label="linear trend")
+ax.plot(gridx, drift_field, "--", label="original field")
+ax.legend()

examples/05_kriging/07_detrended_ordinary_kriging.py

@@ -0,0 +1,30 @@
+"""
+Detrended Ordinary Kriging
+--------------------------
+"""
+import numpy as np
+from gstools import SRF, Gaussian, krige
+
+
+def trend(x):
+    """Example for a simple linear trend."""
+    return x * 0.1 + 1
+
+
+# synthetic condtions with trend/drift
+drift_model = Gaussian(dim=1, var=0.1, len_scale=2)
+drift = SRF(drift_model, seed=101)
+cond_pos = np.linspace(0.1, 8, 10)
+cond_val = drift(cond_pos) + trend(cond_pos)
+# resulting grid
+gridx = np.linspace(0.0, 15.0, 151)
+drift_field = drift(gridx) + trend(gridx)
+# kriging
+model = Gaussian(dim=1, var=0.1, len_scale=2)
+krig_trend = krige.Ordinary(model, cond_pos, cond_val, trend)
+krig_trend(gridx)
+ax = krig_trend.plot()
+ax.scatter(cond_pos, cond_val, color="k", zorder=10, label="Conditions")
+ax.plot(gridx, trend(gridx), ":", label="linear trend")
+ax.plot(gridx, drift_field, "--", label="original field")
+ax.legend()

examples/05_kriging/README.rst

@@ -6,12 +6,6 @@ Tutorial 5: Kriging
 The subpackage :py:mod:`gstools.krige` provides routines for Gaussian process regression, also known as kriging.
 Kriging is a method of data interpolation based on predefined covariance models.
 
-We provide two kinds of kriging routines:
-
-* Simple: The data is interpolated with a given mean value for the kriging field.
-* Ordinary: The mean of the resulting field is unkown and estimated during interpolation.
-
-
 The aim of kriging is to derive the value of a field at some point :math:`x_0`,
 when there are fixed observed values :math:`z(x_1)\ldots z(x_n)` at given points :math:`x_i`.
 
@@ -26,24 +20,38 @@ The weights :math:`W = (w_1,\ldots,w_n)` depent on the given covariance model an
 The different kriging approaches provide different ways of calculating :math:`W`.
 
 
-
 The routines for kriging are almost identical to the routines for spatial random fields.
 First you define a covariance model, as described in :ref:`tutorial_02_cov`,
 then you initialize the kriging class with this model:
 
 .. code-block:: python
 
-    from gstools import Gaussian, krige
+    import gstools as gs
     # condtions
-    cond_pos = ...
-    cond_val = ...
-    model = Gaussian(dim=1, var=0.5, len_scale=2)
-    krig = krige.Simple(model, mean=1, cond_pos=cond_pos, cond_val=cond_val)
+    cond_pos = [...]
+    cond_val = [...]
+    model = gs.Gaussian(dim=1, var=0.5, len_scale=2)
+    krig = gs.krige.Simple(model, cond_pos=cond_pos, cond_val=cond_val, mean=1)
 
-The resulting field instance ``krig`` has the same methods as the :any:`SRF` class.
+The resulting field instance ``krig`` has the same methods as the
+:any:`SRF` class.
 You can call it to evaluate the kriged field at different points,
 you can plot the latest field or you can export the field and so on.
-Have a look at the documentation of :any:`Simple` and :any:`Ordinary`.
+
+Provided Kriging Methods
+------------------------
+
+.. currentmodule:: gstools.krige
+
+The following kriging methods are provided within the
+submodule :any:`gstools.krige`.
+
+.. autosummary::
+    Simple
+    Ordinary
+    Universal
+    ExtDrift
+    Detrended
 
 Gallery
 -------

gstools/_version.py

@@ -1,2 +1,2 @@
 """Provide a central version."""
-__version__ = "1.1.2.dev0"
+__version__ = "1.2.0.dev0"

gstools/covmodel/base.py

@@ -2,7 +2,7 @@
 """
 GStools subpackage providing the base class for covariance models.
 
-.. currentmodule:: gstools.covmodel
+.. currentmodule:: gstools.covmodel.base
 
 The following classes are provided
 
@@ -61,15 +61,22 @@ class CovModel(metaclass=InitSubclassMeta):
         If a single value is given, the same length-scale will be used for
         every direction. If multiple values (for main and transversal
         directions) are given, `anis` will be
-        recalculated accordingly.
+        recalculated accordingly. If only two values are given in 3D,
+        the latter one will be used for both transversal directions.
         Default: ``1.0``
     nugget : :class:`float`, optional
         nugget of the model. Default: ``0.0``
     anis : :class:`float` or :class:`list`, optional
-        anisotropy ratios in the transversal directions [y, z].
+        anisotropy ratios in the transversal directions [e_y, e_z].
+
+            * e_y = l_y / l_x
+            * e_z = l_z / l_x
+
+        If only one value is given in 3D, e_y will be set to 1.
+        This value will be ignored, if multiple len_scales are given.
         Default: ``1.0``
     angles : :class:`float` or :class:`list`, optional
-        angles of rotation:
+        angles of rotation (given in rad):
 
             * in 2D: given as rotation around z-axis
             * in 3D: given by yaw, pitch, and roll (known as Tait–Bryan angles)
@@ -218,7 +225,7 @@ def __init_subclass__(cls):
         Best practice is to use the ``correlation`` function, or the ``cor``
         function. The latter one takes the dimensionles distance h=r/l.
         """
-        # overrid one of these ################################################
+        # override one of these ###############################################
 
         def variogram(self, r):
             r"""Isotropic variogram of the model.
@@ -326,7 +333,8 @@ def _get_iso_rad(self, pos):
         x, y, z = pos2xyz(pos, max_dim=self.dim)
         if self.do_rotation:
             x, y, z = unrotate_mesh(self.dim, self.angles, x, y, z)
-        y, z = make_isotropic(self.dim, self.anis, y, z)
+        if not self.is_isotropic:
+            y, z = make_isotropic(self.dim, self.anis, y, z)
         return np.linalg.norm((x, y, z)[: self.dim], axis=0)
 
     def vario_spatial(self, pos):
@@ -638,11 +646,8 @@ def spectral_rad_pdf(self, r):
 
     def ln_spectral_rad_pdf(self, r):
         """Log radial spectral density of the model."""
-        spec = np.array(self.spectral_rad_pdf(r))
-        spec_gz = np.logical_not(np.isclose(spec, 0))
-        res = np.full_like(spec, -np.inf, dtype=np.double)
-        res[spec_gz] = np.log(spec[spec_gz])
-        return res
+        with np.errstate(divide="ignore"):
+            return np.log(self.spectral_rad_pdf(r))
 
     def _has_cdf(self):
         """State if a cdf is defined with 'spectral_rad_cdf'."""
@@ -866,7 +871,7 @@ def check_arg_bounds(self):
                         + str(val)
                     )
 
-    ### bounds  properties ####################################################
+    ### bounds properties #####################################################
 
     @property
     def var_bounds(self):
@@ -1121,8 +1126,6 @@ def hankel_kw(self, hankel_kw):
         if self.dim is not None:
             self._sft = SFT(ndim=self.dim, **self.hankel_kw)
 
-    ### properties ############################################################
-
     @property
     def dist_func(self):
         """:class:`tuple` of :any:`callable`: pdf, cdf and ppf.
@@ -1211,9 +1214,14 @@ 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)) and not self.is_isotropic
+        )
 
-    ### magic methods #########################################################
+    @property
+    def is_isotropic(self):
+        """:any:`bool`: State if a model is isotropic."""
+        return np.all(np.isclose(self.anis, 1.0))
 
     def __eq__(self, other):
         """Compare CovModels."""

gstools/covmodel/models.py

@@ -2,7 +2,7 @@
 """
 GStools subpackage providing different covariance models.
 
-.. currentmodule:: gstools.covmodel
+.. currentmodule:: gstools.covmodel.models
 
 The following classes and functions are provided
 
@@ -100,9 +100,7 @@ def spectral_rad_ppf(self, u):
 
     def _has_ppf(self):
         # since the ppf is not analytical for dim=3, we have to state that
-        if self.dim == 3:
-            return False
-        return True
+        return False if self.dim == 3 else True
 
     def calc_integral_scale(self):  # noqa: D102
         return self.len_scale
@@ -183,9 +181,7 @@ def spectral_rad_ppf(self, u):
 
     def _has_ppf(self):
         # since the ppf is not analytical for dim=3, we have to state that
-        if self.dim == 3:
-            return False
-        return True
+        return False if self.dim == 3 else True
 
     def calc_integral_scale(self):  # noqa: D102
         return self.len_scale