# Simple Kriging in Python

In this post I will work through an example of Simple Kriging. Kriging is a set of techniques for interpolation. It differs from other interpolation techniques in that it sacrifices smoothness for the integrity of sampled points. Most interpolation techniques will over or undershoot the value of the function at sampled locations, but kriging honors those measurements and keeps them fixed. In future posts I would like to cover other types of kriging, other semivariaogram models, and colocated co-kriging. Until then, I’m keeping relatively up to date code at my GitHub project, geostatsmodels.

# Data

The data used in this exercise is in a zip file at this site. (Click on Geostatistics Resources.) I have adapted my material from the Kriging document on the same site. My approach will focus more on programming. First we will import some modules and then load the data and parse it,

from pylab import *
import numpy as np
from pandas import DataFrame, Series
from scipy.spatial.distance import pdist, squareform

z = [ i.strip().split() for i in z[10:] ]
z = np.array( z, dtype=np.float )
z = DataFrame( z, columns=['x','y','thk','por','perm','lperm','lpermp','lpermr'] )


Next, we will plot the data,

fig, ax = subplots()
ax.scatter( z.x, z.y, c=z.por, cmap='gray' )
ax.set_aspect(1)
xlim(-1500,22000)
ylim(-1500,17500)
xlabel('Easting [m]')
ylabel('Northing [m]')
title('Porosity %') ; In surveys, we generally specify one point in latitude and longitude, and then measure things as North and East of that point, hence the Northing and Easting.

# The Semivariogram

The semivariogram encodes data about spatial variance over the region at a given distance or lag. We generally expect data points that are close together spatially to share other characteristics, and we expect points that are separated by greater distances to have lesser correlation. The semivariogram allows us to model the similarity points in a field as a function of distance. The semivariogram is given by,

(1) Here, is distance specified by the user, and and are two points that are separated spatially by . The term is the number of points we have that are separated by the distance . The semivariogram then is the sum of squared differences between values separated by a distance . As an aside, contrast this with the formulation for variance,

(2) Here, is the number of data points, is the sample mean, and is a data point. For sample variance, we are taking the squared difference between data points and the mean, and in the semivariogram we are taking the squared difference between data points separated by distance . We can write some functions to calculate the semivariogram at one lag, and then at multiple lags as follows.

def SVh( P, h, bw ):
'''
Experimental semivariogram for a single lag
'''
pd = squareform( pdist( P[:,:2] ) )
N = pd.shape
Z = list()
for i in range(N):
for j in range(i+1,N):
if( pd[i,j] >= h-bw )and( pd[i,j] <= h+bw ):
Z.append( ( P[i,2] - P[j,2] )**2.0 )
return np.sum( Z ) / ( 2.0 * len( Z ) )

def SV( P, hs, bw ):
'''
Experimental variogram for a collection of lags
'''
sv = list()
for h in hs:
sv.append( SVh( P, h, bw ) )
sv = [ [ hs[i], sv[i] ] for i in range( len( hs ) ) if sv[i] > 0 ]
return np.array( sv ).T

def C( P, h, bw ):
'''
Calculate the sill
'''
c0 = np.var( P[:,2] )
if h == 0:
return c0
return c0 - SVh( P, h, bw )


The C() function is the covariance function, and will be used later. Let us now calculate and plot the semivariogram,

# part of our data set recording porosity
P = np.array( z[['x','y','por']] )
# bandwidth, plus or minus 250 meters
bw = 500
# lags in 500 meter increments from zero to 10,000
hs = np.arange(0,10500,bw)
sv = SV( P, hs, bw )
plot( sv, sv, '.-' )
xlabel('Lag [m]')
ylabel('Semivariance')
title('Sample Semivariogram') ;
savefig('sample_semivariogram.png',fmt='png',dpi=200) # Modeling

Now that we’ve calculated the semivariogram, we will need to fit a model to the data. There are three popular models, the spherical, exponential, and the Gaussian. Here, we’ll implement the spherical model. First, we will present a function named opt() for determining the optimal value a for the spherical model.

