SPDE¶
This document aims at demonstrating the use of SPDE for performing Estimation. It is based on Scotland Data
In [1]:
import gstlearn as gl
import gstlearn.plot as gp
import gstlearn.document as gdoc
import matplotlib.pyplot as plt
import os
import numpy as np
gdoc.setNoScroll()
Getting the Data Bases from the official website.
In [2]:
temp_nf = gdoc.loadData("Scotland", "Scotland_Temperatures.NF")
dat = gl.Db.createFromNF(temp_nf)
The input Data Base (called temperatures) contains the target variable (January_temp). Note that this data base also contains the elevation variable (called Elevation) which is assigned the locator f (for external drift). We finally add a selection (called sel) which only consider the samples where the temperature is actually calculated.
In [3]:
temperatures = gl.Db.createFromNF(temp_nf)
temperatures.setLocator("January_temp", gl.ELoc.Z)
temperatures.setLocator("Elevation",gl.ELoc.F)
iuid = temperatures.addSelection(np.invert(np.isnan(temperatures["J*"])),"sel")
temperatures.display()
ax = temperatures.plot("January_temp")
ax.decoration(title="Temperature data")
Data Base Characteristics ========================= Data Base Summary ----------------- File is organized as a set of isolated points Space dimension = 2 Number of Columns = 6 Total number of samples = 236 Number of active samples = 151 Variables --------- Column = 0 - Name = rank - Locator = NA Column = 1 - Name = Longitude - Locator = x1 Column = 2 - Name = Latitude - Locator = x2 Column = 3 - Name = Elevation - Locator = f1 Column = 4 - Name = January_temp - Locator = z1 Column = 5 - Name = sel - Locator = sel
The output file is a grid (called grid). It contains an active selection (inshore) as well as the elevation over the field (called Elevation). Note that this variable is assigned the locator f (external drift).
In [4]:
elev_nf = gdoc.loadData("Scotland", "Scotland_Elevations.NF")
grid = gl.DbGrid.createFromNF(elev_nf)
grid.display()
ax = grid.plot("Elevation")
ax.decoration(title="Elevation")
Data Base Grid Characteristics ============================== Data Base Summary ----------------- File is organized as a regular grid Space dimension = 2 Number of Columns = 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 = Longitude - Locator = x1 Column = 1 - Name = Latitude - Locator = x2 Column = 2 - Name = Elevation - Locator = f1 Column = 3 - Name = inshore - Locator = sel
Calculate the omni-directional variogram of the temperature for 18 lags of 25.
In [5]:
vparam = gl.VarioParam.createOmniDirection(npas=18, dpas=25)
vario = gl.Vario(vparam)
vario.compute(temperatures)
vario.display()
ax = gp.variogram(vario)
ax.decoration(title="Variogram of Temperature (Raw)")
Variogram characteristics
=========================
Number of variable(s) = 1
Number of direction(s) = 1
Space dimension = 2
Variance-Covariance Matrix 1.020
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.137
1 512.000 26.905 0.376
2 915.000 50.556 0.614
3 1135.000 74.820 0.855
4 1285.000 100.070 0.935
5 1190.000 124.927 1.050
6 1117.000 149.550 1.099
7 1004.000 175.060 1.218
8 924.000 199.603 1.175
9 769.000 224.493 1.502
10 594.000 249.105 1.377
11 438.000 274.605 1.316
12 363.000 298.685 1.071
13 236.000 323.920 1.190
14 153.000 349.725 1.321
15 105.000 373.088 1.219
16 76.000 399.464 0.802
17 58.000 426.095 0.677
We calculate the variogram (using the same calculation parameters) based on the residuals after a trend has been removed. This trend is considered as a linear combination of the external drift information.
In [6]:
vparam = gl.VarioParam.createOmniDirection(npas=18, dpas=25)
vario = gl.Vario(vparam)
md = gl.Model()
md.setDriftIRF(order=0, nfex=1)
vario.compute(temperatures,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
Fit the variogram of residuals in a model having drifts. Some constraints have been added during the fitting step.
In [7]:
model = md
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)
a = gl.Constraints()
a.addItem(consNug)
a.addItem(cons1P)
a.addItem(cons1Rm)
a.addItem(cons1RM)
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.583 - Theo. Range = 81.575 Total Sill = 0.520 Drift Part ---------- Universality_Condition External_Drift:0
Derive the parameter of a Global Trend (composed of the Elevation as a drift function) using the SPDE formalism.
In [8]:
spde = gl.SPDE(model,grid,temperatures,gl.ESPDECalcMode.KRIGING)
coeffs = spde.getCoeffs()
print("Trend coefficients:", coeffs)
Trend coefficients: [ 3.97004997 -0.00659772]
Represent the scatter plot of the Temperature given the Elevation and add the Global Trend (calculated beforehand)
In [9]:
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])
We perform the Estimation in the SPDE framework (considering the Elevation as the Global Trend)
In [10]:
err = gl.krigingSPDE(temperatures, grid, model)
ax = grid.plot(useSel=True)
ax.decoration(title="Temperature (using Elevation as global Trend)")
We also perform the Estimation by Kriging (using Elevation as External Drift). This estimation is performed in Unique Neighborhood.
In [11]:
neighU = gl.NeighUnique.create();
gl.kriging(temperatures, grid, model, neighU);
ax = grid.plot(useSel=True)
ax.decoration(title="Temperature (with Elevation as External Drift)")
Comparing both estimates
In [12]:
grid
Out[12]:
Data Base Grid Characteristics ============================== Data Base Summary ----------------- File is organized as a regular grid Space dimension = 2 Number of Columns = 7 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 = Longitude - Locator = x1 Column = 1 - Name = Latitude - Locator = x2 Column = 2 - Name = Elevation - Locator = f1 Column = 3 - Name = inshore - Locator = sel Column = 4 - Name = KrigingSPDE.January_temp.estim - Locator = NA Column = 5 - Name = Kriging.January_temp.estim - Locator = z1 Column = 6 - Name = Kriging.January_temp.stdev - Locator = NA
In [13]:
ax = gp.correlation(grid,namex="Kriging.January_temp.estim",namey="KrigingSPDE.January_temp.estim",
asPoint=True, diagLine=True)
ax.decoration(xlabel="Kriging with External Drift", ylabel="SPDE Kriging")
