Kriging¶
Purpose¶
This notebook describes step-by-step guidelines used for applying kriging algorithm using Deeplime's GeoLime geostatistics library. A widespread spatial dataset is used as en entrypoint for applying the GeoLime API for kriging. This datased will first be statistically described prior to apply ordinary and simple kriging (punctual or by block). The outputs will then be displayed numerically and graphically for interpretation.
This notebook mostly details how to achieve kriging outputs and visualization with GeoLime, and doesn't not emphasis on the theoretical fundamentals of kriging formulas and algorithm, although basic knowledge of these estimation techniques are given as insights throughout this notebook.
For starters, on needs to prepare this Python notebook and import GeoLime.
The GeoLime library can be imported as such
import geolime as geo
import plotly.io as pio
pio.renderers.default = "sphinx_gallery"
Dataset¶
The Walker Lake dataset openly available for educational purposes will be used in this notebook. Several geostatistics papers and acedemics use this widespread dataset for applying and analyzing geostatistics concepts, algorithm and models. This dataset is available as a pandas DataFrame within GeoLime and can be loaded with Python using the following commands.
data = geo.datasets.load("walker_lake/walker_lake_sample.csv")
We notice this dataset contains two generic variables of interest (u and v as concentration in ppm, a chemical or contaminate mass measurement per unit mass) and one indicator variable (t).
data.head()
| x | y | z | v | u | t | |
|---|---|---|---|---|---|---|
| 0 | 11.0 | 8.0 | 0.0 | 0.0 | NaN | 2 |
| 1 | 8.0 | 30.0 | 0.0 | 0.0 | NaN | 2 |
| 2 | 9.0 | 48.0 | 0.0 | 224.4 | NaN | 2 |
| 3 | 8.0 | 69.0 | 0.0 | 434.4 | NaN | 2 |
| 4 | 9.0 | 90.0 | 0.0 | 412.1 | NaN | 2 |
The following command displays basic data statistics, that will be used later on specifically for simple kriging.
data.describe()
| x | y | z | v | u | t | |
|---|---|---|---|---|---|---|
| count | 470.000000 | 470.000000 | 470.0 | 470.000000 | 275.000000 | 470.000000 |
| mean | 111.089362 | 141.221277 | 0.0 | 435.298723 | 604.081091 | 1.904255 |
| std | 70.894682 | 77.968954 | 0.0 | 299.882302 | 767.405620 | 0.294554 |
| min | 8.000000 | 8.000000 | 0.0 | 0.000000 | 0.000000 | 1.000000 |
| 25% | 51.000000 | 80.000000 | 0.0 | 184.600000 | 82.150000 | 2.000000 |
| 50% | 89.000000 | 139.500000 | 0.0 | 424.000000 | 319.300000 | 2.000000 |
| 75% | 170.000000 | 208.000000 | 0.0 | 640.850000 | 844.550000 | 2.000000 |
| max | 251.000000 | 291.000000 | 0.0 | 1528.100000 | 5190.100000 | 2.000000 |
Coordinates and variables of interest must first be converted to GeoLime objects before applying kriging API functions.
point_cloud = geo.PointCloud(
name='WalkerLake',
xyz=data[['x', 'y', 'z']]
)
point_cloud.set_property(
name='V',
data=data['v']
)
Definition of the interpolation grid.
grid = geo.Voxel(
name='WalkerLake',
shape=[26, 30],
axis=[[10., 0., 0.], [0., 10., 0.]],
)
Simple kriging¶
Simple kriging is the starting point for spatial data estimation. Kriging is the best linear estimator for interpolation. But the mean used for kriging weights computation, which are used to estimate the unknown points on a target frid, is considered to be known prior to applying the kriging algorithm. It is then given as an input argument for the main kriging API function. This function takes several input arguments, and more particularly a target grid that defines the estimation resolution. It also requires the covariance model computed beneath and a neighborhood. This function outputs estimation values for the input grid and standard deviation for each grid points. Below, a 2D example on the u variable of interest.
Parameters¶
covariance = (
7000 * geo.Nugget(dimension=2)
+ 88000 * geo.Spherical(
dimension=2,
metric=geo.Metric(angles=160,scales=[45,45])
)
)
neighborhood = geo.MinMaxPointsNeighborhood(
dimension=2,
min_n=1,
max_n=25
)
Initialize Solver¶
skriging_estimator = geo.SimpleKriging(
covariance_model=covariance,
neighborhood_model=neighborhood,
axes=[geo.Coord.U, geo.Coord.V],
mean=point_cloud.data('V').mean()
)
Solving¶
skriging_estimator.solve(
obj=point_cloud,
obj_region=None,
obj_attribute='V',
support=grid,
support_region=None,
support_attribute='sk'
)
Plotting Results¶
Display of Kriging results.
geo.plot(grid, property="sk", agg_method="mean", width=600, height=600)
Display of Kriging variance.
geo.plot(grid, property="sk_variance", agg_method="mean", width=600, height=600)
Ordinary kriging¶
Ordinary kriging provides a more general approach compared to simple kriging because no assumptions are provided for the variable of interest mean. More specifically, kriging is also applied to estimate the mean of the variable of interest to be estimated. Ordinary kriging applies to kriging steps, one on the variable's mean and another on the spatial grid using the former result in the same fashion as seen in the simple kriging section. This section shows how to apply punctual and block kriging.
Punctual kriging¶
Using the previously defined target grid, and still considering the same unique neighborhood as input argument and covariance model fitted in the variography section, ordinary kriging is applied on the entire target grid with the same function as seen in the simple kriging section, minus the optional argument. If mean is not supplied during the API call, ordinry kriging will be applied as the default kriging method.
The pattern argument previously defined allows the user to define the block size. This list has been set to be defaulted for punctual kriging, meaning that we don't subdivide the grid for the moment. Ordinary kriging outputs follows the same coding practice seen before, and the output stays similar to simple kriging output.
Initialize Solver¶
okriging_estimator = geo.OrdinaryKriging(
covariance_model=covariance,
neighborhood_model=neighborhood,
axes = [geo.Coord.U, geo.Coord.V]
)
Solving¶
okriging_estimator.solve(
obj=point_cloud,
obj_region=None,
obj_attribute='V',
support=grid,
support_region=None,
support_attribute='ok'
)
Plotting Results¶
geo.plot(grid, property="ok", agg_method="mean", width=600, height=600)
geo.plot(grid, property="ok_variance", agg_method="mean", width=600, height=600)
Block kriging¶
Block kriging uses the same estimator, but is now defined through a discretization parameter. If we consider 10 meters square blocks, block kriging can effectively be computed by modyfing the input discretization. The discretization will define how many points are estimated in a block (here 5x5). Considering 10x10 blocks, the discretization will average 25 estimates in each block.
Here the Block Kriging is performed on the previously defined Ordinary Kriging estimator.
Solving¶
okriging_estimator.solve(
obj=point_cloud,
obj_region=None,
obj_attribute='V',
support=grid,
support_region=None,
support_attribute='bok',
discr=[5,5]
)
geo.plot(grid, property="bok", agg_method="mean", width=600, height=600)