-
Notifications
You must be signed in to change notification settings - Fork 71
/
01_ordinary_kriging.py
48 lines (37 loc) · 1.51 KB
/
01_ordinary_kriging.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
r"""
Ordinary Kriging
----------------
Ordinary kriging will estimate an appropriate mean of the field,
based on the given observations/conditions and the covariance model used.
The resulting system of equations 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/conditions.
The estimated mean can be accessed by ``krig.mean``.
"""
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()