Merge pull request #160 from GeoStat-Framework/pykrige_compat

Krige: Adapting to PyKrige

Author: MuellerSeb

examples/05_kriging/02_pykrige_interface.py

@@ -5,59 +5,58 @@
 To use fancier methods like
 `regression kriging <https://en.wikipedia.org/wiki/Regression-kriging>`__,
 we provide an interface to
-`PyKrige <https://github.com/bsmurphy/PyKrige>`__.
+`PyKrige <https://github.com/GeoStat-Framework/PyKrige>`__ (>v1.5), which means
+you can pass a GSTools covariance model to the kriging routines of PyKrige.
 
-In the future you can pass a GSTools Covariance Model
-to the PyKrige routines as ``variogram_model``.
-
-At the moment we only provide prepared
-keyword arguments for the pykrige routines.
-
-To demonstrate the general workflow, we compare the ordinary kriging of PyKrige
-with GSTools in 2D:
+To demonstrate the general workflow, we compare ordinary kriging of PyKrige
+with the corresponding GSTools routine in 2D:
 """
 import numpy as np
 import gstools as gs
 from pykrige.ok import OrdinaryKriging
 from matplotlib import pyplot as plt
 
 # conditioning data
-data = np.array(
-    [
-        [0.3, 1.2, 0.47],
-        [1.9, 0.6, 0.56],
-        [1.1, 3.2, 0.74],
-        [3.3, 4.4, 1.47],
-        [4.7, 3.8, 1.74],
-    ]
-)
+cond_x = [0.3, 1.9, 1.1, 3.3, 4.7]
+cond_y = [1.2, 0.6, 3.2, 4.4, 3.8]
+cond_val = [0.47, 0.56, 0.74, 1.47, 1.74]
 
 # grid definition for output field
 gridx = np.arange(0.0, 5.5, 0.1)
 gridy = np.arange(0.0, 6.5, 0.1)
 
 ###############################################################################
-# A GSTools based covariance model.
+# A GSTools based :any:`Gaussian` covariance model:
 
-cov_model = gs.Gaussian(
+model = gs.Gaussian(
     dim=2, len_scale=1, anis=0.2, angles=-0.5, var=0.5, nugget=0.1
 )
 
 ###############################################################################
-# Ordinary kriging with pykrige.
-# A dictionary containing keyword arguments for the pykrige routines is
-# provided by the gstools covariance models.
-
-pk_kwargs = cov_model.pykrige_kwargs
-OK1 = OrdinaryKriging(data[:, 0], data[:, 1], data[:, 2], **pk_kwargs)
+# Ordinary Kriging with PyKrige
+# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+#
+# One can pass the defined GSTools model as
+# variogram model, which will `not` be fitted to the given data.
+# By providing the GSTools model, rotation and anisotropy are also
+# automatically defined:
+
+OK1 = OrdinaryKriging(cond_x, cond_y, cond_val, variogram_model=model)
 z1, ss1 = OK1.execute("grid", gridx, gridy)
 plt.imshow(z1, origin="lower")
 plt.show()
 
 ###############################################################################
-# Ordinary kriging with gstools for comparison.
-
-OK2 = gs.krige.Ordinary(cov_model, [data[:, 0], data[:, 1]], data[:, 2])
+# Ordinary Kriging with GSTools
+# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+#
+# The :any:`Ordinary` kriging class is provided by GSTools as a shortcut to
+# define ordinary kriging with the general :any:`Krige` class.
+#
+# PyKrige's routines are using exact kriging by default (when given a nugget).
+# To reproduce this behavior in GSTools, we have to set ``exact=True``.
+
+OK2 = gs.krige.Ordinary(model, [cond_x, cond_y], cond_val, exact=True)
 OK2.structured([gridx, gridy])
 ax = OK2.plot()
 ax.set_aspect("equal")

gstools/krige/base.py

@@ -28,7 +28,7 @@
 __all__ = ["Krige"]
 
 
-P_INV = {1: spl.pinv, 2: spl.pinv2, 3: spl.pinvh}
+P_INV = {"pinv": spl.pinv, "pinv2": spl.pinv2, "pinvh": spl.pinvh}
 """dict: Standard pseudo-inverse routines"""
 
 
@@ -93,16 +93,16 @@ class Krige(Field):
         kriging matrix. If `True`, this leads to more numerical stability
         and redundant points are averaged. But it can take more time.
         Default: True
-    pseudo_inv_type : :class:`int` or :any:`callable`, optional
+    pseudo_inv_type : :class:`str` or :any:`callable`, optional
         Here you can select the algorithm to compute the pseudo-inverse matrix:
 
-            * `1`: use `pinv` from `scipy` which uses `lstsq`
-            * `2`: use `pinv2` from `scipy` which uses `SVD`
-            * `3`: use `pinvh` from `scipy` which uses eigen-values
+            * `"pinv"`: use `pinv` from `scipy` which uses `lstsq`
+            * `"pinv2"`: use `pinv2` from `scipy` which uses `SVD`
+            * `"pinvh"`: use `pinvh` from `scipy` which uses eigen-values
 
         If you want to use another routine to invert the kriging matrix,
         you can pass a callable which takes a matrix and returns the inverse.
-        Default: `1`
+        Default: `"pinv"`
     fit_normalizer : :class:`bool`, optional
         Wheater to fit the data-normalizer to the given conditioning data.
         Default: False
@@ -139,7 +139,7 @@ def __init__(
         exact=False,
         cond_err="nugget",
         pseudo_inv=True,
-        pseudo_inv_type=1,
+        pseudo_inv_type="pinv",
         fit_normalizer=False,
         fit_variogram=False,
     ):
@@ -650,13 +650,13 @@ def pseudo_inv(self):
 
     @property
     def pseudo_inv_type(self):
-        """:class:`int`: Method selector for pseudo inverse calculation."""
+        """:class:`str`: Method selector for pseudo inverse calculation."""
         return self._pseudo_inv_type
 
     @pseudo_inv_type.setter
     def pseudo_inv_type(self, val):
-        if val not in [1, 2, 3] and not callable(val):
-            raise ValueError("Krige: pseudo_inv_type needs to be in [1,2,3]")
+        if val not in P_INV and not callable(val):
+            raise ValueError(f"Krige: pseudo_inv_type not in {sorted(P_INV)}")
         self._pseudo_inv_type = val
 
     @property

gstools/krige/methods.py

@@ -64,16 +64,16 @@ class Simple(Krige):
         kriging matrix. If `True`, this leads to more numerical stability
         and redundant points are averaged. But it can take more time.
         Default: True
-    pseudo_inv_type : :class:`int` or :any:`callable`, optional
+    pseudo_inv_type : :class:`str` or :any:`callable`, optional
         Here you can select the algorithm to compute the pseudo-inverse matrix:
 
-            * `1`: use `pinv` from `scipy` which uses `lstsq`
-            * `2`: use `pinv2` from `scipy` which uses `SVD`
-            * `3`: use `pinvh` from `scipy` which uses eigen-values
+            * `"pinv"`: use `pinv` from `scipy` which uses `lstsq`
+            * `"pinv2"`: use `pinv2` from `scipy` which uses `SVD`
+            * `"pinvh"`: use `pinvh` from `scipy` which uses eigen-values
 
         If you want to use another routine to invert the kriging matrix,
         you can pass a callable which takes a matrix and returns the inverse.
-        Default: `1`
+        Default: `"pinv"`
     fit_normalizer : :class:`bool`, optional
         Wheater to fit the data-normalizer to the given conditioning data.
         Default: False
@@ -96,7 +96,7 @@ def __init__(
         exact=False,
         cond_err="nugget",
         pseudo_inv=True,
-        pseudo_inv_type=1,
+        pseudo_inv_type="pinv",
         fit_normalizer=False,
         fit_variogram=False,
     ):
@@ -159,16 +159,16 @@ class Ordinary(Krige):
         kriging matrix. If `True`, this leads to more numerical stability
         and redundant points are averaged. But it can take more time.
         Default: True
-    pseudo_inv_type : :class:`int` or :any:`callable`, optional
+    pseudo_inv_type : :class:`str` or :any:`callable`, optional
         Here you can select the algorithm to compute the pseudo-inverse matrix:
 
-            * `1`: use `pinv` from `scipy` which uses `lstsq`
-            * `2`: use `pinv2` from `scipy` which uses `SVD`
-            * `3`: use `pinvh` from `scipy` which uses eigen-values
+            * `"pinv"`: use `pinv` from `scipy` which uses `lstsq`
+            * `"pinv2"`: use `pinv2` from `scipy` which uses `SVD`
+            * `"pinvh"`: use `pinvh` from `scipy` which uses eigen-values
 
         If you want to use another routine to invert the kriging matrix,
         you can pass a callable which takes a matrix and returns the inverse.
-        Default: `1`
+        Default: `"pinv"`
     fit_normalizer : :class:`bool`, optional
         Wheater to fit the data-normalizer to the given conditioning data.
         Default: False
@@ -190,7 +190,7 @@ def __init__(
         exact=False,
         cond_err="nugget",
         pseudo_inv=True,
-        pseudo_inv_type=1,
+        pseudo_inv_type="pinv",
         fit_normalizer=False,
         fit_variogram=False,
     ):
@@ -264,16 +264,16 @@ class Universal(Krige):
         kriging matrix. If `True`, this leads to more numerical stability
         and redundant points are averaged. But it can take more time.
         Default: True
-    pseudo_inv_type : :class:`int` or :any:`callable`, optional
+    pseudo_inv_type : :class:`str` or :any:`callable`, optional
         Here you can select the algorithm to compute the pseudo-inverse matrix:
 
-            * `1`: use `pinv` from `scipy` which uses `lstsq`
-            * `2`: use `pinv2` from `scipy` which uses `SVD`
-            * `3`: use `pinvh` from `scipy` which uses eigen-values
+            * `"pinv"`: use `pinv` from `scipy` which uses `lstsq`
+            * `"pinv2"`: use `pinv2` from `scipy` which uses `SVD`
+            * `"pinvh"`: use `pinvh` from `scipy` which uses eigen-values
 
         If you want to use another routine to invert the kriging matrix,
         you can pass a callable which takes a matrix and returns the inverse.
-        Default: `1`
+        Default: `"pinv"`
     fit_normalizer : :class:`bool`, optional
         Wheater to fit the data-normalizer to the given conditioning data.
         Default: False
@@ -296,7 +296,7 @@ def __init__(
         exact=False,
         cond_err="nugget",
         pseudo_inv=True,
-        pseudo_inv_type=1,
+        pseudo_inv_type="pinv",
         fit_normalizer=False,
         fit_variogram=False,
     ):
@@ -366,16 +366,16 @@ class ExtDrift(Krige):
         kriging matrix. If `True`, this leads to more numerical stability
         and redundant points are averaged. But it can take more time.
         Default: True
-    pseudo_inv_type : :class:`int` or :any:`callable`, optional
+    pseudo_inv_type : :class:`str` or :any:`callable`, optional
         Here you can select the algorithm to compute the pseudo-inverse matrix:
 
-            * `1`: use `pinv` from `scipy` which uses `lstsq`
-            * `2`: use `pinv2` from `scipy` which uses `SVD`
-            * `3`: use `pinvh` from `scipy` which uses eigen-values
+            * `"pinv"`: use `pinv` from `scipy` which uses `lstsq`
+            * `"pinv2"`: use `pinv2` from `scipy` which uses `SVD`
+            * `"pinvh"`: use `pinvh` from `scipy` which uses eigen-values
 
         If you want to use another routine to invert the kriging matrix,
         you can pass a callable which takes a matrix and returns the inverse.
-        Default: `1`
+        Default: `"pinv"`
     fit_normalizer : :class:`bool`, optional
         Wheater to fit the data-normalizer to the given conditioning data.
         Default: False
@@ -398,7 +398,7 @@ def __init__(
         exact=False,
         cond_err="nugget",
         pseudo_inv=True,
-        pseudo_inv_type=1,
+        pseudo_inv_type="pinv",
         fit_normalizer=False,
         fit_variogram=False,
     ):
@@ -461,16 +461,16 @@ class Detrended(Krige):
         kriging matrix. If `True`, this leads to more numerical stability
         and redundant points are averaged. But it can take more time.
         Default: True
-    pseudo_inv_type : :class:`int` or :any:`callable`, optional
+    pseudo_inv_type : :class:`str` or :any:`callable`, optional
         Here you can select the algorithm to compute the pseudo-inverse matrix:
 
-            * `1`: use `pinv` from `scipy` which uses `lstsq`
-            * `2`: use `pinv2` from `scipy` which uses `SVD`
-            * `3`: use `pinvh` from `scipy` which uses eigen-values
+            * `"pinv"`: use `pinv` from `scipy` which uses `lstsq`
+            * `"pinv2"`: use `pinv2` from `scipy` which uses `SVD`
+            * `"pinvh"`: use `pinvh` from `scipy` which uses eigen-values
 
         If you want to use another routine to invert the kriging matrix,
         you can pass a callable which takes a matrix and returns the inverse.
-        Default: `1`
+        Default: `"pinv"`
     fit_variogram : :class:`bool`, optional
         Wheater to fit the given variogram model to the data.
         This is done by using isotropy settings of the given model,
@@ -488,7 +488,7 @@ def __init__(
         exact=False,
         cond_err="nugget",
         pseudo_inv=True,
-        pseudo_inv_type=1,
+        pseudo_inv_type="pinv",
         fit_variogram=False,
     ):
         super().__init__(