def opt( fct, x, y, C0, parameterRange=None, meshSize=1000 ):
if parameterRange == None:
parameterRange = [ x, x[-1] ]
mse = np.zeros( meshSize )
a = np.linspace( parameterRange, parameterRange, meshSize )
for i in range( meshSize ):
mse[i] = np.mean( ( y - fct( x, a[i], C0 ) )**2.0 )
return a[ mse.argmin() ]


The opt() function finds the optimal parameter for fitting a spherical model to the semivariogram data. The spherical model is given by the function spherical(). On the last line we see that spherical() returns itself in a map() function, which seems odd. The idea is that the input h can be a single float value, or list or NumPy array of floats. If h is a single value, then line 9 is called. If h is a list or an array (an iterable) then line 17 is called, which applies line 9 to each value of h.

def spherical( h, a, C0 ):
'''
Spherical model of the semivariogram
'''
# if h is a single digit
if type(h) == np.float64:
# calculate the spherical function
if h <= a:
return C0*( 1.5*h/a - 0.5*(h/a)**3.0 )
else:
return C0
# if h is an iterable
else:
# calcualte the spherical function for all elements
a = np.ones( h.size ) * a
C0 = np.ones( h.size ) * C0
return map( spherical, h, a, C0 )


Next, cvmodel() fits a model to the semivariogram data and returns a covariance method named covfct().

def cvmodel( P, model, hs, bw ):
'''
Input:  (P)      ndarray, data
(model)  modeling function
- spherical
- exponential
- gaussian
(hs)     distances
(bw)     bandwidth
Output: (covfct) function modeling the covariance
'''
# calculate the semivariogram
sv = SV( P, hs, bw )
# calculate the sill
C0 = C( P, hs, bw )
# calculate the optimal parameters
param = opt( model, sv, sv, C0 )
# return a covariance function
covfct = lambda h, a=param: C0 - model( h, a, C0 )
return covfct


At this point we’ll plot our model and see if it represents our data well.

sp = cvmodel( P, model=spherical, hs=np.arange(0,10500,500), bw=500 )
plot( sv, sv, '.-' )
plot( sv, sp( sv ) ) ;
title('Spherical Model')
ylabel('Semivariance')
xlabel('Lag [m]')
savefig('semivariogram_model.png',fmt='png',dpi=200) # Simple Kriging

Now that we have a model for the semivariogram, we can write a function to perform the kriging. The fundamental relationship is a matrix equation,

(3) Here, is a matrix of covariances calculated using the spherical model, is a vector of simple kriging weights, and is the vector of covariances between the data points and an unsampled point. Our kriging function takes the data set P, the model, the distances hs, the bandwidth bw, the coordinates of the unsampled point u, and the number of surrounding points N to use in the calculation.

def krige( P, model, hs, bw, u, N ):
'''
Input  (P)     ndarray, data
(model) modeling function
- spherical
- exponential
- gaussian
(hs)    kriging distances
(bw)    kriging bandwidth
(u)     unsampled point
(N)     number of neighboring
points to consider
'''

# covariance function
covfct = cvmodel( P, model, hs, bw )
# mean of the variable
mu = np.mean( P[:,2] )

# distance between u and each data point in P
d = np.sqrt( ( P[:,0]-u )**2.0 + ( P[:,1]-u )**2.0 )
# add these distances to P
P = np.vstack(( P.T, d )).T
# sort P by these distances
# take the first N of them
P = P[d.argsort()[:N]]

# apply the covariance model to the distances
k = covfct( P[:,3] )
# cast as a matrix
k = np.matrix( k ).T

# form a matrix of distances between existing data points
K = squareform( pdist( P[:,:2] ) )
# apply the covariance model to these distances
K = covfct( K.ravel() )
# re-cast as a NumPy array -- thanks M.L.
K = np.array( K )
# reshape into an array
K = K.reshape(N,N)
# cast as a matrix
K = np.matrix( K )

# calculate the kriging weights
weights = np.linalg.inv( K ) * k
weights = np.array( weights )

# calculate the residuals
residuals = P[:,2] - mu

# calculate the estimation
estimation = np.dot( weights.T, residuals ) + mu

return float( estimation )


# Calculation

Here, we’ll calculate the kriging estimate at a number of unsampled points.

X0, X1 = P[:,0].min(), P[:,0].max()
Y0, Y1 = P[:,1].min(), P[:,1].max()
Z = np.zeros((80,100))
dx, dy = (X1-X0)/100.0, (Y1-Y0)/80.0
for i in range( 80 ):
print i,
for j in range( 100 ):
Z[i,j] = krige( P, model, hs, bw, (dy*j,dx*i), 16 )


# Plotting

Finally, for fun, we’ll try to reproduce one of the graphs found on the site with the data. I kind of guessed on the colormap.

cdict = {'red':   ((0.0, 1.0, 1.0),
(0.5, 225/255., 225/255. ),
(0.75, 0.141, 0.141 ),
(1.0, 0.0, 0.0)),
'green': ((0.0, 1.0, 1.0),
(0.5, 57/255., 57/255. ),
(0.75, 0.0, 0.0 ),
(1.0, 0.0, 0.0)),
'blue':  ((0.0, 0.376, 0.376),
(0.5, 198/255., 198/255. ),
(0.75, 1.0, 1.0 ),
(1.0, 0.0, 0.0)) }

my_cmap = matplotlib.colors.LinearSegmentedColormap('my_colormap',cdict,256)

fig, ax = subplots()
H = np.zeros_like( Z )
for i in range( Z.shape ):
for j in range( Z.shape ):
H[i,j] = np.round( Z[i,j]*3 )

ax.matshow( H, cmap=my_cmap, interpolation='nearest' )
ax.scatter( z.x/200.0, z.y/200.0, facecolor='none', linewidths=0.75, s=50 )
xlim(0,99) ; ylim(0,80)
xticks( [25,50,75], [5000,10000,15000] )
yticks( [25,50,75], [5000,10000,15000] )

savefig( 'krigingpurple.png', fmt='png', dpi=200 ) ## 27 thoughts on “Simple Kriging in Python”

1. awallin55 says:

https://wordpress.org/plugins/wp-syntax/

also whitespace matters in python, and at least my browser shows some additional unwanted line-breaks in the code-windows on this page 🙁

1. cjohnson318 says:

I have syntax highlighting, and all the code looks alright to me. What browser are you using?

2. manejandodatos (@manejandodatos) says:

Brilliant!

3. Patrick Li says:

great. Thank you.

4. Ngoc An says:

Could you please show me how to find the file ZoneA.dat?

1. cjohnson318 says:

Right, that was obscure. Go to http://people.ku.edu/~gbohling/geostats/index.html, and then click on My tutorial on reservoir modeling with WinGslib, with example data. This will download a ZIP file with the data sets in it.

5. luisblacuttluis says:

Hi, I’m having problems running the
fig, ax = subplots()
line

1. cjohnson318 says:

You might need to execute something like “from pylab import *”. I use the IPython notebook interface with the “–pylab=inline” flag at startup, and I forget that that imports a lot of plotting functionality automatically. Let me know if the “from pylab import *” line works for you.

1. luis says:

Hi, the instruction went through, however I get the following error message

Traceback (most recent call last):
File “<pyshell#10>”, line 1, in
ax.scatter( z.x, z.y, c=z.por, cmap=’gray’ )
AttributeError: ‘list’ object has no attribute ‘x’

I’m working on a 64-bit machine with OS Ubuntu 14.04

1. cjohnson318 says:

Okay, that says that (lowercase) z is a list, but at the top we declare (lowercase) z to be a pandas DataFrame. Using a DataFrame isn’t essential, I just like being able to name columns so I can have more readable code. If you want, you can email me your code and data and I can look at it. I’m connor (dot) labaume (dot) johnson (at gmail dot com).

