Skip to content

Commit

Permalink
add missing tutorials
Browse files Browse the repository at this point in the history
  • Loading branch information
MuellerSeb committed Aug 28, 2019
1 parent 2599541 commit 0822a89
Show file tree
Hide file tree
Showing 14 changed files with 432 additions and 2 deletions.
12 changes: 10 additions & 2 deletions docs/source/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -30,8 +30,16 @@ To get the latest development version you can install it directly from GitHub:
pip install https://github.com/GeoStat-Framework/GSTools/archive/master.zip
To enable the OpenMP support, you have to provide a C compiler, Cython and OpenMP.
To get all other dependencies, it is recommended to first install gstools once
in the standard way just decribed.
Then use the following command:

.. code-block:: none
pip install --global-option="--openmp" gstools
Or for the development version:

.. code-block:: none
pip install --global-option="--openmp" https://github.com/GeoStat-Framework/GSTools/archive/master.zip
Expand Down Expand Up @@ -82,7 +90,7 @@ with a :any:`Gaussian` covariance model.
:align: center

A similar example but for a three dimensional field is exported to a
`VTK <https://vtk.org/>`_ file, which can be visualized with
`VTK <https://vtk.org/>`__ file, which can be visualized with
`ParaView <https://www.paraview.org/>`_.

.. code-block:: python
Expand Down Expand Up @@ -275,7 +283,7 @@ VTK Export
==========

After you have created a field, you may want to save it to file, so we provide
a handy `VTK <https://www.vtk.org/>`_ export routine (:any:`vtk_export`):
a handy `VTK <https://www.vtk.org/>`__ export routine:

.. code-block:: python
Expand Down
Binary file added docs/source/pics/05_ordinary.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added docs/source/pics/05_simple.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added docs/source/pics/06_ensemble.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added docs/source/pics/07_00_std.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added docs/source/pics/07_01_lognormal.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added docs/source/pics/07_02_binary.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added docs/source/pics/07_03_zinnharvey.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added docs/source/pics/07_04_arcsin.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added docs/source/pics/07_05_combine.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
153 changes: 153 additions & 0 deletions docs/source/tutorial_05_kriging.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,153 @@
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.


Theoretical Background
----------------------

Aim of kriging is to derive the value of a field at some point :math:`x_0`,
when there are fix observed values :math:`z(x_1)\ldots z(x_n)` at given points :math:`x_i`.

The resluting value :math:`z_0` at :math:`x_0` is calculated as a weighted mean:

.. math::
z_0 = \sum_{i=1}^n w_i(x_0) \cdot z_i
The weights :math:`W = (w_1,\ldots,w_n)` depent on the given covariance model and the location of the target point.

The different kriging approaches provide different ways of calculating :math:`W`.


Implementation
--------------

The routines for kriging are almost identical to the routines for spatial random fields.
First you define a covariance model, as described in the SRF Tutorial,
then you initialize the kriging class with this model:

