Variography¶
In this preamble, we load the gstlearn library.
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 the necessary data set is downloaded and named dat: the target variable is January_temp
In [3]:
!wget -q https://soft.minesparis.psl.eu/gstlearn/data/Scotland/Scotland_Temperatures.NF
In [4]:
fileNF = "Scotland_Temperatures.NF"
dat = gl.Db.createFromNF(fileNF)
Experimental variogram¶
Data is a \alert{regionalized} variable
$$z_i = z(x_i)$$
The experimental variogram is a (dicrete) function:
$$\gamma(h)=\frac{1}{2N(h)}\sum_{i=1}^{N(h)}[z(x_i+h)-z(x_i)]^2$$
where $N(h)$ is the number of pairs of points distant by $h$
Variogram Cloud¶
In [5]:
varioParamOmni = gl.VarioParam.createOmniDirection(100)
grid_cloud = gl.db_variogram_cloud(dat, varioParamOmni)
grid_cloud
Out[5]:
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 = 10000 Grid characteristics: --------------------- Origin : 0.000 0.000 Mesh : 7.789 0.068 Number : 100 100 Variables --------- Column = 0 - Name = rank - Locator = NA Column = 1 - Name = x1 - Locator = x1 Column = 2 - Name = x2 - Locator = x2 Column = 3 - Name = Cloud.January_temp - Locator = NA
Experimental Variogram¶
We calculate the omni-directional variogram of the temperatures.
In [6]:
varioParamOmni = gl.VarioParam.createOmniDirection(40, 10)
varioexp = gl.Vario(varioParamOmni, dat)
err = varioexp.compute()
Print the variogram contents
In [7]:
varioexp
Out[7]:
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 = 40
Direction coefficients = 1.000 0.000
Direction angles (degrees) = 0.000 0.000
Tolerance on direction = 90.000 (degrees)
Calculation lag = 10.000
Tolerance on distance = 50.000 (Percent of the lag value)
For variable 1
Rank Npairs Distance Value
0 19.000 3.118 0.042
1 89.000 10.690 0.172
2 168.000 20.346 0.280
3 244.000 30.323 0.456
4 316.000 40.429 0.459
5 385.000 50.162 0.660
6 399.000 60.296 0.729
7 463.000 70.062 0.874
8 450.000 79.807 0.800
9 473.000 90.115 0.977
10 549.000 100.141 0.879
11 484.000 109.866 0.970
12 491.000 119.865 1.082
13 458.000 130.001 1.074
14 474.000 139.888 0.983
15 472.000 149.918 1.086
16 390.000 159.896 1.218
17 410.000 169.913 1.303
18 385.000 179.840 1.212
19 425.000 189.849 1.079
20 354.000 200.117 1.168
21 341.000 209.881 1.229
22 325.000 219.710 1.497
23 296.000 229.902 1.597
24 277.000 239.927 1.419
25 236.000 250.111 1.353
26 201.000 259.921 1.449
27 177.000 270.161 1.335
28 166.000 279.560 1.245
29 177.000 289.873 1.089
30 148.000 299.866 0.992
31 116.000 310.335 1.130
32 95.000 319.854 1.161
33 92.000 329.857 1.330
34 64.000 339.936 1.151
35 57.000 350.172 1.551
36 61.000 359.473 1.395
37 51.000 369.274 1.350
38 33.000 380.087 0.847
39 31.000 389.553 1.014
Plot the omni-directional variogram
In [8]:
gp.setDefault(dims=[6,6])
ax = gp.varmod(varioexp)
Represent the variogram with the number of pairs
In [9]:
ax = gp.varmod(varioexp,show_pairs=True)
Automatic Fitting procedure
In [10]:
fitmod = gl.Model()
err = fitmod.fit(varioexp)
ax = gp.varmod(varioexp, fitmod)
Print the Model
In [11]:
fitmod
Out[11]:
Model characteristics ===================== Space dimension = 2 Number of variable(s) = 1 Number of basic structure(s) = 1 Number of drift function(s) = 0 Number of drift equation(s) = 0 Covariance Part --------------- Spherical - Sill = 1.155 - Range = 135.117 Total Sill = 1.155
List of Basic structures
In [12]:
gl.ECov.printAll()
-2 - UNKNOWN : Unknown covariance
-1 - FUNCTION : External covariance function
0 - NUGGET : Nugget effect
1 - EXPONENTIAL : Exponential
2 - SPHERICAL : Spherical
3 - GAUSSIAN : Gaussian
4 - CUBIC : Cubic
5 - SINCARD : Sine Cardinal
6 - BESSEL_J : Bessel J
7 - BESSEL_K : Bessel K
8 - GAMMA : Gamma
9 - CAUCHY : Cauchy
10 - STABLE : Stable
11 - LINEAR : Linear
12 - POWER : Power
13 - ORDER1_GC : First Order Generalized covariance
14 - SPLINE_GC : Spline Generalized covariance
15 - ORDER3_GC : Third Order Generalized covariance
16 - ORDER5_GC : Fifth Order Generalized covariance
17 - COSINUS : Cosine
18 - TRIANGLE : Triangle
19 - COSEXP : Cosine Exponential
20 - REG1D : 1-D Regular
21 - PENTA : Pentamodel
22 - SPLINE2_GC : Order-2 Spline
23 - STORKEY : Storkey covariance in 1-D
24 - WENDLAND0 : Wendland covariance (2,0)
25 - WENDLAND1 : Wendland covariance (3,1)
26 - WENDLAND2 : Wendland covariance (4,2)
27 - MARKOV : Markovian covariances
Automatic Fitting (with given basic structures)
In [13]:
types=[gl.ECov.NUGGET, gl.ECov.CUBIC, gl.ECov.SPHERICAL]
err = fitmod.fit(varioexp, types=types)
ax = gp.varmod(varioexp, fitmod)
The resulting Model
In [14]:
fitmod
Out[14]:
Model characteristics ===================== Space dimension = 2 Number of variable(s) = 1 Number of basic structure(s) = 2 Number of drift function(s) = 0 Number of drift equation(s) = 0 Covariance Part --------------- Cubic - Sill = 0.413 - Range = 75.995 Spherical - Sill = 0.893 - Range = 240.635 Total Sill = 1.305
Model Fitting with Inequality constraints
In [15]:
constraints = gl.Constraints()
err = constraints.addItemFromParamId(gl.EConsElem.RANGE,1,0,0,gl.EConsType.UPPER,20.)
err = constraints.addItemFromParamId(gl.EConsElem.SILL,1,0,0,gl.EConsType.LOWER,0.03)
err = fitmod.fit(varioexp, types, constraints, gl.Option_VarioFit(True))
ax = gp.varmod(varioexp, fitmod)
The resulting Model
In [16]:
fitmod
Out[16]:
Model characteristics ===================== Space dimension = 2 Number of variable(s) = 1 Number of basic structure(s) = 3 Number of drift function(s) = 0 Number of drift equation(s) = 0 Covariance Part --------------- Nugget Effect - Sill = 0.000 Cubic - Sill = 0.109 - Range = 20.000 Spherical - Sill = 1.056 - Range = 155.372 Total Sill = 1.165
Model Fitting with Equality constraints
In [17]:
constraints = gl.Constraints()
err = constraints.addItemFromParamId(gl.EConsElem.RANGE,1,0,0,gl.EConsType.EQUAL,1000.)
err = constraints.addItemFromParamId(gl.EConsElem.SILL,1,0,0,gl.EConsType.EQUAL,0.4)
err = fitmod.fit(varioexp, types, constraints, gl.Option_VarioFit(True))
ax = gp.varmod(varioexp, fitmod)
The resulting Model
In [18]:
fitmod
Out[18]:
Model characteristics ===================== Space dimension = 2 Number of variable(s) = 1 Number of basic structure(s) = 3 Number of drift function(s) = 0 Number of drift equation(s) = 0 Covariance Part --------------- Nugget Effect - Sill = 0.053 Cubic - Sill = 0.400 - Range = 1000.000 Spherical - Sill = 1.003 - Range = 130.078 Total Sill = 1.457
Directional Variograms
In [19]:
varioParamMulti = gl.VarioParam.createMultiple(4, 15, 15.)
vario_4dir = gl.Vario(varioParamMulti, dat)
err = vario_4dir.compute()
ax = gp.varmod(vario_4dir)
Fitting a Multi-directional variogram
In [20]:
model_4dir = gl.Model()
err = model_4dir.fit(vario_4dir,types=types)
ax = gp.varmod(vario_4dir, model_4dir)
Calculating Variogram Map
In [21]:
grid_vmap = gl.db_vmap_compute(dat, gl.ECalcVario.VARIOGRAM)
ax = grid_vmap.plot()
Automatic Model Fitting from Variogram Map
In [22]:
modelVM = gl.Model()
err = modelVM.fitFromVMap(grid_vmap, types=types)
modelVM
Out[22]:
Model characteristics
=====================
Space dimension = 2
Number of variable(s) = 1
Number of basic structure(s) = 2
Number of drift function(s) = 0
Number of drift equation(s) = 0
Covariance Part
---------------
Nugget Effect
- Sill = 0.257
Cubic
- Sill = 0.943
- Ranges = 158.007 214.374
- Angles = 334.704 0.000
- Rotation Matrix
[, 0] [, 1]
[ 0,] 0.904 0.427
[ 1,] -0.427 0.904
Total Sill = 1.200
Drawing the Fitted Model as a Variogram Map
In [23]:
err = gl.dbgrid_model(grid_vmap, modelVM)
ax = grid_vmap.plot()
Compare Directional Variograms and Variogram Map
In [24]:
ax = gp.varmod(vario_4dir, modelVM)
In [ ]: