-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathKrigingExample.py
124 lines (75 loc) · 2.94 KB
/
KrigingExample.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
# -*- coding: utf-8 -*-
# <nbformat>3.0</nbformat>
# <markdowncell>
# This notebook is the reproduction of an exercise found at http://people.ku.edu/~gbohling/cpe940/Kriging.pdf
import pdb
#pdb.set_trace()
import utilities
import variograms
import model
import kriging
import geoplot
import pandas
import timeit
from matplotlib.pyplot import *
# <markdowncell>
# We'll read the data from `ZoneA.dat`.
# <codecell>
z = utilities.readGeoEAS('ZoneA.dat')
# <markdowncell>
# We want the first, second and fourth columns of the data set, representing the x and y spatial coordinates, and the porosity.
# <codecell>
P = z[:,[0,1,3]]
# <markdowncell>
# We'll be interested in determining the porosity at a point (2000,4700).
# <codecell>
pt = [ 2000, 4700 ]
# <markdowncell>
# We can plot our region of interest as follows:
# <codecell>
scatter( P[:,0], P[:,1], c=P[:,2], cmap=geoplot.YPcmap )
title('Zone A Subset % Porosity')
colorbar()
xmin, xmax = 0, 4250
ymin, ymax = 3200, 6250
xlim(xmin,xmax)
ylim(ymin,ymax)
for i in range( len( P[:,2] ) ):
x, y, por = P[i]
if( x < xmax )&( y > ymin )&( y < ymax ):
text( x+100, y, '{:4.2f}'.format( por ) )
scatter( pt[0], pt[1], marker='x', c='k' )
text( pt[0]+100 , pt[1], '?')
xlabel('Easting (m)')
ylabel('Northing (m)') ;
# <markdowncell>
# We can determine the parameters for our model by looking at the semivariogram and trying to determine the appropriate range and sill.
# <codecell>
tolerance = 250
lags = np.arange( tolerance, 10000, tolerance*2 )
sill = np.var( P[:,2] )
# <markdowncell>
# The semivariogram plotting function, `svplot()`, plots sill as a dashed line, and the empirical semivariogram as determined from the data. It optionally plots a semivariance model.
# <codecell>
geoplot.semivariogram( P, lags, tolerance )
# <markdowncell>
# We can pass a model to this function using the optional `model` argument and see it plotted in red.
# <codecell>
svm = model.semivariance( model.spherical, ( 4000, sill ) )
geoplot.semivariogram( P, lags, tolerance, model=svm )
# <markdowncell>
# The covariance modeling function function will return a spherical covariance model that takes a distance as input, and returns an covariance estimate. We've used the global variance of the porosity in `ZoneA.dat` as the sill.
# <codecell>
covfct = model.covariance( model.spherical, ( 4000, sill ) )
# <markdowncell>
# We can then krige the data, using the covariance model, the point we are interested in, (2000,47000), and `N=6` signifying that we only want to use the six nearest points. The output of the simple and ordinary kriging functions below is the krigin estimate, and the standard deviation of the kriging estimate.
# <codecell>
start = timeit.timeit()
kriging.simple( P, covfct, pt, N=16 )
end = timeit.timeit()
print end - start
# <codecell>
# <codecell>
#est, kstd = kriging.krige( P, covfct, [[2000,4700],[2100,4700],[2000,4800],[2100,4800]], 'simple', N=6 )
# <codecell>
# <codecell>