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))]), np.concatenate([np.ones((S.shape, 1)), np.atleast_2d(0)])], axis=1) Sinv = np.linalg.inv(S) # Semivariance for point to estimate vec_k = np.asarray([q - point_unknown 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, ]) weights = np.matmul(Sinv, D) # should be either [.24, .44, .05, .09, .07, .02, .09] or [.11, .29, .16, .05, .29, .04, .06] print(weights[:-1])