Python version of test_SPDEDRIFT.cpp¶
In [1]:
%%javascript
IPython.OutputArea.prototype._should_scroll = function(lines) {
return false;
}
In [2]:
import gstlearn as gl
import matplotlib.pyplot as plt
import gstlearn.plot as gp
import os
verbose = False
flagSPDE = True
ndim = 2
gp.setDefaultGeographic(dims=[7,7], aspect=1)
In [3]:
!wget -q https://soft.minesparis.psl.eu/gstlearn/data/Scotland/temperatures.ascii
In [4]:
!wget -q https://soft.minesparis.psl.eu/gstlearn/data/Scotland/temperatures.asci
In [5]:
!wget -q https://soft.minesparis.psl.eu/gstlearn/data/Scotland/grid.ascii
In [6]:
filename = "temperatures.ascii"
temperatures = gl.Db.createFromNF(filename,verbose)
temperatures.setLocator("January_temp", gl.ELoc.Z)
temperatures.display()
ax = temperatures.plot("January_temp")
Data Base Characteristics ========================= Data Base Summary ----------------- File is organized as a set of isolated points Space dimension = 2 Number of Columns = 5 Maximum Number of UIDs = 5 Total number of samples = 151 Number of active samples = 151 Variables --------- Column = 0 - Name = Longitude - Locator = x1 Column = 1 - Name = Latitude - Locator = x2 Column = 2 - Name = Elevation - Locator = f1 Column = 3 - Name = January_temp - Locator = z1 Column = 4 - Name = sel - Locator = sel
In [7]:
filename = "grid.ascii"
grid = gl.DbGrid.createFromNF(filename,verbose)
grid.display()
ax = grid.plot("Elevation")
Data Base Grid Characteristics ============================== Data Base Summary ----------------- File is organized as a regular grid Space dimension = 2 Number of Columns = 4 Maximum Number of UIDs = 4 Total number of samples = 11097 Number of active samples = 3092 Grid characteristics: --------------------- Origin : 65.000 535.000 Mesh : 4.938 4.963 Number : 81 137 Variables --------- Column = 0 - Name = x1 - Locator = x1 Column = 1 - Name = x2 - Locator = x2 Column = 2 - Name = Elevation - Locator = f1 Column = 3 - Name = inshore - Locator = sel
Calculate the variogram
In [8]:
v = gl.VarioParam()
dir1 = gl.DirParam(18,25)
v.addDir(dir1)
vario = gl.Vario(v,temperatures)
md = gl.Model.createFromDb(temperatures)
md.addDrift(gl.Drift1(md.getContext()))
md.addDrift(gl.DriftF(md.getContext()))
vario.compute(model=md)
vario.display()
ax = gp.variogram(vario)
Variogram characteristics
=========================
Number of variable(s) = 1
Number of direction(s) = 1
Space dimension = 2
Variance-Covariance Matrix 0.363
Direction #1
------------
Number of lags = 18
Direction coefficients = 1.000 0.000
Direction angles (degrees) = 0.000 0.000
Tolerance on direction = 90.000 (degrees)
Calculation lag = 25.000
Tolerance on distance = 50.000 (Percent of the lag value)
For variable 1
Rank Npairs Distance Value
0 78.000 7.681 0.102
1 512.000 26.905 0.163
2 915.000 50.556 0.211
3 1135.000 74.820 0.264
4 1285.000 100.070 0.265
5 1190.000 124.927 0.321
6 1117.000 149.550 0.363
7 1004.000 175.060 0.444
8 924.000 199.603 0.466
9 769.000 224.493 0.547
10 594.000 249.105 0.561
11 438.000 274.605 0.582
12 363.000 298.685 0.472
13 236.000 323.920 0.565
14 153.000 349.725 0.424
15 105.000 373.088 0.345
16 76.000 399.464 0.410
17 58.000 426.095 0.411
Load the (same) variogram from a neutral file
In [9]:
!wget -q https://soft.minesparis.psl.eu/gstlearn/data/Scotland/vario.ascii
In [10]:
!wget -q https://soft.minesparis.psl.eu/gstlearn/data/Scotland/model.ascii
In [11]:
filename = "vario.ascii"
vario = gl.Vario.createFromNF(filename,verbose)
vario.display()
ax = gp.variogram(vario)
Variogram characteristics
=========================
Number of variable(s) = 1
Number of direction(s) = 1
Space dimension = 2
Variance-Covariance Matrix 0.363
Direction #1
------------
Number of lags = 18
Direction coefficients = 1.000 0.000
Direction angles (degrees) = 0.000 0.000
Tolerance on direction = 90.000 (degrees)
Calculation lag = 25.000
Tolerance on distance = 50.000 (Percent of the lag value)
For variable 1
Rank Npairs Distance Value
0 78.000 7.681 0.102
1 512.000 26.905 0.163
2 915.000 50.556 0.211
3 1135.000 74.820 0.264
4 1285.000 100.070 0.265
5 1190.000 124.927 0.321
6 1117.000 149.550 0.363
7 1004.000 175.060 0.444
8 924.000 199.603 0.466
9 769.000 224.493 0.547
10 594.000 249.105 0.561
11 438.000 274.605 0.582
12 363.000 298.685 0.472
13 236.000 323.920 0.565
14 153.000 349.725 0.424
15 105.000 373.088 0.345
16 76.000 399.464 0.410
17 58.000 426.095 0.411
Fit the variogram of residuals in a model having drifts
In [12]:
model = gl.Model.createFromDb(temperatures)
model.setDriftList(md.getDriftList())
structs = [gl.ECov.NUGGET,gl.ECov.BESSEL_K,gl.ECov.BESSEL_K]
structs = [gl.ECov.NUGGET,gl.ECov.BESSEL_K]
consNug = gl.ConsItem.define(gl.EConsElem.SILL,0, type = gl.EConsType.UPPER,value = 0.1)
cons1P = gl.ConsItem.define(gl.EConsElem.PARAM,1, type = gl.EConsType.EQUAL,value = 1)
cons1Rm = gl.ConsItem.define(gl.EConsElem.RANGE,1, type = gl.EConsType.LOWER,value = 100)
cons1RM = gl.ConsItem.define(gl.EConsElem.RANGE,1, type = gl.EConsType.UPPER,value = 350)
cons2P = gl.ConsItem.define(gl.EConsElem.PARAM,2, type = gl.EConsType.EQUAL,value = 2)
cons2Rm = gl.ConsItem.define(gl.EConsElem.RANGE,2, type = gl.EConsType.LOWER,value = 100)
cons2RM = gl.ConsItem.define(gl.EConsElem.RANGE,2, type = gl.EConsType.UPPER,value = 400)
a = gl.Constraints()
a.addItem(consNug)
a.addItem(cons1P)
a.addItem(cons1Rm)
a.addItem(cons1RM)
a.addItem(cons2P)
a.addItem(cons2Rm)
a.addItem(cons2RM)
err = model.fit(vario,structs,constraints=a)
model.display()
ax = gp.varmod(vario,model)
ax.decoration(title="Vario of the residuals")
Model characteristics ===================== Space dimension = 2 Number of variable(s) = 1 Number of basic structure(s) = 2 Number of drift function(s) = 2 Number of drift equation(s) = 2 Covariance Part --------------- Nugget Effect - Sill = 0.100 K-Bessel (Third Parameter = 1) - Sill = 0.420 - Range = 282.602 - Theo. Range = 81.580 Total Sill = 0.520 Drift Part ---------- Universality Condition External Drift - Rank=0
Reload the (same) model from a neutral file
In [13]:
filename = "model.ascii"
model = gl.Model.createFromNF(filename,verbose)
model.display()
ax = gp.varmod(vario,model)
ax.decoration(title="Vario of the residuals")
Model characteristics ===================== Space dimension = 2 Number of variable(s) = 1 Number of basic structure(s) = 2 Number of drift function(s) = 2 Number of drift equation(s) = 2 Covariance Part --------------- Nugget Effect - Sill = 0.100 K-Bessel (Third Parameter = 1) - Sill = 0.420 - Range = 282.574 - Theo. Range = 81.572 Total Sill = 0.520 Drift Part ---------- Universality Condition External Drift - Rank=0
In [14]:
spde = gl.SPDE(model,grid,temperatures,gl.ESPDECalcMode.KRIGING)
coeffs = spde.getCoeffs()
ax = gp.correlation(temperatures, namex="Elevation", namey="*temp", asPoint=True)
if len(coeffs)>1:
plt.plot([0,400], [coeffs[0],coeffs[0]+coeffs[1]*400])
In [15]:
spde.compute()
spde.query(grid)
ax = grid.plot(usesel=True)
In [16]:
neighU = gl.NeighUnique.create();
gl.kriging(temperatures, grid, model, neighU);
ax = grid.plot(usesel=True)
In [17]:
ax = gp.correlation(grid,namex="*estim",namey="spde*",asPoint=True, diagLine=True)
In [ ]: