Feature: GSTools covariance model as variogram_model input (#125)

* added GSTools CovModel as possible variogram_model input; added example

* README GSTools example

* raise ValueError for wrong dim; remove links in comments

* add test for gstools interface

* update README gstools example

* update gstools example; remove png

* README: remove link to GSTools before PyKrige Documentation

* GSTools example: apply CI hack for plotting

* GSTools: update link in DOC

* CI: add gstools as dep

* GSTools: fix example

* CI: gstools test only on travis and py36

* CI: skip additional installation if sklearn and gstools

Author: MuellerSeb

.travis.yml

@@ -5,10 +5,10 @@ services:
 matrix:
     include:
     - python: "3.6"
-      env: DEPS="numpy scipy cython nose matplotlib scikit-learn"
-    - python: "3.5"                                                                                         
+      env: DEPS="numpy scipy cython nose matplotlib scikit-learn gstools"
+    - python: "3.5"
       env: DEPS="numpy scipy=0.18 cython nose matplotlib scikit-learn=0.18.1"
-    - python: "2.7"                                                                                         
+    - python: "2.7"
       env: DEPS="numpy=1.10.4 scipy=0.17 cython nose matplotlib scikit-learn=0.17.1"
 
 # setup adapted from https://github.com/soft-matter/trackpy/blob/master/.travis.yml
@@ -32,13 +32,13 @@ before_script:
   - "export DISPLAY=:99.0"
 
 # command to run tests
-script: 
+script:
     - source activate krige-env
     - pytest -sv ./pykrige/ --cov=pykrige
     - |
         # run examples
         # scikit-learn needed by the regression kriging example
-        conda install -y scikit-learn
+        # conda install -y scikit-learn gstools
         set -x
         for f in examples/*.py; do
            python $f 2>&1 | tee -a ~/log.txt
@@ -47,4 +47,3 @@ script:
 
 #after_success:
 #    coveralls
-

README.rst

@@ -62,34 +62,35 @@ First we will create a 2D dataset together with the associated x, y grids,
     import numpy as np
     import pykrige.kriging_tools as kt
     from pykrige.ok import OrdinaryKriging
-    
+
     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]])
-    
+
     gridx = np.arange(0.0, 5.5, 0.5)
     gridy = np.arange(0.0, 5.5, 0.5)
-    
+
     # Create the ordinary kriging object. Required inputs are the X-coordinates of
     # the data points, the Y-coordinates of the data points, and the Z-values of the
     # data points. If no variogram model is specified, defaults to a linear variogram
     # model. If no variogram model parameters are specified, then the code automatically
-    # calculates the parameters by fitting the variogram model to the binned 
+    # calculates the parameters by fitting the variogram model to the binned
     # experimental semivariogram. The verbose kwarg controls code talk-back, and
     # the enable_plotting kwarg controls the display of the semivariogram.
     OK = OrdinaryKriging(data[:, 0], data[:, 1], data[:, 2], variogram_model='linear',
                          verbose=False, enable_plotting=False)
-    					 
+
     # Creates the kriged grid and the variance grid. Allows for kriging on a rectangular
     # grid of points, on a masked rectangular grid of points, or with arbitrary points.
     # (See OrdinaryKriging.__doc__ for more information.)
     z, ss = OK.execute('grid', gridx, gridy)
-    
+
     # Writes the kriged grid to an ASCII grid file.
     kt.write_asc_grid(gridx, gridy, z, filename="output.asc")
 
+
 Universal Kriging Example
 ^^^^^^^^^^^^^^^^^^^^^^^^^
 
@@ -111,16 +112,17 @@ Universal Kriging Example
     # the data points, the Y-coordinates of the data points, and the Z-values of the
     # data points. Variogram is handled as in the ordinary kriging case.
     # drift_terms is a list of the drift terms to include; currently supported terms
-    # are 'regional_linear', 'point_log', and 'external_Z'. Refer to 
+    # are 'regional_linear', 'point_log', and 'external_Z'. Refer to
     # UniversalKriging.__doc__ for more information.
     UK = UniversalKriging(data[:, 0], data[:, 1], data[:, 2], variogram_model='linear',
                           drift_terms=['regional_linear'])
-                                             
+
     # Creates the kriged grid and the variance grid. Allows for kriging on a rectangular
     # grid of points, on a masked rectangular grid of points, or with arbitrary points.
     # (See UniversalKriging.__doc__ for more information.)
     z, ss = UK.execute('grid', gridx, gridy)
 
+
 Three-Dimensional Kriging Example
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
@@ -140,44 +142,84 @@ Three-Dimensional Kriging Example
     gridy = np.arange(0.0, 0.6, 0.01)
     gridz = np.arange(0.0, 0.6, 0.1)
 
-    # Create the 3D ordinary kriging object and solves for the three-dimension kriged 
+    # Create the 3D ordinary kriging object and solves for the three-dimension kriged
     # volume and variance. Refer to OrdinaryKriging3D.__doc__ for more information.
     ok3d = OrdinaryKriging3D(data[:, 0], data[:, 1], data[:, 2], data[:, 3],
                                                      variogram_model='linear')
     k3d, ss3d = ok3d.execute('grid', gridx, gridy, gridz)
 
-    # Create the 3D universal kriging object and solves for the three-dimension kriged 
+    # Create the 3D universal kriging object and solves for the three-dimension kriged
     # volume and variance. Refer to UniversalKriging3D.__doc__ for more information.
-    uk3d = UniversalKriging3D(data[:, 0], data[:, 1], data[:, 2], data[:, 3], 
+    uk3d = UniversalKriging3D(data[:, 0], data[:, 1], data[:, 2], data[:, 3],
                                                       variogram_model='linear', drift_terms=['regional_linear'])
     k3d, ss3d = uk3d.execute('grid', gridx, gridy, gridz)
 
-    # To use the generic 'specified' drift term, the user must provide the drift values 
-    # at each data point and at every grid point. The following example is equivalent to 
+    # To use the generic 'specified' drift term, the user must provide the drift values
+    # at each data point and at every grid point. The following example is equivalent to
     # using a linear drift in all three spatial dimensions. Refer to
     # UniversalKriging3D.__doc__ for more information.
     zg, yg, xg = np.meshgrid(gridz, gridy, gridx, indexing='ij')
-    uk3d = UniversalKriging3D(data[:, 0], data[:, 1], data[:, 2], data[:, 3], 
+    uk3d = UniversalKriging3D(data[:, 0], data[:, 1], data[:, 2], data[:, 3],
                                                       variogram_model='linear', drift_terms=['specified'],
                                                       specified_drift=[data[:, 0], data[:, 1]])
     k3d, ss3d = uk3d.execute('grid', gridx, gridy, gridz, specified_drift_arrays=[xg, yg, zg])
 
-    # To use the generic 'functional' drift term, the user must provide a callable 
-    # function that takes only the spatial dimensions as arguments. The following example 
-    # is equivalent to using a linear drift only in the x-direction. Refer to 
+    # To use the generic 'functional' drift term, the user must provide a callable
+    # function that takes only the spatial dimensions as arguments. The following example
+    # is equivalent to using a linear drift only in the x-direction. Refer to
     # UniversalKriging3D.__doc__ for more information.
     func = lambda x, y, z: x
-    uk3d = UniversalKriging3D(data[:, 0], data[:, 1], data[:, 2], data[:, 3], 
+    uk3d = UniversalKriging3D(data[:, 0], data[:, 1], data[:, 2], data[:, 3],
                                                       variogram_model='linear', drift_terms=['functional'],
                                                       functional_drift=[func])
     k3d, ss3d = uk3d.execute('grid', gridx, gridy, gridz)
 
-    # Note that the use of the 'specified' and 'functional' generic drift capabilities is 
-    # essentially identical in the two-dimensional universal kriging class (except for a 
-    # difference in the number of spatial coordinates for the passed drift functions). 
+    # Note that the use of the 'specified' and 'functional' generic drift capabilities is
+    # essentially identical in the two-dimensional universal kriging class (except for a
+    # difference in the number of spatial coordinates for the passed drift functions).
     # See UniversalKriging.__doc__ for more information.
 
 
+GSTools covariance models
+^^^^^^^^^^^^^^^^^^^^^^^^^
+
+You can also use `GSTools <https://github.com/GeoStat-Framework/GSTools>`_
+covariance models as input for the ``variogram_model`` in the
+PyKrige routines:
+
+.. code:: python
+
+    import numpy as np
+    from gstools import Gaussian
+    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]])
+    # 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
+    cov_model = Gaussian(dim=2, len_scale=1, anis=0.2, angles=-0.5, var=0.5, nugget=0.1)
+    # ordinary kriging with pykrige
+    OK1 = OrdinaryKriging(data[:, 0], data[:, 1], data[:, 2], cov_model)
+    z1, ss1 = OK1.execute('grid', gridx, gridy)
+    plt.imshow(z1, origin="lower")
+    plt.show()
+
+Which gives:
+
+.. image:: https://raw.githubusercontent.com/GeoStat-Framework/GSTools/master/docs/source/pics/20_pykrige.png
+   :width: 400px
+   :align: center
+
+Have a look at the `documentation about the Covariance Model of GSTools <https://geostat-framework.readthedocs.io/projects/gstools/en/latest/tutorial_02_cov.html>`_.
+
+
 Kriging Parameters Tuning
 ^^^^^^^^^^^^^^^^^^^^^^^^^
 
@@ -194,8 +236,8 @@ Regression Kriging
 with `pykrige.rk.RegressionKriging <http://pykrige.readthedocs.io/en/latest/examples/regression_kriging2d.html>`_.
 This class takes as parameters a scikit-learn regression model, and details of either either
 the ``OrdinaryKriging`` or the ``UniversalKriging`` class, and performs a correction steps on the ML regression prediction.
- 
-A demonstration of the regression kriging is provided in the 
+
+A demonstration of the regression kriging is provided in the
 `corresponding example <http://pykrige.readthedocs.io/en/latest/examples/regression_kriging2d.html#sphx-glr-examples-regression-kriging2d-py>`_.
 
 License

appveyor.yml

@@ -32,9 +32,9 @@ install:
 
 test_script:
   - pytest -sv --cov=pykrige pykrige/
-  - ps: | 
+  - ps: |
           cd ".\\examples"
-          Get-ChildItem ".\\" -Filter *.py | 
+          Get-ChildItem ".\\" -Filter *.py |
           Foreach-Object {
               python $_.FullName 2>&1 | Tee-Object -Append -FilePath "C:\\log.txt"
           }

examples/gstools_covmodel.py

@@ -0,0 +1,41 @@
+# -*- coding: utf-8 -*-
+"""
+GSTools Interface
+=================
+
+Example how to use the PyKrige routines with a GSTools CovModel.
+"""
+import os
+
+import numpy as np
+from pykrige.ok import OrdinaryKriging
+from matplotlib import pyplot as plt
+try:
+    from gstools import Gaussian
+    GS_IMP = True
+except ImportError:
+    GS_IMP = False
+
+# 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]])
+# 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
+if GS_IMP:
+    cov_model = Gaussian(
+        dim=2, len_scale=4, anis=.2, angles=-.5, var=.5, nugget=.1
+    )
+else:
+    cov_model = "gaussian"
+# ordinary kriging with pykrige
+OK1 = OrdinaryKriging(data[:, 0], data[:, 1], data[:, 2], cov_model)
+z1, ss1 = OK1.execute('grid', gridx, gridy)
+plt.imshow(z1, origin="lower")
+if 'CI' not in os.environ:
+    # skip in continous integration
+    plt.show()