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
"""


