Kriging¶

In [1]:
%%javascript
IPython.OutputArea.prototype._should_scroll = function(lines) {
    return false;
}
In [2]:
import gstlearn as gl
import gstlearn.plot as gp
import matplotlib.pyplot as plt
import numpy as np
import os

Then we download the data base dat.

In [3]:
!wget -q https://soft.minesparis.psl.eu/gstlearn/data/Scotland/Scotland_Temperatures.NF
In [4]:
!wget -q https://soft.minesparis.psl.eu/gstlearn/data/Scotland/Scotland_Elevations.NF
In [5]:
fileNF = "Scotland_Temperatures.NF"
dat = gl.Db.createFromNF(fileNF)

Calculate the experimental variogram vario2dir (in 2 directions)

In [6]:
varioParamMulti = gl.VarioParam.createMultiple(2, 15, 15.)
vario2dir = gl.Vario(varioParamMulti, dat)
err = vario2dir.compute()

Calculate the fitted model fitmodOK (add the Universality Condition)

In [7]:
fitmodOK = gl.Model()
err = fitmodOK.fit(vario2dir,types=[gl.ECov.NUGGET, gl.ECov.EXPONENTIAL, gl.ECov.GAUSSIAN])
err = fitmodOK.addDrift(gl.Drift1())
ax = gp.varmod(vario2dir, fitmodOK)
No description has been provided for this image

Define the Unique Neighborhood unique.neigh.

In [8]:
uniqueNeigh = gl.NeighUnique.create()

Get the extension of the Data:

In [9]:
dat.getExtremas()
Out[9]:
array([[  78.2,  460.7],
       [ 530.4, 1208.9]])

Create the Target file grid.

In [10]:
grid = gl.DbGrid.create(x0=[65,535],dx=[4.94, 4.96],nx=[81,137])
dbfmt = gl.DbStringFormat.createFromFlags(flag_resume=False, flag_vars=True,
                                          flag_extend=True, flag_stats=False,
                                          flag_array=False)
grid.display(dbfmt)

Data Base Grid Characteristics
==============================

Data Base Extension
-------------------
Coor #1 - Min =     65.000 - Max =    460.200 - Ext = 395.2
Coor #2 - Min =    535.000 - Max =   1209.560 - Ext = 674.56

Variables
---------
Column = 0 - Name = rank - Locator = NA
Column = 1 - Name = x1 - Locator = x1
Column = 2 - Name = x2 - Locator = x2

Perform Ordinary Kriging

In [11]:
err = gl.kriging(dat,grid,fitmodOK,uniqueNeigh,namconv=gl.NamingConvention("OK"))

Plotting the results

In [12]:
gp.setDefaultGeographic(dims=[8,8])
ax = grid.plot()
ax.decoration(title="Ordinary Kriging over whole Grid")
No description has been provided for this image

Reading the Elevation Grid

In [13]:
fileNF = "Scotland_Elevations.NF"
grid = gl.DbGrid.createFromNF(fileNF)
grid.display(dbfmt)

Data Base Grid Characteristics
==============================

Data Base Extension
-------------------
Coor #1 - Min =     65.000 - Max =    455.123 - Ext = 390.123
Coor #2 - Min =    535.000 - Max =   1200.109 - Ext = 665.109

Variables
---------
Column = 0 - Name = Longitude - Locator = x1
Column = 1 - Name = Latitude - Locator = x2
Column = 2 - Name = Elevation - Locator = f1
Column = 3 - Name = inshore - Locator = sel

The output grid now contains the selection inshore. Estimation is restricted to the active cells only.

In [14]:
err = gl.kriging(dat,grid,fitmodOK,uniqueNeigh,namconv=gl.NamingConvention("OK"))

Plotting the Estimation

In [15]:
ax = grid.plot("OK*estim")
ax = dat.plot(name_size="January_temp")
ax.decoration(title="Estimation by Ordinary Kriging")
No description has been provided for this image

Plotting the Standard deviation of Estimation Error

In [16]:
ax = grid.plot("OK*stdev")
ax = dat.plot(name_size="January_temp")
ax.decoration(title="St. dev. by Ordinary Kriging")
No description has been provided for this image

The Model fitmodOK is first duplicated into fitmodSK. Then the Universality Condition is deleted.

In [17]:
fitmodSK = fitmodOK.clone()
err = fitmodSK.delDrift(0)
err = fitmodSK.setMean(20.)

Simple Kriging is performed

In [18]:
err = gl.kriging(dat,grid,fitmodSK,uniqueNeigh, namconv=gl.NamingConvention("SK"))

Plotting the Estimation

In [19]:
ax = grid.plot("SK*estim")
ax = dat.plot(name_size="January_temp")
ax.decoration(title="Estimation by Simple Kriging")
No description has been provided for this image

Plotting the Standard deviation of the Estimation error

In [20]:
ax = grid.plot("SK*stdev")
ax = dat.plot(name_size="January_temp",flagCst=True)
ax.decoration(title="St. dev. by Simple Kriging")
No description has been provided for this image

Comparing Ordinary and Simple Kriging - Estimations

