Tutorial on SPDE
gstlearn Team
Preamble
In this preamble, we load the gstlearn library, clean the workspace.
rm(list=ls())
library(gstlearn)
library(ggplot2)
library(ggpubr)
library(ggnewscale)Getting the Data Bases from the repository.
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.
fileNF = loadData("Scotland", "Scotland_Temperatures.NF")
temperatures = Db_createFromNF(fileNF)
temperatures$setLocator("Elevation", ELoc_F())## NULLiuid = temperatures$addSelection(as.numeric(! is.na(temperatures["J*"])))
temperatures$display()##
## 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 = NewSel - Locator = sel## NULLThe 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).
fileNF = loadData("Scotland", "Scotland_Elevations.NF")
grid = DbGrid_createFromNF(fileNF)
grid$display()##
## 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## NULLVariography
Calculate the omni-directional variogram of the temperature for 18 lags of 25.
vparam = VarioParam_createOmniDirection(npas=18, dpas=25)
vario = Vario(vparam)
err = vario$compute(temperatures)
vario$display()##
## Variogram characteristics
## =========================
## Number of variable(s) = 1
## Number of direction(s) = 1
## Space dimension = 2
## Variable(s) = [January_temp]
##
## Variance-Covariance Matrix 1.020
##
## Direction #1
## ------------
## Number of lags = 18
## Direction coefficients = 1.000 0.000
## Direction angles (degrees) = 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## NULLp = ggDefault()
p = p + plot.varmod(vario)
p = p + plot.decoration(title="Variogram of Temperature (Raw)")
ggPrint(p)
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.
vparam = VarioParam_createOmniDirection(npas=18, dpas=25)
vario = Vario(vparam)
md = Model()
md$setDriftIRF(order=0, nfex=1)## NULLvario$compute(temperatures,model=md)## [1] 0vario$display()##
## Variogram characteristics
## =========================
## Number of variable(s) = 1
## Number of direction(s) = 1
## Space dimension = 2
## Variable(s) = [January_temp]
##
## Variance-Covariance Matrix 0.363
##
## Direction #1
## ------------
## Number of lags = 18
## Direction coefficients = 1.000 0.000
## Direction angles (degrees) = 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## NULLp = ggDefault()
p = p + plot.varmod(vario)
p = p + plot.decoration(title="Variogram of Temperature (Residuals)")
ggPrint(p)
Fit the variogram of residuals in a model having drifts. Some constraints have been added during the fitting step.
model = md
structs = ECov_fromKeys(c("NUGGET", "MATERN"))
const = Constraints()
err = const$addItemFromParamId(EConsElem_SILL(),icov=0,type=EConsType_UPPER(),value=0.1)
err = const$addItemFromParamId(EConsElem_PARAM(),icov=1, type = EConsType_EQUAL(),value = 1)
err = const$addItemFromParamId(EConsElem_RANGE(),icov=1, type = EConsType_LOWER(),value = 100)
err = const$addItemFromParamId(EConsElem_RANGE(),icov=1, type = EConsType_UPPER(),value = 350)
err = md$fit(vario,structs,constraints=const)p = ggDefault()
p = p + plot.varmod(vario, model)
p = p + plot.decoration(title="Variogram of Temperature (Residuals)")
ggPrint(p)
Identify the Drift
Derive the parameter of a Global Trend (composed of the Elevation as a drift function) using the SPDE formalism.
spde = SPDE(model,grid,temperatures,ESPDECalcMode_KRIGING())
coeffs = spde$getCoeffs()
cat("Trend coefficients:", coeffs)## Trend coefficients: 3.970026 -0.00659772Represent the scatter plot of the Temperature given the Elevation and add the Global Trend (caledlated beforehand)
p = ggplot()
p = p + plot.correlation(temperatures,namex="Elevation",namey="*temp", asPoint=TRUE)
p = p + geom_segment(aes(x = 0, y = coeffs[1], xend = 400, yend = coeffs[1] + 400 * coeffs[2]))
ggPrint(p)
Estimation
We perform the Estimation in the SPDE framework (considering the Elevation as the Global Trend)
err = krigingSPDE(temperatures, grid, model)p = ggDefaultGeographic()
p = p + plot.grid(grid)
p = p + plot.decoration(title="Temperature (using Elevation as global Trend)")
ggPrint(p)
We also perform the Estimation by Kriging (using Elevation as External Drift). This estimation is performed in Unique Neighborhood.
neighU = NeighUnique_create()
err = kriging(temperatures, grid, model, neighU)p = ggDefaultGeographic()
p = p + plot.grid(grid)
p = p + plot.decoration(title="Temperature (with Elevation as External Drift)")
ggPrint(p)
Comparing both estimates
p = ggplot()
p = p + plot.correlation(grid,namex="Kriging.January_temp.estim",
namey="KrigingSPDE.January_temp.estim", asPoint=TRUE)
p = p + plot.decoration(xlab="Kriging with External Drift", ylab="SPDE Kriging")
ggPrint(p)