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 [ ]: