Krige unification

examples/05_kriging/01_ordinary_kriging.py

@@ -10,14 +10,14 @@
 .. math::
 
    \begin{pmatrix}W\\\mu\end{pmatrix} = \begin{pmatrix}
-   \gamma(x_1,x_1) & \cdots & \gamma(x_1,x_n) &1 \\
+   c(x_1,x_1) & \cdots & c(x_1,x_n) &1 \\
    \vdots & \ddots & \vdots  & \vdots \\
-   \gamma(x_n,x_1) & \cdots & \gamma(x_n,x_n) & 1 \\
+   c(x_n,x_1) & \cdots & c(x_n,x_n) & 1 \\
    1 &\cdots& 1 & 0
    \end{pmatrix}^{-1}
-   \begin{pmatrix}\gamma(x_1,x_0) \\ \vdots \\ \gamma(x_n,x_0) \\ 1\end{pmatrix}
+   \begin{pmatrix}c(x_1,x_0) \\ \vdots \\ c(x_n,x_0) \\ 1\end{pmatrix}
 
-Thereby :math:`\gamma(x_i,x_j)` is the semi-variogram of the given observations
+Thereby :math:`c(x_i,x_j)` is the covariance of the given observations
 and :math:`\mu` is a Lagrange multiplier to minimize the kriging error and estimate the mean.
 
 

examples/05_kriging/05_universal_kriging.py

@@ -1,6 +1,17 @@
 """
 Universal Kriging
 -----------------
+
+You can give a polynomial order or a list of self defined
+functions representing the internal drift of the given values.
+This drift will be fitted internally during the kriging interpolation.
+
+In the following we are creating artificial data, where a linear drift
+was added. The resulting samples are then used as input for Universal kriging.
+
+The "linear" drift is then estimated during the interpolation.
+To access only the estimated mean/drift, we provide a switch `only_mean`
+in the call routine.
 """
 import numpy as np
 from gstools import SRF, Gaussian, krige
@@ -21,4 +32,8 @@
 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")
+
+mean, mean_err = krig(gridx, only_mean=True)
+ax.plot(gridx, mean, label="estimated drift")
+
 ax.legend()

examples/05_kriging/08_measurement_errors.py

@@ -0,0 +1,55 @@
+r"""
+Incorporating measurement errors
+--------------------------------
+
+To incorporate the nugget effect and/or given measurement errors,
+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.
+
+In the following we will show the influence of the nugget and
+measurement errors.
+"""
+
+import numpy as np
+import gstools as gs
+
+# condtions
+cond_pos = [0.3, 1.1, 1.9, 3.3, 4.7]
+cond_val = [0.47, 0.74, 0.56, 1.47, 1.74]
+cond_err = [0.01, 0.0, 0.1, 0.05, 0]
+# resulting grid
+gridx = np.linspace(0.0, 15.0, 151)
+# spatial random field class
+model = gs.Gaussian(dim=1, var=0.9, len_scale=1, nugget=0.1)
+
+###############################################################################
+# Here we will use Simple kriging (`unbiased=False`) to interpolate the given
+# conditions.
+
+krig = gs.krige.Krige(
+    model=model,
+    cond_pos=cond_pos,
+    cond_val=cond_val,
+    mean=1,
+    unbiased=False,
+    exact=False,
+    cond_err=cond_err,
+)
+krig(gridx)
+
+###############################################################################
+# Let's plot the data. You can see, that the estimated values differ more from
+# the input, when the given measurement errors get bigger.
+# In addition we plot the standard deviation.
+
+ax = krig.plot()
+ax.scatter(cond_pos, cond_val, color="k", zorder=10, label="Conditions")
+ax.fill_between(
+    gridx,
+    # plus/minus standard deviation (70 percent confidence interval)
+    krig.field - np.sqrt(krig.krige_var),
+    krig.field + np.sqrt(krig.krige_var),
+    alpha=0.3,
+    label="Standard deviation",
+)
+ax.legend()

examples/05_kriging/09_pseudo_inverse.py

@@ -0,0 +1,38 @@
+r"""
+Redundant data and pseudo-inverse
+---------------------------------
+
+It can happen, that the kriging system gets numerically unstable.
+One reason could be, that the input data contains redundant conditioning points
+that hold different values.
+
+To smoothly deal with such situations, you can switch on using the pseuodo
+inverse for the kriging matrix.
+
+This will result in the average value for the redundant data.
+
+Example
+^^^^^^^
+
+In the following we have two different values at the same location.
+The resulting kriging field will hold the average at this point.
+"""
+import numpy as np
+from gstools import Gaussian, krige
+
+# condtions
+cond_pos = [0.3, 1.9, 1.1, 3.3, 1.1]
+cond_val = [0.47, 0.56, 0.74, 1.47, 1.14]
+# resulting grid
+gridx = np.linspace(0.0, 8.0, 81)
+# spatial random field class
+model = Gaussian(dim=1, var=0.5, len_scale=1)
+
+###############################################################################
+krig = krige.Ordinary(model, cond_pos=cond_pos, cond_val=cond_val)
+krig(gridx)
+
+###############################################################################
+ax = krig.plot()
+ax.scatter(cond_pos, cond_val, color="k", zorder=10, label="Conditions")
+ax.legend()

examples/05_kriging/README.rst

@@ -3,7 +3,8 @@
 Tutorial 5: Kriging
 ===================
 
-The subpackage :py:mod:`gstools.krige` provides routines for Gaussian process regression, also known as 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.
 
 The aim of kriging is to derive the value of a field at some point :math:`x_0`,
@@ -19,8 +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.
+The :any:`Krige` class provides everything in one place and you can switch on/off
+the features you want:
+
+* `unbiased`: the weights have to sum up to `1`. If true, this results in
+  :any:`Ordinary` kriging, where the mean is estimated, otherwise it will result in
+  :any:`Simple` kriging, where the mean has to be given.
+* `drift_functions`: you can give a polynomial order or a list of self defined
+  functions representing the internal drift of the given values. This drift will
+  be fitted internally during the kriging interpolation. This results in :any:`Universal` kriging.
+* `ext_drift`: You can also give an external drift per point to the routine.
+  In contrast to the internal drift, that is evaluated at the desired points with
+  the given functions, the external drift has to given for each point form an "external"
+  source. This results in :any:`ExtDrift` kriging.
+* `trend_function`: If you already have fitted a trend model, that is provided as a
+  callable, you can give it to the kriging routine. This trend is subtracted from the
+  conditional values before the kriging is done, meaning, that only the residuals are
+  used for kriging. This can be used with separate regression of your data.
+  This results in :any:`Detrended` kriging.
+* `exact` and `cond_err`: To incorporate the nugget effect and/or measurement errors,
+  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
+  use the pseudo-inverse of the matrix. Then redundant conditional values will automatically
+  be averaged.
+
+All mentioned features can be combined within the :any:`Krige` class.
+All other kriging classes are just shortcuts to this class with a limited list
+of input parameters.
+
+The routines for kriging are almost identical to the routines for spatial random fields,
+with regard to their handling.
 First you define a covariance model, as described in :ref:`tutorial_02_cov`,
 then you initialize the kriging class with this model:
 
@@ -47,6 +78,7 @@ The following kriging methods are provided within the
 submodule :any:`gstools.krige`.
 
 .. autosummary::
+    Krige
     Simple
     Ordinary
     Universal

examples/06_conditioned_fields/00_condition_ensemble.py

@@ -18,7 +18,7 @@
 ###############################################################################
 
 # spatial random field class
-model = gs.Gaussian(dim=1, var=0.5, len_scale=2)
+model = gs.Gaussian(dim=1, var=0.5, len_scale=1.5)
 srf = gs.SRF(model)
 srf.set_condition(cond_pos, cond_val, "ordinary")
 
@@ -34,6 +34,15 @@
 plt.plot(gridx, np.mean(fields, axis=0), linestyle=":", label="Ensemble mean")
 plt.plot(gridx, srf.krige_field, linestyle="dashed", label="kriged field")
 plt.scatter(cond_pos, cond_val, color="k", zorder=10, label="Conditions")
+# 99 percent confidence interval
+conf = gs.tools.confidence_scaling(0.99)
+plt.fill_between(
+    gridx,
+    srf.krige_field - conf * np.sqrt(srf.krige_var),
+    srf.krige_field + conf * np.sqrt(srf.krige_var),
+    alpha=0.3,
+    label="99% confidence interval",
+)
 plt.legend()
 plt.show()
 

gstools/field/base.py

@@ -185,7 +185,7 @@ def _pre_pos(self, pos, mesh_type="unstructured", make_unstruct=False):
         axis_lens : :class:`tuple` or :any:`None`
             axis lengths of the structured mesh if mesh type was changed.
         """
-        x, y, z = pos2xyz(pos, max_dim=self.model.dim)
+        x, y, z = pos2xyz(pos, dtype=np.double, max_dim=self.model.dim)
         pos = xyz2pos(x, y, z)
         mesh_type_gen = mesh_type
         # format the positional arguments of the mesh

gstools/field/condition.py

@@ -54,7 +54,7 @@ def ordinary(srf):
     )
     err_field, __ = err_ok(srf.pos, srf.mesh_type)
     cond_field = srf.raw_field + krige_field - err_field
-    info = {"mean": krige_ok.mean}
+    info = {"mean": krige_ok.get_mean()}
     return cond_field, krige_field, err_field, krige_var, info
 
 

gstools/krige/__init__.py

@@ -10,6 +10,7 @@
 .. autosummary::
    :toctree: generated
 
+   Krige
    Simple
    Ordinary
    Universal
@@ -18,6 +19,7 @@
 
 ----
 """
+from gstools.krige.base import Krige
 from gstools.krige.methods import (
     Simple,
     Ordinary,
@@ -26,4 +28,4 @@
     Detrended,
 )
 
-__all__ = ["Simple", "Ordinary", "Universal", "ExtDrift", "Detrended"]
+__all__ = ["Krige", "Simple", "Ordinary", "Universal", "ExtDrift", "Detrended"]