6. martin says:

I have a problem with your script:

sp.sv (from sp = cvmodel( P, model=spherical, hs=np.arange(0,10500,500), bw=500 ))

so.sv
Traceback (most recent call last):
File “”, line 1, in
AttributeError: ‘function’ object has no attribute ‘sv’

For me, the correct formulation is :

sp(sv)

After:with model=spherical

sp.model( sp.sv ) ) is sp(model(sp.sv)) = sp(spherical(sp.sv))

TypeError: spherical() takes exactly 3 arguments (1 given)

effectively from def spherical( h, a, C0 ):

Thanks

1. cjohnson318 says:

Hi Martin, I simplified the code a few weeks ago from a more OOP approach to the present state in order to save time. (It didn’t save any time.) I accidentally forgot to update some of the examples. I’m apologize for the inconvenience. I also have code up at https://github.com/cjohnson318/geostatistics

1. Martin says:

Thanks,I know, but how to plot plot( sv, sp( sv ) ) ?

1. cjohnson318 says:

Oh, okay. Try adding “from pylab import *” before you plot.

1. Martin says:

Sorry, but It is not the problem of matplotlib, it is the problem of your new solution:

plot( sv, sp( sv ) )

The resulting curve is horizontally reflected

http://i.imgur.com/lnlfFfp.png

2. cjohnson318 says:

Oh, right. So, the model fits the empirical semivariogram, but the output is flipped because it is modeling the covariance described by the semivariogram. For small distances the semivariogram is low, but the covariance is high. For greater distances, the semivariogram is high, but the covariance is low.

7. Martin says:

Many thanks, it works now with spherical (an also with exponential and gaussian from your https://github.com/cjohnson318/geostatistics/blob/master/krige.py script)

A little correction in the function krige:
you need to put K = np.array(K) between the lines 29 and 31 ,

K = covfct( K.ravel() ) K =np.array(K) #reshape into an array K = K.reshape(N,N) 

otherwise:

K = K.reshape(N,N)
AttributeError: ‘list’ object has no attribute ‘reshape’

because originally: type(K) = list()

1. cjohnson318 says:

Excellent! Thanks! I’ll definitely fix that!

2. Javier Jofré says:

What happened with the Python code?

3. guido says:

i have the same problem of the flipped semivariogram, the link https://github.com/cjohnson318/geostatistics/blob/master/krige.py it’s broken, where can I find the files kriege.py? can you help me?

1. cjohnson318 says:

Hi, I’m sorry, try this link: https://github.com/cjohnson318/geostatsmodels. There’s some more tutorial material there also.

1. Guido says:

Hi, tank you, I have studied the intro to geostatistcs and variogram analysis, I may have found the problem, line 19 of cvmodel function
covfct = lambda h, a=param: C0 – model( h, a, C0 )
should be changed in
covfct = lambda h, a=param: model( h, a, C0 )
what do you think about? may be the solution?

1. cjohnson318 says:

I think the semivariogram and the covariance function are always flipped versions of each other. The semivariogram shows how the variance increases with distance, and the covariance function describes how covariance decreases with distance. You examine the semivariogram to decide on a model, and the sill and nugget parameters, but internally kriging uses a covariance function, which is a flipped version of the semivariogram.

8. Jesús Gómez says:

Great. Just to bring you back something: there are little changes to make all that code work on Python3:

the map call used in the spherical function definition should be wrapped with a list call: return list(map( spherical, h, a, C0 )). This is due Python3 map function returns a iterator and not a list directly.
The print statement when calculating the kriging estimate should be substituted with a print function call: print(i). In Python3, the print statement is gone and substituted with a print function.

Additionally, you should point that before “calculating the kriging estimate” the variable model must be set to any of the model functions: in this case, the spherical model. model = spherical.

Thank you.

9. Patrick Ng says:

A nice and easy to follow geostats example. Attaboy Mr. Johnson.

10. John says:

It is easy to understand, excellent explanation,
great tool. Thank you very much !!!