In [21]:
ax = gp.correlation(grid,namex="OK*estim",namey="SK*estim", bissLine=True, bins=100, flagSameAxes=True)
ax.decoration(title="Estimation Simple vs. Ordinary", 
              xlabel="Ordinary Kriging", ylabel="Simple Kriging")
No description has been provided for this image

Comparing Ordinary and Simple Kriging - St. dev.

In [22]:
ax = gp.correlation(grid,namex = "OK*stdev",namey="SK*stdev", bissLine=True, bins=100,flagSameAxes=True)
ax.decoration(title="St. dev. Simple vs. Ordinary", 
              xlabel="Ordinary Kriging", ylabel="Simple Kriging")
No description has been provided for this image

Cross-Validation

In [23]:
err = gl.xvalid(dat,fitmodOK,uniqueNeigh, 
                namconv=gl.NamingConvention("Xvalid", True, True, False, gl.ELoc.UNKNOWN))

Cross-validation - Histogram of Errors

In [24]:
ax = gp.histogram(dat,name="*esterr*",bins=30,fill="blue")
ax.decoration(xlabel="Estimation Errors",title="Cross-Validation")
No description has been provided for this image

Cross-validation - Histogram of Standardized Errors

In [25]:
ax = gp.histogram(dat,name="*stderr*",bins=30,fill="blue")
ax.decoration(xlabel="Standardized Errors", title="Cross-Validation")
No description has been provided for this image

Cross-validation Scores

In [26]:
print(round(np.nanmean(dat.getColumn("*esterr*")),4))
print(round(np.nanmean(np.square(dat.getColumn("*esterr*"))),4))
print(round(np.nanmean(np.square(dat.getColumn("*stderr*"))),4))

-0.0042
0.2394
0.9118

Display Errors of the Cross Validation

In [27]:
ax = grid.plot("inshore")
ax = dat.plot(name_size="*esterr",sizmax=300)
ax.decoration(title="Cross-Validation scores")
No description has been provided for this image

Display Errors of the Cross Validation (in absolute value)

In [28]:
ax = grid.plot("inshore")
ax = dat.plot(name_size="*esterr",sizmax=300,flagAbsSize=True)
ax.decoration(title="Cross-Validation scores (abs. value)")
No description has been provided for this image

We design a small Moving Neighborhood small.neigh with only 1 sample per neighborhood.

In [29]:
smallNeigh = gl.NeighMoving.create(flag_xvalid=False, nmini=1, nmaxi=1, radius=1000000)

We perform Ordinary Kriging

In [30]:
err = gl.kriging(dat,grid,fitmodOK,smallNeigh,namconv=gl.NamingConvention("Small"))

Graphic representation

In [31]:
ax = grid.plot("Small*estim")
ax = dat.plot(name_size="January_temp")
ax.decoration(title="Estimation by Ordinary Kriging (Small Moving Neigh.)")
No description has been provided for this image

Building a reasonable Moving Neighborhood, although with a limited extension (radius)

In [32]:
movingNeigh = gl.NeighMoving.create(nmini=1, nmaxi=10, radius=20)

Running the Ordinary Kriging

In [33]:
err = gl.kriging(dat,grid,fitmodOK,movingNeigh, namconv=gl.NamingConvention("Reduced"))

Plotting the results

In [34]:
ax = grid.plot("Reduced*estim")
ax = dat.plot(name_size="January_temp")
ax.decoration(title="Estimation by Ordinary Kriging (Reduced Moving Neigh.)")
No description has been provided for this image

Lots of target sites are not estimated as no sample is found within the neighborhood.

Building a reasonable Moving Neighborhood correctly tuned: 10 samples (maximum) selected in a radius of 150 around the target site.

In [35]:
movingNeigh = gl.NeighMoving.create(nmini=1, nmaxi=10, radius=150)

Running the Ordinary Kriging

In [36]:
err = gl.kriging(dat,grid,fitmodOK,movingNeigh, namconv=gl.NamingConvention("Moving"))

Graphic representation

In [37]:
ax = grid.plot("Moving*estim")
ax = dat.plot(name_size="January_temp")
ax.decoration(title="Estimation by Ordinary Kriging (Moving Neigh.)")
No description has been provided for this image

For the standard deviation of Estimation error

In [38]:
ax = grid.plot("Moving*stdev")
ax = dat.plot(name_size="January_temp")
ax.decoration(title="St. dev. by Ordinary Kriging (Moving Neigh.)")
No description has been provided for this image

Comparing Unique and Moving Neighborhoods: Estimations

In [39]:
ax = gp.correlation(grid,namex = "OK*estim",namey="Moving*estim", bins=100, bissLine=True, flagSameAxes=True)
ax.decoration(title="Comapring Estimations", xlabel="Unique", ylabel="Moving")
No description has been provided for this image

Comparing Unique and Moving Neighborhoods: Standard deviations

In [40]:
ax = gp.correlation(grid,namex = "OK*stdev",namey="Moving*stdev", bins=100, bissLine=True, flagSameAxes=True)
ax.decoration(title="Comparing St. Deviations", xlabel="Unique", ylabel="Moving")
No description has been provided for this image