B.2.1 Directional Ordinary Kriging

Table of Contents:

1. Read meuse dataset, transform data,

2. Create directional semivariograms,

3. Create directional kriging models,

4. Compare results of interpolation.

Introduction

In this tutorial, we will learn about Ordinary Kriging with directions. The classic Kriging works with an isotropic semivariogram. The similarity between points has the same weights in all directions. Actual geographical data is rarely isotropic, and we distinguish leading directions in spatial processes. The example may be a temperature, where differences in the North-South axis are greater than in the West-East axis.

1. Read meuse dataset and transform data

We use meuse dataset and interpolate zinc concentrations. Dataset is provided in:

Pebesma, Edzer. (2009). The meuse data set: a tutorial for the gstat R package -> link to the publication

import geopandas as gpd
import numpy as np
import pandas as pd
import pyinterpolate as ptp
from pyinterpolate.variogram.empirical.experimental_variogram import DirectionalVariogram
import matplotlib.pyplot as plt

MEUSE_FILE = 'samples/point_data/csv/meuse/meuse.csv'
MEUSE_GRID = 'samples/point_data/csv/meuse/meuse_grid.csv'

# Variogram parameters
STEP_SIZE = 100
MAX_RANGE = 1600

df = pd.read_csv(MEUSE_FILE, usecols=['x', 'y', 'zinc'])

df.head()

	x	y	zinc
0	181072	333611	1022
1	181025	333558	1141
2	181165	333537	640
3	181298	333484	257
4	181307	333330	269

Let’s check the distribution of zinc values:

df['zinc'].plot(kind='hist')

<Axes: ylabel='Frequency'>

Our zinc concentrations are highly skewed. We should take a natural logarithm of values to make the distribution closer to the normal.

df['log-zinc'] = np.log(df['zinc'])

df['log-zinc'].plot(kind='hist')

<Axes: ylabel='Frequency'>

The histogram is still not close to the normal distribution; it has two modes, but the tail is shorter. We will build models based on this histogram.

2. Create directional variograms

Data is prepared, we can model it’s semivariance. We will perform an isotropic modeling and a full directional modeling in N-S, W-E, NE-SW, and NW-SE axes.

dirvar = DirectionalVariogram(input_array=df[['x', 'y', 'log-zinc']],
                              step_size=STEP_SIZE,
                              max_range=MAX_RANGE,
                              tolerance=0.1,
                              method='e')
dirvar.show()

We have five experimental variograms for further modeling. Let’s build theoretical models from those.

MODEL_NAME = 'spherical'
DIRECTIONS = ['ISO', 'NS', 'WE', 'NE-SW', 'NW-SE']

theoretical_models = {}
variograms = dirvar.get()

for direction in DIRECTIONS:
    theoretical_model = ptp.TheoreticalVariogram()
    theoretical_model.autofit(experimental_variogram=variograms[direction],
                              model_name=MODEL_NAME,
                              nugget=0)
    theoretical_models[direction] = theoretical_model

3. Create directional Ordinary Kriging models

With a set of directional variograms we will perform spatial interpolation over a defined grid.

grid = pd.read_csv(MEUSE_GRID, usecols=['x', 'y'])
grid.head()

	x	y
0	181180	333740
1	181140	333700
2	181180	333700
3	181220	333700
4	181100	333660

kriged_results = {}

for direction in DIRECTIONS:
    result = ptp.kriging(observations=df[['x', 'y', 'log-zinc']].values,
                         theoretical_model=theoretical_models[direction],
                         points=grid.values,
                         no_neighbors=16,
                         use_all_neighbors_in_range=True)
    result = pd.DataFrame(result, columns=['pred', 'err', 'x', 'y'])
    rgeometry = gpd.points_from_xy(result['x'], result['y'])
    rgeometry.name = 'geometry'
    result = gpd.GeoDataFrame(result, geometry=rgeometry)

    kriged_results[direction] = result

