Variography
gstlearn Team
In this preamble, we load the gstlearn library and clean the workspace.
Then the necessary data set is downloaded and named dat: the target variable is January_temp
dlfile = "https://soft.minesparis.psl.eu/gstlearn/data/Scotland/Scotland_Temperatures.NF"
fileNF = "Scotland_Temperatures.NF"
download.file(dlfile, fileNF)
dat = Db_createFromNF(fileNF)\frametitle{Experimental Variograms} \begin{itemize} \item Data is a \alert{regionalized} variable $$z_i = z(x_i)$$ \item 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$ \end{itemize}
\frametitle {Variogram Cloud}
varioParamOmni = VarioParam_createOmniDirection(npas = 100)
grid.cloud = db_variogram_cloud(dat, varioParamOmni)
grid.cloud$display()##
## 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## NULLp = ggplot()
p = p + plot.grid(grid.cloud, "Cloud.January*")
p = p + plot.geometry(asp=0)
ggPrint(p)
\frametitle {Experimental Variogram}
We calculate the omni-directional variogram of the temperatures.
varioParamOmni = VarioParam_createOmniDirection(npas=40, dpas=10)
varioexp = Vario(varioParamOmni, dat)
err = varioexp$compute()Print the variogram contents
varioexp\frametitle {Experimental Variogram}
Plot the omni-directional variogram
ggplot() + plot.varmod(varioexp)
\frametitle {Experimental Variogram}
ggplot() + plot.varmod(varioexp,draw_psize=TRUE)
\frametitle {Automatic Model Fitting}
fitmod = Model()
err = fitmod$fit(varioexp)
ggplot() + plot.varmod(varioexp, fitmod)
\frametitle {Automatic Model Fitting}
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\frametitle {List of Basic structures}
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## NULL\frametitle {Automatic Fitting (with given basic structures)}
types = ECov_fromKeys(c("NUGGET","CUBIC","SPHERICAL"))
err = fitmod$fit(varioexp, types=types)
ggplot() + plot.varmod(varioexp, fitmod)
\frametitle {Automatic Fitting (with given basic structures)}
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\frametitle{Model Fitting with Inequality constraints}
constraints = Constraints()
err = constraints$addItemFromParamId(EConsElem_RANGE(),icov=1,type=EConsType_UPPER(),value=20.)
err = constraints$addItemFromParamId(EConsElem_SILL(),icov=1,type=EConsType_LOWER(),value=0.03)
err = fitmod$fit(varioexp, types=types, constraints, Option_VarioFit(TRUE))
ggplot() + plot.varmod(varioexp, fitmod)
\frametitle{Model Fitting with Inequality constraints}
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\frametitle{Model Fitting with Equality constraints}
constraints = Constraints()
err = constraints$addItemFromParamId(EConsElem_RANGE(),icov=1,type=EConsType_EQUAL(),value=1000.)
err = constraints$addItemFromParamId(EConsElem_SILL(),icov=1,type=EConsType_EQUAL(),value=0.4)
err = fitmod$fit(varioexp, types=types, constraints, Option_VarioFit(flag_noreduce=TRUE))
ggplot() + plot.varmod(varioexp, fitmod)
\frametitle{Model Fitting with Equality constraints}
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\frametitle{Directional Variograms}
varioParamMulti = VarioParam_createMultiple(ndir=4, npas=15, dpas=15.)
vario.4dir = Vario(varioParamMulti, dat)
err = vario.4dir$compute()
ggplot() + plot.varmod(vario.4dir)
\frametitle{Directional Variograms}
model.4dir = Model()
err = model.4dir$fit(vario.4dir,types=types)
ggplot() + plot.varmod(vario.4dir, model.4dir)
\frametitle{Calculating Variogram Map}
grid.vmap = db_vmap_compute(dat, ECalcVario_VARIOGRAM())
ggplot() + plot.grid(grid.vmap)
\frametitle{Automatic Model Fitting from Variogram Map}
modelVM = 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\frametitle{Automatic Model Fitting from Variogram Map}
err = dbgrid_model(grid.vmap, modelVM)
ggplot() + plot.grid(grid.vmap)
\frametitle{Compare Directional Variograms and Variogram Map}
ggplot() + plot.varmod(vario.4dir, modelVM)