link管理

链接快照平台

输入网页链接，自动生成快照
标签化管理网页链接

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,

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

Here, $N$ is the number of data points, $\hat{\mu}$ is the sample mean, and $z_{k}$ 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 $h$ . 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[0] 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[0], sv[1], '.-' ) 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[1], x[-1] ] mse = np.zeros( meshSize ) a = np.linspace( parameterRange[0], parameterRange[1], 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[0], bw ) # calculate the optimal parameters param = opt( model, sv[0], sv[1], 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[0], sv[1], '.-' ) plot( sv[0], sp( sv[0] ) ) ; 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,

Here, $K$ is a matrix of covariances calculated using the spherical model, $\lambda$ is a vector of simple kriging weights, and $k$ 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[0] )**2.0 + ( P[:,1]-u[1] )**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[0] ): for j in range( Z.shape[1] ): 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 )

how about using a syntax-highlighting plugin!? makes the code more readable:
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 🙁

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

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).

After:with model=spherical

sp.model( sp.sv[0] ) ) is sp(model(sp.sv[0])) = sp(spherical(sp.sv[0]))

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

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

Thanks

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

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

plot( sv[0], sp( sv[0] ) )

The resulting curve is horizontally reflected

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

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.

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()
   
   
    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?
    

    thanks in advance
   
   
    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?
    

    thanks in advance
   
   
    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.
   
   
    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.