kriging¶
Kriging¶
In [1]:
%matplotlib inline
In [2]:
# loading packages
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
geolime package can be imported as such:
In [3]:
import geolime as glm
Data¶
geolime provides random generated data for kriging purposes.
In [4]:
random_points_path = "../../tests/data/random_points.csv"
# load input dataset
data = pd.read_csv(random_points_path)
data.head()
Out[4]:
| Unnamed: 0 | x | y | z | input | |
|---|---|---|---|---|---|
| 0 | 0 | 575.3 | 143.5 | 433.9 | 111 |
| 1 | 1 | 453.6 | 769.6 | 774.6 | 89 |
| 2 | 2 | 866.8 | 245.5 | 353.7 | 17 |
| 3 | 3 | 773.5 | 268.7 | 908.0 | 145 |
| 4 | 4 | 658.7 | 735.3 | 533.4 | 146 |
In [5]:
coords_df = data[['x', 'y', 'z']]
coords = coords_df.to_numpy()
values_df = data[['input']]
values = values_df.to_numpy()
In [6]:
# create 3d grid
grid = glm.create_grid_3d(origin=[-0.75, -0.75, -0.75], cell_size=[1.5, 1.5, 1.5], n_dim=[2, 2, 2], angle=[0, 0, 0])
Here we define a covariance model
In [7]:
cov1 = glm.Spherical(dim=3, angles=[0, 0, 0], scales=[5000, 5000, 5000])
cov2 = glm.Spherical(dim=3, angles=[0, 0, 0], scales=[50, 50, 50])
cov3 = glm.Nugget()
cov = cov3 + 4 * cov1 + 3 * cov2
Here we define a unique neighborhood
In [8]:
neigh = glm.Neighborhood()
We have all the arguments to apply a ponctual kriging (based on the discretization parameter)
In [9]:
result_ok = glm.block_kriging(coords, values, grid.get_corner(), cov, neigh)
result_ok_df = pd.DataFrame(data=result_ok, columns=['estim', 'stdev'])
In [10]:
result_ok_df.head()
Out[10]:
| estim | stdev | |
|---|---|---|
| 0 | 81.428991 | 2.193899 |
| 1 | 81.444129 | 2.193474 |
| 2 | 81.406035 | 2.193527 |
| 3 | 81.421201 | 2.193102 |
| 4 | 81.426982 | 2.193552 |
Block kriging¶
Block kriging works in a similar fashion. One needs to specify the block size using a discretization parameter
In [11]:
pattern_generator = glm.GridPatternGenerator(grid, [5, 5, 5])
result_bk = glm.block_kriging(coords, values, grid.get_center(), cov, neigh, pattern_generator)
result_bk_df = pd.DataFrame(data=result_bk, columns=['estim', 'stdev'])
In [12]:
result_bk_df.head()
Out[12]:
| estim | stdev | |
|---|---|---|
| 0 | 81.424123 | 1.928741 |
| 1 | 81.439338 | 1.928257 |
| 2 | 81.401047 | 1.928320 |
| 3 | 81.416290 | 1.927835 |
| 4 | 81.422113 | 1.928348 |
Simple kriging¶
Simple kriging is operated with the same function. But the known mean is given as an input argument. Below, a 2D example.
In [13]:
# define a covariance model
cov = 3 * glm.Nugget() + 5 * glm.Spherical(2, -45, [30, 10])
# create a 2D grid
target_grid = glm.create_grid_2d(origin=[0, 0], n_dim=[10, 10], cell_size=[10, 10], angle=[0.0, 0.0, 0.0])
# define a unique neighborhood
neigh = glm.Neighborhood()
We use the same random dataset, but without the third dimension
In [14]:
coords_df = data[['x', 'y']]
coords = coords_df.to_numpy()
We can then apply punctual simple kriging
In [15]:
result_sk = glm.block_kriging(coords, values, target_grid.get_corner()[:1000], cov, neigh, mean=10)
result_sk_df = pd.DataFrame(data=result_sk, columns=['estim', 'stdev'])
The resulting dataframe provides the same columns as seen before
In [16]:
result_sk_df.head()
Out[16]:
| estim | stdev | |
|---|---|---|
| 0 | 10.0 | 2.828427 |
| 1 | 10.0 | 2.828427 |
| 2 | 10.0 | 2.828427 |
| 3 | 10.0 | 2.828427 |
| 4 | 10.0 | 2.828427 |
It is possible to display the kriging estimation (only in 2D) using the following call:
In [17]:
glm.display_grid(target_grid, result_sk_df['estim'].to_numpy(), coords, z=None, title="Simple kriging")
Out[17]:
(<Figure size 432x288 with 2 Axes>,
<AxesSubplot:title={'center':'Simple kriging'}>)
Discrete Gaussian Kriging¶
This kriging method uses alternative API functions. We will retain the same dataset for demonstration purposes.
In [18]:
# define a covariance model
cov = glm.Spherical(2, 0, [30, 30])
In [19]:
# create a 2D grid
target_grid = glm.create_grid_2d(origin=[0, 0], n_dim=[10, 10], cell_size=[10, 10], angle=[0.0, 0.0, 0.0])
In [20]:
# define a unique neighborhood
neigh = glm.Neighborhood()
In [21]:
# compute change of support coefficient
pattern_generator = glm.GridPatternGenerator(target_grid, [5, 5])
r = np.sqrt(glm.cvv(cov, pattern_generator))
In [22]:
# retrieve center coordinates, number of points in block and mean by block
xcoordsc, n_points_in_blocks, mean_by_block = glm.group_data(target_grid, coords, values)
WARNING: some points are not included in the union of the blocks. np.nan is returned as index in such case.
In [23]:
# apply Discrete Gaussian Kriging
result_dgk = glm.discrete_gaussian_kriging_by_index(target_grid.get_center(), xcoordsc, mean_by_block, n_points_in_blocks, cov, neigh, r)
result_dgk_df = pd.DataFrame({'estim': result_dgk[:, 0], 'stdev': result_dgk[:, 1]})
In [24]:
result_dgk_df.head()
Out[24]:
| estim | stdev | |
|---|---|---|
| 0 | 0.000000 | 1.000000 |
| 1 | 0.000000 | 1.000000 |
| 2 | 0.000000 | 1.000000 |
| 3 | 0.000000 | 1.000000 |
| 4 | 8.153204 | 0.997037 |