100%|██████████| 3103/3103 [00:03<00:00, 776.43it/s]
100%|██████████| 3103/3103 [00:00<00:00, 4160.84it/s]
100%|██████████| 3103/3103 [00:00<00:00, 3911.17it/s]
100%|██████████| 3103/3103 [00:00<00:00, 3740.61it/s]
100%|██████████| 3103/3103 [00:00<00:00, 4165.74it/s]

kriged_results['ISO'].head()

	pred	err	x	y	geometry
0	6.523941	0.283231	181180.0	333740.0	POINT (181180.000 333740.000)
1	6.664006	0.198841	181140.0	333700.0	POINT (181140.000 333700.000)
2	6.541937	0.224790	181180.0	333700.0	POINT (181180.000 333700.000)
3	6.407741	0.254611	181220.0	333700.0	POINT (181220.000 333700.000)
4	6.809450	0.109191	181100.0	333660.0	POINT (181100.000 333660.000)

With all models calculated, we can compare the results of interpolation. We will show interpolation and variance errors across the models.

4. Compare interpolated results

fig, ax = plt.subplots(nrows=5, ncols=2, figsize=(15, 20))

kriged_results['ISO'].plot(ax=ax[0, 0], column='pred', cmap='cividis')
kriged_results['ISO'].plot(ax=ax[0, 1], column='err', cmap='Reds')
ax[0, 0].set_title('Prediction - isotropic')
ax[0, 1].set_title('Error - isotropic')

kriged_results['NS'].plot(ax=ax[1, 0], column='pred', cmap='cividis')
kriged_results['NS'].plot(ax=ax[1, 1], column='err', cmap='Reds')
ax[1, 0].set_title('Prediction - NS')
ax[1, 1].set_title('Error - NS')

kriged_results['WE'].plot(ax=ax[2, 0], column='pred', cmap='cividis')
kriged_results['WE'].plot(ax=ax[2, 1], column='err', cmap='Reds')
ax[2, 0].set_title('Prediction - WE')
ax[2, 1].set_title('Error - WE')

kriged_results['NE-SW'].plot(ax=ax[3, 0], column='pred', cmap='cividis')
kriged_results['NE-SW'].plot(ax=ax[3, 1], column='err', cmap='Reds')
ax[3, 0].set_title('Prediction - NE-SW')
ax[3, 1].set_title('Error - NE-SW')

kriged_results['NW-SE'].plot(ax=ax[4, 0], column='pred', cmap='cividis')
kriged_results['NW-SE'].plot(ax=ax[4, 1], column='err', cmap='Reds')
ax[4, 0].set_title('Prediction - NW-SE')
ax[4, 1].set_title('Error - NW-SE')

plt.show()

Directional Kriging clearly shows the leading direction and is even more pronounced in the variance error maps. At the end, let’s compare the mean of directional interpolations to the isotropic interpolation result:

mean_directional_results = (kriged_results['NS']['pred'] +
                            kriged_results['WE']['pred'] +
                            kriged_results['NE-SW']['pred'] +
                            kriged_results['NW-SE']['pred']) / 4

kriged_results['ISO']['mean_dir_pred'] = mean_directional_results

_, ax = plt.subplots(nrows=1, ncols=2, figsize=(15, 20))

kriged_results['ISO'].plot(ax=ax[0], column='pred', cmap='cividis')
ax[0].set_title('Prediction - isotropic')
kriged_results['ISO'].plot(ax=ax[1], column='mean_dir_pred', cmap='cividis')
ax[1].set_title('Prediction - mean directional results')
plt.show()

---

Where to go from here?

• C.1.1 Blocks to Points Interpolation with Ordinary Kriging

• C.1.2 Semivariogram Regularization

Changelog

Date	Change description	Author	Package Version
2023-08-23	The tutorial was refreshed and set along with the 0.5.0 version of the package	@SimonMolinsky	0.5.0
2023-04-15	Update: version 0.4.1	@SimonMolinsky	0.4.1
2022-11-19	The first version of tutorial	@SimonMolinsky	0.3.6