<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 This notebook is the reproduction of an exercise found at http
Step1:  We'll read the data from ZoneA.dat.
Step2:  We want the first, second and fourth columns of the data set, representing the x and y spatial coordinates, and the porosity.
Step3:  We'll be interested in determining the porosity at a point (2000,4700).
Step4:  We can plot our region of interest as follows
Step5:  We can determine the parameters for our model by looking at the semivariogram and trying to determine the appropriate range and sill.
Step6:  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.
Step7:  We can pass a model to this function using the optional model argument and see it plotted in red.
Step8:  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.
Step9:  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.
<ASSISTANT_TASK:>
Python Code:
import sys
sys.path.append('..')
sys.path.append('../geostatsmodels')
from geostatsmodels import utilities, variograms, model, kriging, geoplot
import matplotlib.pyplot as plt
import numpy as np
import pandas
Explanation: This notebook is the reproduction of an exercise found at http://people.ku.edu/~gbohling/cpe940/Kriging.pdf
End of explanation
z = utilities.readGeoEAS('../data/ZoneA.dat')
Explanation: We'll read the data from ZoneA.dat.
End of explanation
P = z[:,[0,1,3]]
Explanation: We want the first, second and fourth columns of the data set, representing the x and y spatial coordinates, and the porosity.
End of explanation
pt = [2000, 4700]
Explanation: We'll be interested in determining the porosity at a point (2000,4700).
End of explanation
plt.scatter(P[:,0], P[:,1], c=P[:,2], cmap=geoplot.YPcmap)
plt.title('Zone A Subset % Porosity')
plt.colorbar()
xmin, xmax = 0, 4250
ymin, ymax = 3200, 6250
plt.xlim(xmin,xmax)
plt.ylim(ymin,ymax)
for i in range(len(P[:,2])):
    x, y, por = P[i]
    if (x < xmax) & (y > ymin) & (y < ymax):
        plt.text( x+100, y, '{:4.2f}'.format( por ) ) 
plt.scatter(pt[0], pt[1], marker='x', c='k')
plt.text(pt[0] + 100 , pt[1], '?')
plt.xlabel('Easting (m)')
plt.ylabel('Northing (m)');
Explanation: We can plot our region of interest as follows:
End of explanation
tolerance = 250
lags = np.arange(tolerance, 10000, tolerance*2)
sill = np.var(P[:,2])
Explanation: We can determine the parameters for our model by looking at the semivariogram and trying to determine the appropriate range and sill.
End of explanation
geoplot.semivariogram(P, lags, tolerance)
Explanation: 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.
End of explanation
svm = model.semivariance(model.spherical, (4000, sill))
geoplot.semivariogram(P, lags, tolerance, model=svm)
Explanation: We can pass a model to this function using the optional model argument and see it plotted in red.
End of explanation
covfct = model.covariance(model.spherical, (4000, sill))
Explanation: 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.
End of explanation
kriging.simple(P, covfct, pt, N=6)
kriging.ordinary(P, covfct, pt, N=6)
est, kstd = kriging.krige(P, covfct, [[2000,4700],[2100,4700],[2000,4800],[2100,4800]], 'simple', N=6)
est
kstd
Explanation: 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.
End of explanation
<END_TASK>