In [5]:
ax = grid_cloud.plot("Cloud*")
plt.gca().set_aspect('100')
ax.decoration(title="Variogram Cloud")
In this preamble, we load the gstlearn library.
%%javascript
IPython.OutputArea.prototype._should_scroll = function(lines) {
return false;
}
import gstlearn as gl
import gstlearn.plot as gp
import matplotlib.pyplot as plt
import numpy as np
import os
import urllib.request
Then the necessary data set is downloaded and named dat: the target variable is January_temp
url = 'https://soft.minesparis.psl.eu/gstlearn/data/Scotland/Scotland_Temperatures.NF'
temp_nf, head = urllib.request.urlretrieve(url)
dat = gl.Db.createFromNF(temp_nf)
Data is a regionalized variable
$$z_i = z(x_i)$$The experimental variogram is a (discrete) 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$
varioParamOmni = gl.VarioParam.createOmniDirection(100)
grid_cloud = gl.db_variogram_cloud(dat, varioParamOmni)
grid_cloud
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
ax = grid_cloud.plot("Cloud*")
plt.gca().set_aspect('100')
ax.decoration(title="Variogram Cloud")
We calculate the omni-directional variogram of the temperatures.
varioParamOmni = gl.VarioParam.createOmniDirection(40, 10)
varioexp = gl.Vario(varioParamOmni, dat)
err = varioexp.compute()
Print the variogram contents
varioexp
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
gp.setDefault(dims=[6,6])
gp.varmod(varioexp)
<AxesSubplot:>
Represent the variogram with the number of pairs. By storing the Axes in a variable ax returned by the last instruction, we prevent jupyter from displaying '\AxesSubplot:\' before the plot.
ax = gp.varmod(varioexp,show_pairs=True)
Automatic Fitting procedure (with default fitting parameters).
fitmod = gl.Model()
err = fitmod.fit(varioexp)
gp.varmod(varioexp, fitmod)
plt.show() # This last instruction also prevents jupyter from displaying a text before the plot
Print the Fittied Model content
fitmod
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 Covariances structures known by gstlearn
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 (using given basic structures)
types = [gl.ECov.NUGGET, gl.ECov.CUBIC, gl.ECov.SPHERICAL]
err = fitmod.fit(varioexp, types=types)
ax = gp.varmod(varioexp, fitmod)
The resulting Model
fitmod
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
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
fitmod
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
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
fitmod
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
varioParamMulti = gl.VarioParam.createMultiple(4, 15, 15.)
vario_4dir = gl.Vario(varioParamMulti, dat)
err = vario_4dir.compute()
ax = gp.varmod(vario_4dir, flagLegend=True)
Fitting a Multi-directional variogram
model_4dir = gl.Model()
err = model_4dir.fit(vario_4dir,types=types)
ax = gp.varmod(vario_4dir, model_4dir)
Calculating Variogram Map (the result is returned as a Grid data base displayed using raster function). At the right is displayed the number of pairs used for the variogram calculation in each grid cell.
grid_vmap = gl.db_vmap_compute(dat, gl.ECalcVario.VARIOGRAM)
grid_vmap.display()
fig, ax = plt.subplots(1,2,figsize=[14,10])
fig.tight_layout(pad=5.0)
ax[0].raster(grid_vmap, flagLegend=True)
ax[1].raster(grid_vmap, name="*Nb", flagLegend=True)
plt.show()
Data Base Grid Characteristics ============================== Data Base Summary ----------------- File is organized as a regular grid Space dimension = 2 Number of Columns = 5 Maximum Number of UIDs = 5 Total number of samples = 1681 Grid characteristics: --------------------- Origin : -382.500 -678.500 Mesh : 19.125 33.925 Number : 41 41 Variables --------- Column = 0 - Name = rank - Locator = NA Column = 1 - Name = x1 - Locator = x1 Column = 2 - Name = x2 - Locator = x2 Column = 3 - Name = VMAP.January_temp.Var - Locator = z1 Column = 4 - Name = VMAP.January_temp.Nb - Locator = NA
Automatic Model Fitting from Variogram Map
modelVM = gl.Model()
err = modelVM.fitFromVMap(grid_vmap, types=types)
modelVM
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
err = gl.dbgrid_model(grid_vmap, modelVM)
grid_vmap.display()
ax = gp.raster(grid_vmap, flagLegend=True)
Data Base Grid Characteristics ============================== Data Base Summary ----------------- File is organized as a regular grid Space dimension = 2 Number of Columns = 6 Maximum Number of UIDs = 6 Total number of samples = 1681 Grid characteristics: --------------------- Origin : -382.500 -678.500 Mesh : 19.125 33.925 Number : 41 41 Variables --------- Column = 0 - Name = rank - Locator = NA Column = 1 - Name = x1 - Locator = x1 Column = 2 - Name = x2 - Locator = x2 Column = 3 - Name = VMAP.January_temp.Var - Locator = NA Column = 4 - Name = VMAP.January_temp.Nb - Locator = NA Column = 5 - Name = VMAP.Model - Locator = z1
Compare Directional Variograms and fitted model from Variogram Map
ax = gp.varmod(vario_4dir, modelVM)