.. code-block:: python
from gstools import Gaussian, krige
# 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)
The resulting field instance ``krig`` has the same methods as the :any:`SRF` class.
You can call it to evaluate the kriging 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`.


Simple Kriging
--------------

Simple kriging assumes a known mean of the data.
For simplicity we assume a mean of 0,
which can be achieved by subtracting the mean from the observed values and
subsequently adding it to the resulting data.

The resulting equation system for :math:`W` is given by:

.. math::
W = \begin{pmatrix}c(x_1,x_1) & \cdots & c(x_1,x_n) \\
\vdots & \ddots & \vdots \\
c(x_n,x_1) & \cdots & c(x_n,x_n)
\end{pmatrix}^{-1}
\begin{pmatrix}c(x_1,x_0) \\ \vdots \\ c(x_n,x_0) \end{pmatrix}
Thereby :math:`c(x_i,x_j)` is the covariance of the given observations.


Example
^^^^^^^

Here we use simple kriging in 1D (for plotting reasons) with 5 given observations/condtions.
The mean of the field has to be given beforehand.

.. code-block:: python
import numpy as np
from gstools import Gaussian, krige
# condtions
cond_pos = [0.3, 1.9, 1.1, 3.3, 4.7]
cond_val = [0.47, 0.56, 0.74, 1.47, 1.74]
# resulting grid
gridx = np.linspace(0.0, 15.0, 151)
# spatial random field class
model = Gaussian(dim=1, var=0.5, len_scale=2)
krig = krige.Simple(model, mean=1, 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()
.. image:: pics/05_simple.png
:width: 600px
:align: center


Ordinary Kriging
----------------

Ordinary kriging will estimate an appropriate mean of the field,
based on the given observations/conditions and the covariance model in use.

The resulting equation system for :math:`W` is given by:

.. math::
\begin{pmatrix}W\\\mu\end{pmatrix} = \begin{pmatrix}
\gamma(x_1,x_1) & \cdots & \gamma(x_1,x_n) &1 \\
\vdots & \ddots & \vdots & \vdots \\
\gamma(x_n,x_1) & \cdots & \gamma(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}
Thereby :math:`\gamma(x_i,x_j)` is the semi-variogram of the given observations
and :math:`\mu` is a Lagrange multiplier to minimize the kriging error and estimate the mean.


Example
^^^^^^^

Here we use ordinary kriging in 1D (for plotting reasons) with 5 given observations/condtions.
The estimated mean can be accessed by ``krig.mean``.

.. code-block:: python
import numpy as np
from gstools import Gaussian, krige
# condtions
cond_pos = [0.3, 1.9, 1.1, 3.3, 4.7]
cond_val = [0.47, 0.56, 0.74, 1.47, 1.74]
# resulting grid
gridx = np.linspace(0.0, 15.0, 151)
# spatial random field class
model = Gaussian(dim=1, var=0.5, len_scale=2)
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()
.. image:: pics/05_ordinary.png
:width: 600px
:align: center


.. raw:: latex

\clearpage
79 changes: 79 additions & 0 deletions docs/source/tutorial_06_conditioning.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,79 @@
Tutorial 6: Conditioned Fields
==============================

Kriging fields tend to approach the field mean outside the area of observations.
To generate random fields, that coincide with given observations, but are still
random in the area outside, we provide the generation of conditioned random fields.


Theoretical Background
----------------------

The idea behind conditioned random fields builts up on kriging.
First we generate a field with kriging, then we generate a random field
and finally we generate another kriging field to eliminate the error between
the random field and teh kriging field of the given observations.

To do so, you can choose between ordinary and simple kriging.
In case of ordinary kriging, the mean of the SRF will be overwritten by the
estimated mean.

The setup of the spatial random field is the same like described in the
srf tutorial.
You just need to add the conditions like described in the kriging tutorial:

.. code-block:: python
srf.set_condition([cond_pos], cond_val, "simple")
or:

.. code-block:: python
srf.set_condition([cond_pos], cond_val, "ordinary")
Example: Conditioning with Ordinary Kriging
-------------------------------------------

Here we use ordinary kriging in 1D (for plotting reasons) with 5 given observations/condtions,
to generate an ensemble of conditioned random fields.
The estimated mean can be accessed by ``srf.mean``.

.. code-block:: python
import numpy as np
from gstools import Gaussian, SRF
import matplotlib.pyplot as plt
# condtions
cond_pos = [0.3, 1.9, 1.1, 3.3, 4.7]
cond_val = [0.47, 0.56, 0.74, 1.47, 1.74]
gridx = np.linspace(0.0, 15.0, 151)
# spatial random field class
model = Gaussian(dim=1, var=0.5, len_scale=2)
srf = SRF(model)
srf.set_condition([cond_pos], cond_val, "ordinary")
fields = []
for i in range(100):
if i % 10 == 0: print(i)
fields.append(srf([gridx], seed=i))
label = "Conditioned ensemble" if i == 0 else None
plt.plot(gridx, fields[i], color="k", alpha=0.1, label=label)
plt.plot(gridx, np.full_like(gridx, srf.mean), label="estimated mean")
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")
plt.legend()
plt.show()
.. image:: pics/06_ensemble.png
:width: 600px
:align: center

As you can see, the kriging field coincides with the ensemble mean of the
conditioned random fields and the estimated mean is the mean of the far-field.


.. raw:: latex

\clearpage
Loading

0 comments on commit 0822a89

Please sign in to comment.