Menu
▾
▴
Kriging
Attachments
image.png
(39906 bytes)
image2.png
(39906 bytes)
Kriging
In statistics, originally in geostatistics, kriging or Gaussian process regression is a method of interpolation for which the interpolated values are modeled by a Gaussian process governed by prior covariances, as opposed to a piecewise-polynomial spline chosen to optimize smoothness of the fitted values. Under suitable assumptions on the priors, kriging gives the best linear unbiased prediction of the intermediate values. Interpolating methods based on other criteria such as smoothness need not yield the most likely intermediate values. The method is widely used in the domain of spatial analysis and computer experiments. The technique is also known as Kolmogorov Wiener prediction.
in wikipedia - en
see also wikipedia - pt
Algorithm of simple kriging in Python
For demonstrations purposes the following algorithm does not take into account a search engine for points (all are used for all node) and the variogram is very simplified (no rotations, you only give ranges in X and Y directions and also the model used is always exponential).
That said simple kriging (SK) can be done in the following manor. We calculate the angles and distance between all points. From the angle we take the ranges required for those directions, and with that we calculate the specific variogram value for that pair of points. When finished we should have a table with all variogram values between all points. From now on we do:
1) Check the distance between node and samples.
2) Check the angle between node and sample.
3) Calculate the M array for variogram values between node and sample.
4) Get the K matrix with the variogram values of all points envolved.
5) Solve the system K*w=M. This will give you the weights.
6) Multiply the weights by the result beetween values minus the mean of values (it could be an user input). Sum the result.
7) Subtract the mean of values to that sum and you have the value for the node.
8) Repeat the steps for all nodes.
Notice the following examples were adapted from here.
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
import scipy.linalg as LA
def SK(x,y,v,variogram,grid):
cov_angulos = np.zeros((x.shape[0],x.shape[0]))
cov_distancias = np.zeros((x.shape[0],x.shape[0]))
K = np.zeros((x.shape[0],x.shape[0]))
for i in xrange(x.shape[0]-1):
cov_angulos[i,i:]=np.arctan2((y[i:]-y[i]),(x[i:]-x[i]))
cov_distancias[i,i:]=np.sqrt((x[i:]-x[i])**2+(y[i:]-y[i])**2)
for i in xrange(x.shape[0]):
for j in xrange(x.shape[0]):
if cov_distancias[i,j]!=0:
amp=np.sqrt((variogram[1]*np.cos(cov_angulos[i,j]))**2+(variogram[0]*np.sin(cov_angulos[i,j]))**2)
K[i,j]=v[:].var()*(1-np.e**(-3*cov_distancias[i,j]/amp))
K = K + K.T
for i in xrange(grid.shape[0]):
for j in xrange(grid.shape[1]):
distancias = np.sqrt((i-x[:])**2+(j-y[:])**2)
angulos = np.arctan2(i-y[:],j-x[:])
amplitudes = np.sqrt((variogram[1]*np.cos(angulos[:]))**2+(variogram[0]*np.sin(angulos[:]))**2)
M = v[:].var()*(1-np.e**(-3*distancias[:]/amplitudes[:]))
W = LA.solve(K,M)
grid[i,j] = np.sum(W*(v[:]-v[:].mean()))+v[:].mean()
return grid
np.random.seed(123433789) # GIVING A SEED NUMBER FOR THE EXPERIENCE TO BE REPRODUCIBLE
grid = np.zeros((100,100),dtype='float32') # float32 gives us a lot precision
x,y = np.random.randint(0,100,10),np.random.randint(0,100,10) # CREATE POINT SET.
v = np.random.randint(0,10,10) # THIS IS MY VARIABLE
grid = SK(x,y,v,(50,30),grid)
plt.imshow(grid.T,origin='lower',interpolation='nearest',cmap='jet')
plt.scatter(x,y,c=v,cmap='jet',s=120)
plt.xlim(0,grid.shape[0])
plt.ylim(0,grid.shape[1])
plt.grid()
plt.show()
And the result is:
For ordinary kriging (OK) we do pretty much the same but add a new row to the K matrix (and a new element to M) as can be seen in the following example (also notice the result is the direct multiplication of weigths with values):
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
import scipy.linalg as LA
def OK(x,y,v,variogram,grid):
cov_angulos = np.zeros((x.shape[0],x.shape[0]))
cov_distancias = np.zeros((x.shape[0],x.shape[0]))
K = np.zeros((x.shape[0]+1,x.shape[0]+1))
for i in xrange(x.shape[0]-1):
cov_angulos[i,i:]=np.arctan2((y[i:]-y[i]),(x[i:]-x[i]))
cov_distancias[i,i:]=np.sqrt((x[i:]-x[i])**2+(y[i:]-y[i])**2)
for i in xrange(x.shape[0]):
for j in xrange(x.shape[0]):
if cov_distancias[i,j]!=0:
amp=np.sqrt((variogram[1]*np.cos(cov_angulos[i,j]))**2+(variogram[0]*np.sin(cov_angulos[i,j]))**2)
K[i,j]=v[:].var()*(1-np.e**(-3*cov_distancias[i,j]/amp))
K = K + K.T
K[-1,:] = 1
K[:,-1] = 1
K[-1,-1] = 0
for i in xrange(grid.shape[0]):
for j in xrange(grid.shape[1]):
distancias = np.sqrt((i-x[:])**2+(j-y[:])**2)
angulos = np.arctan2(i-y[:],j-x[:])
amplitudes = np.sqrt((variogram[1]*np.cos(angulos[:]))**2+(variogram[0]*np.sin(angulos[:]))**2)
M = np.ones((x.shape[0]+1,1))
M[:-1,0] = v[:].var()*(1-np.e**(-3*distancias[:]/amplitudes[:]))
W = LA.solve(K,M)
grid[i,j] = np.sum(W[:-1,0]*(v[:]))
return grid
grid = OK(x,y,v,(50,30),grid)
plt.imshow(grid.T,origin='lower',interpolation='nearest',cmap='jet')
plt.scatter(x,y,c=v,cmap='jet',s=120)
plt.xlim(0,grid.shape[0])
plt.ylim(0,grid.shape[1])
plt.grid()
plt.show()
And the result is:
Important considerations
The variogram table is normally given by three rotations: azimuth, dip and rake. The examples above do not take that into consideration.
Also there are quite a few methods to search the nearest samples. The examples above just use all available.
Finally the variogram matrices should be correlogram matrices. Once again for simplicity we've used variograms (can cause system instability).
See also
Related
Wiki: Code snippets
Wiki: Inverse weighted distance
Wiki: Nearest neighbor