I am reading the book from Isaaks and Srivastava "An Introduction to Applied Geostatistics" and I am coding a small library in Python to perform kriging. I know there are other tools but I wanted to do it my way in order to fully understand the method.

In the part 12, they explain how the kriging sytem works and how to solve it for different scenarios with a toy example of 7 samples and 1 point to evaluate. This is perfect for my library for it enables to have a test which asserts that my code is working correctly.

I have therefore been reproducing their example when you change the shape / range / sill / nugget of a uniform variogram and would like to do it also for the anisotropic scenario, but I can't succeed in reproducing their results.

I have linked the step by step example I have written to match their experiment but even when changing the theta angle, I can't have the corresponding kriging weights given in the book.

I was wondering if someone has a small toy example of directional variogram to be sure my library was doing the proper calculations.

import numpy as np

# values in the book
points_lake = np.asarray([[61, 139],
                          [63, 140],
                          [64, 129],
                          [68, 128],
                          [71, 140],
                          [73, 141],
                          [75, 128]])
values_lake = np.asarray([477, 696, 227, 646, 606, 791, 783])
point_unknown = np.asarray([[65, 137]])

# Exponential variogram model
f = lambda x: 1-np.exp(-3*x)

# Anisotropy axis
theta = np.pi/4

# Rotation and translation matrices. The anisotropy is modeled as a range of 10 for the principal axis and 5 for the perpendicular
R = np.asarray([[np.cos(theta), np.sin(theta)], [-np.sin(theta), np.cos(theta)]])
T = np.asarray([[.1, 0], [0, .2]])

# Semivariance matrix
vecs = np.asarray([p-q for p in points_lake for q in points_lake])
tr_vecs = T.dot(R.dot(vecs.T)).T
lags = np.linalg.norm(tr_vecs, axis=1)
S = np.asarray([f(x) for x in lags]).reshape(7, 7)
S = np.concatenate([np.concatenate([S, np.ones((1, S.shape[0]))]),
                                    np.concatenate([np.ones((S.shape[0], 1)),
                                        np.atleast_2d(0)])], axis=1)
Sinv = np.linalg.inv(S)

# Semivariance for point to estimate
vec_k = np.asarray([q - point_unknown[0] for q in points_lake])
tr_vec_k = T.dot(R.dot(vec_k.T)).T
lags_k = np.linalg.norm(tr_vec_k, axis=1)
D = np.asarray([f(x) for x in lags_k])
D = np.concatenate([D, [1]])

weights = np.matmul(Sinv, D)
# should be either [.24, .44, .05, .09, .07, .02, .09] or [.11, .29, .16, .05, .29, .04, .06]

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.