Kriging linear combinations
gstlearn core team
January 31, 2024
Introduction
This notebook presents examples for the co-kriging of linear combinations. Starting from a set of variables \(\mathbf{Z}(\alpha) = (Z_i(\alpha))_{i \in 1:n}\), the task is to compute the co-kriging of a linear combination of the variable at location \(\mathbf{o}\): \[ (\beta^T \mathbf{Z}(\mathbf{o}))^{K} = \sum_{i \in 1:n} \sum_{\alpha \in A_i} \lambda_i^{\alpha} Z_i(\alpha) \] Kriging is linear, hence \((\beta^T \cdot \mathbf{Y}(\mathbf{o}))^{K} = \sum_{i \in 1:p} \beta_i \times Y_i(\mathbf{o})^{K}\) but the variance of the kriging error is not easily computable. However, the matLC option of the kriging function offers the possibility to compute directly the estimate and its kriging variance of any linear combinations. The method is evaluated on a simulated data set.
The following cases are considered:
- Ordinary cokriging using a IRF-0 model (Intrinsic Random Function of order 0 with unknown means)
- Simple cokriging using a SOS model (Second Order Stationary random function with known means)
- Homotopic data set (both variables are defined for each location)
- Heterotopic data set (variables may be defined or not at each location)
- Direct kriging of the linear combination or linear combination of the estimates
In all cases, the results seem consistent with the theory.
Kriging systems when kriging a linear combination are given below:
Simple kriging with a second order stationary model (SK/SOS)
The expression of the estimator is: \[ (\beta^T \, \mathbf{Z_{(o)}})^{SK} = \beta^T \mathbf{m} +\sum_{i \in 1:n} \sum_{\alpha \in A_i} \lambda^{\alpha} (Z_i(\alpha) - m_i) \]
There is no bias as: \[ \mathbb{E}\{(\beta^T \, \mathbf{Z_{(o)}})^{SK}\} = \beta^T \mathbf{m} = \mathbb{E}\{\beta^T \, \mathbf{Z_{(o)}}\}. \]
The variance of the kriging error is: \[ \begin{array}{cl} \text{Var}\{\beta^T \, \mathbf{Z_{(o)}} - (\beta^T \, \mathbf{Z_{(o)}})^{SK}\} = & \beta^T \, \text{Var}\{\mathbf{Z_{(o)}}\} \, \beta \\ & - 2 \sum_{k,i \in 1:n} \sum_{\alpha \in A_i} \beta_k \, \lambda^{\alpha} \text{Cov}(Z_k(o),Z_i(\alpha)) \\ & + \sum_{i,j \in 1:n} \sum_{\alpha \in A_i, \alpha' \in A_j} \beta_k \, \lambda^{\alpha} \lambda^{\alpha'} \text{Cov}(Z_i(\alpha),Z_j(\alpha')) \end{array} \]
Minimizing this variance, we gets the kriging system of for the estimation of the linear combination: \[ \left( \text{Cov}(Z_{i}(\alpha),Z_{j}(\alpha')) \right) \Lambda = \sum_{k\in 1:n} \beta_k \text{Cov}(Z_k(o), Z_i(\alpha)) \]
The variance of the kriging error boils down to: \[ \text{Var}\{\beta^T \, \mathbf{Z_{(o)}} - (\beta^T \, \mathbf{Z_{(o)}})^{SK}\} = \beta^T \, \text{Var}\{\mathbf{Z_{(o)}}\} \, \beta - \sum_{k,i \in 1:n} \sum_{\alpha \in A_i} \beta_k \text{Cov}(Z_k(o),Z_i(\alpha)) \lambda^{\alpha} \]
ordinary kriging with a Intrinsic Random Fonction of order 0 (OK/IRF0)
The expression of the estimator is: \[ (\beta^T \, \mathbf{Z_{(o)}})^{OK} = \sum_{i \in 1:n} \sum_{\alpha \in A_i} \lambda^{\alpha} Z_i(\alpha) \]
The no-bias conditions are: \[ \sum_{\alpha \in A_i} \lambda^{\alpha} = \beta_i, \text{ for } \forall i \in 1:n. \]
The variance of the kriging error is: \[ \begin{array}{cl} \text{Var}\{\beta^T \, \mathbf{Z_{(o)}} - (\beta^T \, \mathbf{Z_{(o)}})^{OK}\} = & \beta^T \, \text{Var}\{\mathbf{Z_{(o)}}\} \, \beta \\ & - 2 \sum_{k\in 1:n} \sum_{i, \alpha \in A_i} \beta_k \text{Cov}(Z_k(o),Z_i(\alpha)) \, \lambda^{\alpha} \\ & + \sum_{i, \alpha \in A_i} \sum_{j, \alpha' \in A_j} \lambda^{\alpha} \text{Cov}(Z_i(\alpha),Z_j(\alpha')) \lambda^{\alpha'} \end{array} \]
Minimizing the Lagrangian, \[ \mathcal{L}(\Lambda, \mu) = \text{Var}\{\beta^T \, \mathbf{Z_{(o)}} - (\beta^T \, \mathbf{Z_{(o)}})^{OK}\} + \sum_{k \in 1:n} \mu_k (\sum_{\alpha \in A_k} \lambda^{\alpha} - \beta_k) \]
we gets the kriging system of for the estimation of the linear combination: \[ \left[ \begin{array}{cc} \text{Cov}(Z_i(\alpha), Z_j(\alpha')) & \mathbf{1}_{\alpha \in A_j} \\ \mathbf{1}_{\alpha' \in A_i} & \mathbf{0} \end{array} \right] \left[ \begin{array}{c} \lambda \\ \mu \end{array} \right] = \left[ \begin{array}{c} \sum_{k \in 1:n} \beta_k\text{Cov}(Z_k(o), Z_i(\alpha)) \\ \beta_i \end{array} \right] \]
The variance of the kriging error boils down to: \[ \text{Var}\{\beta^T \, \mathbf{Z_{(o)}} - (\beta^T \, \mathbf{Z_{(o)}})^{OK}\} = \beta^T \, \text{Var}\{\mathbf{Z_{(o)}}\} \, \beta - \sum_{k \in 1:n} \sum_{i,\alpha \in A_i} \beta_k \text{Cov}(Z_k(o),Z_i(\alpha))\, \lambda^{\alpha} - \sum_{k \in 1:n} \beta_k \, \mu_k \]
Initialisation
Simulation of a reference data set
We create a reference data set using the the Turning Bands method (simtub):
the distribution is log normal,
the stationary covariance function is exponential.
5% of the two dimensional grid is sampled and considered as the reference data set.
# parameters for the simulation
m = c(1, 2)
sig = c(0.25, 0.75)
rho = 0.7
CO = matrix(c(1, rho, rho, 1), 2, 2)
# initialization of the grid
grd = DbGrid_create(x0=c(0.0,0.0), dx=c(0.01,0.01), nx=c(100,100))
# simulation from a model
model = Model_createFromParam(type = ECov_EXPONENTIAL(), range=0.2, sills = CO)
err = simtub(dbin = NULL, dbout = grd, model = model, nbsimu = 1)
err = grd$setName("Simu.V1", "Y1")
err = grd$setName("Simu.V2", "Y2")
grd["Z1"] = m[1] * exp(sig[1] * grd["Y1"] - sig[1]**2 / 2)
grd["Z2"] = m[2] * exp(sig[2] * grd["Y2"] - sig[2]**2 / 2)
# Data set (ssp proportion of the grid)
data = Db_createFillRandom(ndat = ssp*grd$getNSample(), nvar = 0)
err = data$setName("x-1", "x1")
err = data$setName("x-2", "x2")
err = migrate(dbin = grd, name = "Z1", dbout = data, dmax = c(0.1, 0.1),
namconv = NamingConvention(""))
err = migrate(dbin = grd, name = "Z2", dbout = data, dmax = c(0.1, 0.1),
namconv = NamingConvention(""))
# plots
p1 = plot.init(asp=1) +
plot.raster(grd, name = "Z1", palette = "Spectral") +
plot.symbol(data, colour = "red") +
plot.decoration(title="Z1")
p2 = plot.init() +
plot.hist(grd, name='Z1', bins = 25, bg = "orange", color= "gray") +
plot.decoration(xlab = "Raw variable", title="Histogram of the initial variable Z1")
p3 = plot.init(asp=1) +
plot.raster(grd, name = "Z2", palette = "Spectral") +
plot.symbol(data, colour = "red") +
plot.decoration(title="Z2")
p4 = plot.init() +
plot.hist(grd, name='Z2', bins = 25, bg = "orange", color= "gray") +
plot.decoration(xlab = "Raw variable", title="Histogram of the initial variable Z2")
ggarrange(p1, p2, p3, p4, ncol = 2, nrow = 2)
# defining the linear combination
beta = c(1, -3) # coefficients and computing the linear combination
nm_var = c("Z1", "Z2")
nm_LC = c("LC", "Z1")
matLC = matrix(c(beta, 1, 0), nrow = 2, ncol = 2, byrow = TRUE)
colnames(matLC) <- nm_var
rownames(matLC) <- nm_LC
# computing the linear computation
data["LC"] = matrix(data$getColumns(colnames(matLC)), nrow = data$getNSample(),
ncol = dim(matLC)[2]) %*% t(matLC)[,1]
grd["LC"] = matrix(grd$getColumns(colnames(matLC)), nrow = grd$getNSample(),
ncol = dim(matLC)[2]) %*% t(matLC)[,1]
# statistics
tab = rbind(
dbStatisticsMono(grd, names = c("Z*", "LC"), opers = opers)$toTL(),
dbStatisticsMono(data, names = c("Z*", "LC"), opers = opers)$toTL())
rownames(tab) <- c(
"grid-Z1", "grid-Z2", "grid-LC",
"samples-Z1", "samples-Z2", "samples-LC")
knitr::kable(tab, caption = "Statistics on the reference data set", digits = 4)| Number | Minimum | Maximum | Mean | St. Dev. | |
|---|---|---|---|---|---|
| grid-Z1 | 10000 | 0.4630 | 2.0321 | 0.9811 | 0.2220 |
| grid-Z2 | 10000 | 0.1315 | 20.7945 | 2.1241 | 1.6525 |
| grid-LC | 10000 | -60.5413 | 0.2857 | -5.3913 | 4.8201 |
| samples-Z1 | 50 | 0.6299 | 1.5488 | 1.0274 | 0.2347 |
| samples-Z2 | 50 | 0.2570 | 10.3052 | 2.3650 | 1.9771 |
| samples-LC | 50 | -29.4314 | -0.1410 | -6.0675 | 5.7941 |
Global parameters for spatial modelling
# Variogram parameters (omnidirectional)
varioParam = VarioParam_createOmniDirection(nlag=10, dlag=0.05)
# Variogram basic structures
types = ECov_fromKeys(c("NUGGET", "EXPONENTIAL", "EXPONENTIAL"))
# types = ECov_fromKeys(c("EXPONENTIAL"))
# weighted proportional to the number of pairs and inverse proportional to the distance
opt = Option_AutoFit()
err = opt$setWmode(2)
#define neighborhood
neigh = NeighMoving_create(nmaxi=30, radius=0.5)
# Display the kriging system for a fixed target when using kriging
OptDbg_setReference(0)## NULL# function to compute the model of the linear combination (using a square matrix)
# TODO: generalize to rectangular matrix?
model_of_LC <- function(model, A) {
stopifnot(dim(A)[1] == dim(A)[2])
nv = model$getNVar()
stopifnot(nv == dim(A)[1])
ns = model$getNCov()
model_LC = model$clone()
for (k in 1:ns) {
CC = A %*% model$getSills(k-1)$toTL() %*% t(A)
for (i in 1:nv) {
for (j in 1:nv) {
err = model_LC$setSill(icov = k-1, ivar=i-1, jvar=j-1, value = CC[i, j])
}
}
}
model_LC
}Ordinary kriging (OK)
Intrinsic Random Function of order 0
# Locate Z
err = data$setLocators("Z*", locatorType = ELoc_Z(), cleanSameLocator = TRUE)
# Compute variogram
var_Z = Vario_computeFromDb(varioParam, data)
# fit model
mod_Z_IRF0 = Model()
err = mod_Z_IRF0$fit(var_Z, types = types, mauto = opt)
# setting the drift "unknown mean" for Ordinary Kriging
err = mod_Z_IRF0$setDriftIRF(order = 0)
# plot
multi.varmod(var_Z, mod_Z_IRF0, flagLegend=TRUE)
if(flag.verbose) {
mod_Z_IRF0
}Ordinary kriging of the raw variables
err = grd$deleteColumns(paste("OK", nm_var, "*", sep = "."))
err = grd$deleteColumns(paste("LC", "OK", "*", sep = "."))
if(flag.verbose) {
print(paste0(">>> Ordinary kriging of the raw variables"))
}
err = data$setLocators(c("Z1", "Z2"), ELoc_Z(), cleanSameLocator=TRUE)
err = kriging(data, grd, mod_Z_IRF0, neigh,
flag_est = TRUE, flag_std = TRUE,
namconv = NamingConvention("OK"))
# computing the LC of kriging
grd["LC.OK.estim"] = matrix(grd$getColumns(names = paste("OK", colnames(matLC), "estim", sep = ".")),
nrow = grd$getNSample(), ncol = dim(matLC)[2]) %*% t(matLC)[,1]
# statistics
vn = c("Z1", paste("OK", "Z1", c("estim", "stdev"), sep = "."),
"Z2", paste("OK", "Z2", c("estim", "stdev"), sep = "."))
tab = dbStatisticsMono(grd, vn, opers = opers)$toTL()
knitr::kable(tab, caption = "Statistics on the grid for Ordinary Kriging with the initial model", digits = 4)| Number | Minimum | Maximum | Mean | St. Dev. | |
|---|---|---|---|---|---|
| Z1 | 10000 | 0.4630 | 2.0321 | 0.9811 | 0.2220 |
| OK.Z1.estim | 10000 | 0.7421 | 1.5329 | 1.0235 | 0.0986 |
| OK.Z1.stdev | 10000 | 0.1327 | 0.2430 | 0.2112 | 0.0168 |
| Z2 | 10000 | 0.1315 | 20.7945 | 2.1241 | 1.6525 |
| OK.Z2.estim | 10000 | 1.6435 | 3.2778 | 2.3509 | 0.3262 |
| OK.Z2.stdev | 10000 | 2.0403 | 2.1457 | 2.0527 | 0.0142 |
# plot of estimated values
p1 = plot.init(asp=1) +
plot.raster(grd, name = paste("OK", "Z1", "estim", sep = "."),
palette = "Spectral", flagLegend = TRUE, legendName = "estim") +
plot.decoration(xlab = "", ylab = "", title = paste0("Ordinary coK for Z1"))
p2 = plot.init(asp=1) +
plot.raster(grd, name = paste("OK", "Z1", "stdev", sep = "."),
palette = "Spectral", flagLegend = TRUE, legendName = "stdev") +
plot.decoration(xlab = "", ylab = "", title = paste0("Std. of Ordinary coK for Z1"))
p3 = plot.init(asp=1) +
plot.raster(grd, name = paste("OK", "Z2", "estim", sep = "."),
palette = "Spectral", flagLegend = TRUE, legendName = "estim") +
plot.decoration(xlab = "", ylab = "", title = paste0("Ordinary coK for Z2"))
p4 = plot.init(asp=1) +
plot.raster(grd, name = paste("OK", "Z2", "stdev", sep = "."),
palette = "Spectral", flagLegend = TRUE, legendName = "stdev") +
plot.decoration(xlab = "", ylab = "", title = paste0("Std. of Ordinary coK for Z2"))
ggarrange(p1, p2, p3, p4, ncol = 2, nrow = 2)
# plot Z vs. Z*
p1 = plot.init() +
plot.correlation(grd, namex = paste("OK", "Z1", "estim", sep = "."), namey = "Z1", asPoint = TRUE,
flagDiag = TRUE, diagColor = "red", flagSameAxes = TRUE) +
plot.decoration(xlab = "Estimated values", ylab = "Actual values", title = "Ordinary coK of Z1")
p2 = plot.init() +
plot.correlation(grd, namex = paste("OK", "Z2", "estim", sep = "."), namey = "Z2", asPoint = TRUE,
flagDiag = TRUE, diagColor = "red", flagSameAxes = TRUE) +
plot.decoration(xlab = "Estimated values", ylab = "Actual values", title = "Ordinary coK of Z2")
ggarrange(p1, p2, ncol = 2, nrow = 1)
Ordinary kriging of a linear combination of the raw variables
err = grd$deleteColumns("OK-LC*")
if(flag.verbose) {
print(paste0(">>> Ordinary kriging of a linear combination of the variables"))
}
mod_LC_IRF0 = model_of_LC(mod_Z_IRF0, A = matLC)
# modelling the linear combination
# Locate Z
err = data$setLocators(nm_LC, locatorType = ELoc_Z(), cleanSameLocator = TRUE)
# computing the model of the LC from the initial model
var_LC = Vario_computeFromDb(varioParam, data)
# setting the drift "unknown mean" for Ordinary Kriging
err = mod_LC_IRF0$setDriftIRF(order = 0)
# plot
multi.varmod(var_LC, mod_LC_IRF0, flagLegend=TRUE)
if(flag.verbose) {
mod_LC_IRF0
}
# kriging the LC variables
err = data$setLocators(nm_LC, ELoc_Z(), cleanSameLocator=TRUE)
err = kriging(data, grd, mod_LC_IRF0, neigh,
flag_est = TRUE, flag_std = TRUE,
namconv = NamingConvention("OK-LC"))
# statistics
vn = c(
paste(c("OK.Z1", "OK-LC.Z1", "OK-LC.LC", "LC.OK"), "estim", sep = "."),
paste(c("OK.Z1", "OK-LC.Z1", "OK-LC.LC"), "stdev", sep = ".")
)
tab = dbStatisticsMono(grd, names = vn, opers = opers)$toTL()
knitr::kable(tab, caption = "Statistics on the grid for Ordinary Kriging with the LC model",
digits = 4)| Number | Minimum | Maximum | Mean | St. Dev. | |
|---|---|---|---|---|---|
| OK.Z1.estim | 10000 | 0.7421 | 1.5329 | 1.0235 | 0.0986 |
| OK-LC.Z1.estim | 10000 | 0.7421 | 1.5329 | 1.0235 | 0.0986 |
| OK-LC.LC.estim | 10000 | -8.5591 | -3.7712 | -6.0292 | 0.9327 |
| LC.OK.estim | 10000 | -8.5591 | -3.7712 | -6.0292 | 0.9327 |
| OK.Z1.stdev | 10000 | 0.1327 | 0.2430 | 0.2112 | 0.0168 |
| OK-LC.Z1.stdev | 10000 | 0.1327 | 0.2430 | 0.2112 | 0.0168 |
| OK-LC.LC.stdev | 10000 | 5.9915 | 6.3041 | 6.0301 | 0.0417 |
# plot Z vs. Z*
p1 = plot.init() +
plot.correlation(grd,
namex = paste("OK", "Z1", "estim", sep = "."),
namey = paste("OK-LC", "Z1", "estim", sep = "."),
asPoint = TRUE, flagDiag = TRUE, diagColor = "red", flagSameAxes = TRUE) +
plot.decoration(xlab = "Initial estimate", ylab = "Estimate with the LC model",
title = "Ordinary coK of Z1")
p2 = plot.init() +
plot.correlation(grd,
namey = paste("OK-LC", "LC", "estim", sep = "."),
namex = paste("LC", "OK", "estim", sep = "."),
asPoint = TRUE, flagDiag = TRUE, diagColor = "red", flagSameAxes = TRUE) +
plot.decoration(xlab = "LC of the estimates", ylab = "Estimate with the LC model",
title = "Ordinary coK of CL")
ggarrange(p1, p2, ncol = 2, nrow = 1)
Direct computation of linear combination kriging
err = grd$deleteColumns("dOK*")
if(flag.verbose) {
print(paste0(">>> Direct computation of linear combination kriging"))
}
# the R matrix matLC has to be converted into a gstlearn MatrixDense:
# MatrixDense_createFromVD(beta, nrow=1, ncol=length(beta))
# OR
# beta |> matrix(nrow=1, ncol=length(beta)) |> fromTL()
# OR
# fromTL(matrix(beta, nrow=1, ncol=length(beta)))
# kriging LC variables
err = data$setLocators(nm_var, ELoc_Z(), cleanSameLocator = TRUE)
krigopt = KrigOpt()
krigopt$setMatLC(fromTL(matrix(beta, nrow = 1, ncol = length(beta))))## [1] 0err = kriging(data, grd, mod_Z_IRF0, neigh,
flag_est = TRUE, flag_std = TRUE, krigopt=krigopt,
namconv = NamingConvention("dOK"))
# statistics
vn = c(
paste(c("OK-LC.LC", "dOK.LC", "LC.OK"), "estim", sep = "."),
paste(c("OK-LC.LC", "dOK.LC"), "stdev", sep = ".")
)
tab = dbStatisticsMono(grd, names = vn, opers = opers, flagIso = FALSE)$toTL()
knitr::kable(tab, caption = "Statistics on the grid for Ordinary Kriging of LC", digits = 4)| Number | Minimum | Maximum | Mean | St. Dev. | |
|---|---|---|---|---|---|
| OK-LC.LC.estim | 10000 | -8.5591 | -3.7712 | -6.0292 | 0.9327 |
| dOK.LC.estim | 10000 | -8.5591 | -3.7712 | -6.0292 | 0.9327 |
| LC.OK.estim | 10000 | -8.5591 | -3.7712 | -6.0292 | 0.9327 |
| OK-LC.LC.stdev | 10000 | 5.9915 | 6.3041 | 6.0301 | 0.0417 |
| dOK.LC.stdev | 10000 | 5.9915 | 6.3041 | 6.0301 | 0.0417 |
# plot direct LC estimate vs. estimate using the LC model
p1 = plot.init() +
plot.correlation(grd,
namex = paste("OK-LC", "LC", "estim", sep = "."),
namey = paste("dOK", "LC", "estim", sep = "."),
asPoint = TRUE, flagDiag = TRUE, diagColor = "red", flagSameAxes = TRUE) +
plot.decoration(xlab = "Estimate using LC model", ylab = "Direct estimate of the LC",
title = "Ordinary coK: estimated value")
p2 = plot.init() +
plot.correlation(grd,
namex = paste("OK-LC", "LC", "stdev", sep = "."),
namey = paste("dOK", "LC", "stdev", sep = "."),
asPoint = TRUE, flagDiag = TRUE, diagColor = "red", flagSameAxes = TRUE) +
plot.decoration(xlab = "Estimate using LC model", ylab = "Direct estimate of the LC",
title = "Ordinary coK: Std. of kriging error")
ggarrange(p1, p2, ncol = 2, nrow = 1)
Simple kriging (SK)
Second Order Stationary model of the raw variable Z
# copying the IRF0 model into SOS model
mod_Z_SOS = mod_Z_IRF0$clone()
err = mod_Z_SOS$setDriftIRF(-1)
err = mod_Z_SOS$setMeans(m)
if(flag.verbose) {
err = mod_Z_SOS$display()
}
# defining the model of the LC
mod_LC_SOS = mod_LC_IRF0$clone()
err = mod_LC_SOS$setDriftIRF(-1)
err = mod_LC_SOS$setMeans(as.numeric(matLC %*% m))
if(flag.verbose) {
err = mod_LC_SOS$display()
}Simple kriging of the raw variables
err = grd$deleteColumns(paste("SK", nm_var, "*", sep = "."))
err = grd$deleteColumns(paste("LC", "SK", "*", sep = "."))
if(flag.verbose) {
print(paste0(">>> Simple kriging of the raw variables"))
}
err = data$setLocators(c("Z1", "Z2"), ELoc_Z(), cleanSameLocator=TRUE)
err = kriging(data, grd, mod_Z_SOS, neigh,
flag_est = TRUE, flag_std = TRUE,
namconv = NamingConvention("SK"))
# computing the LC of kriging
grd["LC.SK.estim"] = matrix(grd$getColumns(names = paste("SK", colnames(matLC), "estim", sep = ".")),
nrow = grd$getNSample(), ncol = dim(matLC)[2]) %*% t(matLC)[,1]
# statistics
vn = c("Z1", paste("SK", "Z1", c("estim", "stdev"), sep = "."),
"Z2", paste("SK", "Z2", c("estim", "stdev"), sep = "."))
tab = dbStatisticsMono(grd, vn, opers = opers)$toTL()
knitr::kable(tab, caption = "Statistics on the grid for Simple Kriging with the initial model",
digits = 4)| Number | Minimum | Maximum | Mean | St. Dev. | |
|---|---|---|---|---|---|
| Z1 | 10000 | 0.4630 | 2.0321 | 0.9811 | 0.2220 |
| SK.Z1.estim | 10000 | 0.7416 | 1.5086 | 0.9972 | 0.0824 |
| SK.Z1.stdev | 10000 | 0.1301 | 0.2240 | 0.2078 | 0.0159 |
| Z2 | 10000 | 0.1315 | 20.7945 | 2.1241 | 1.6525 |
| SK.Z2.estim | 10000 | 1.9990 | 2.0019 | 2.0000 | 0.0003 |
| SK.Z2.stdev | 10000 | 2.0071 | 2.0071 | 2.0071 | 0.0000 |
# plot of estimated values
p1 = plot.init(asp=1) +
plot.raster(grd, name = paste("SK", "Z1", "estim", sep = "."),
palette = "Spectral", flagLegend = TRUE, legendName = "estim") +
plot.decoration(xlab = "", ylab = "", title = paste0("Simple coK for Z1"))
p2 = plot.init(asp=1) +
plot.raster(grd, name = paste("SK", "Z1", "stdev", sep = "."),
palette = "Spectral", flagLegend = TRUE, legendName = "stdev") +
plot.decoration(xlab = "", ylab = "", title = paste0("Std. of Simple coK for Z1"))
p3 = plot.init(asp=1) +
plot.raster(grd, name = paste("SK", "Z2", "estim", sep = "."),
palette = "Spectral", flagLegend = TRUE, legendName = "estim") +
plot.decoration(xlab = "", ylab = "", title = paste0("Simple coK for Z2"))
p4 = plot.init(asp=1) +
plot.raster(grd, name = paste("SK", "Z2", "stdev", sep = "."),
palette = "Spectral", flagLegend = TRUE, legendName = "stdev") +
plot.decoration(xlab = "", ylab = "", title = paste0("Std. of Simple coK for Z2"))
ggarrange(p1, p2, p3, p4, ncol = 2, nrow = 2)
# plot Z vs. Z*
p1 = plot.init() +
plot.correlation(grd, namex = paste("SK", "Z1", "estim", sep = "."), namey = "Z1", asPoint = TRUE,
flagDiag = TRUE, diagColor = "red", flagSameAxes = TRUE) +
plot.decoration(xlab = "Estimated values", ylab = "Actual values", title = "Simple coK of Z1")
p2 = plot.init() +
plot.correlation(grd, namex = paste("SK", "Z2", "estim", sep = "."), namey = "Z2", asPoint = TRUE,
flagDiag = TRUE, diagColor = "red", flagSameAxes = TRUE) +
plot.decoration(xlab = "Estimated values", ylab = "Actual values", title = "Simple coK of Z2")
ggarrange(p1, p2, ncol = 2, nrow = 1)
Simple kriging of a linear combination of the raw variables
err = grd$deleteColumns("SK-LC*")
if(flag.verbose) {
print(paste0(">>> Simple kriging of a linear combination of the variables"))
}
# modelling the linear combination
# Locate Z
err = data$setLocators(nm_LC, locatorType = ELoc_Z(), cleanSameLocator = TRUE)
# computing the model of the LC from the initial model
var_LC = Vario_computeFromDb(varioParam, data)
# plot
multi.varmod(var_LC, mod_LC_SOS, flagLegend=TRUE)
if(flag.verbose) {
mod_LC_SOS
}
# kriging the LC variables
err = data$setLocators(nm_LC, ELoc_Z(), cleanSameLocator=TRUE)
err = kriging(data, grd, mod_LC_SOS, neigh,
flag_est = TRUE, flag_std = TRUE,
namconv = NamingConvention("SK-LC"))
# statistics
vn = c(
paste(c("SK.Z1", "SK-LC.Z1", "SK-LC.LC", "LC.SK"), "estim", sep = "."),
paste(c("SK.Z1", "SK-LC.Z1", "SK-LC.LC"), "stdev", sep = ".")
)
tab = dbStatisticsMono(grd, names = vn, opers = opers)$toTL()
knitr::kable(tab, caption = "Statistics on the grid for Simple Kriging with the LC model", digits = 4)| Number | Minimum | Maximum | Mean | St. Dev. | |
|---|---|---|---|---|---|
| SK.Z1.estim | 10000 | 0.7416 | 1.5086 | 0.9972 | 0.0824 |
| SK-LC.Z1.estim | 10000 | 0.7416 | 1.5086 | 0.9972 | 0.0824 |
| SK-LC.LC.estim | 10000 | -5.2554 | -4.4973 | -5.0028 | 0.0814 |
| LC.SK.estim | 10000 | -5.2554 | -4.4973 | -5.0028 | 0.0814 |
| SK.Z1.stdev | 10000 | 0.1301 | 0.2240 | 0.2078 | 0.0159 |
| SK-LC.Z1.stdev | 10000 | 0.1301 | 0.2240 | 0.2078 | 0.0159 |
| SK-LC.LC.stdev | 10000 | 5.8940 | 5.8968 | 5.8962 | 0.0005 |
# plot Z vs. Z*
p1 = plot.init() +
plot.correlation(grd,
namex = paste("SK", "Z1", "estim", sep = "."),
namey = paste("SK-LC", "Z1", "estim", sep = "."),
asPoint = TRUE, flagDiag = TRUE, diagColor = "red", flagSameAxes = TRUE) +
plot.decoration(xlab = "Initial estimate", ylab = "Estimate with the LC model",
title = "Simple coK of Z1")
p2 = plot.init() +
plot.correlation(grd,
namey = paste("SK-LC", "LC", "estim", sep = "."),
namex = paste("LC", "SK", "estim", sep = "."),
asPoint = TRUE, flagDiag = TRUE, diagColor = "red", flagSameAxes = TRUE) +
plot.decoration(xlab = "LC of the estimates", ylab = "Estimate with the LC model",
title = "Simple coK of CL")
ggarrange(p1, p2, ncol = 2, nrow = 1)
Direct computation of linear combination kriging
err = grd$deleteColumns("dSK*")
if(flag.verbose) {
print(paste0(">>> Direct computation of linear combination kriging"))
}
# the R matrix matLC has to be converted into a gstlearn MatrixDense:
# MatrixDense_createFromVD(beta, nrow=1, ncol=length(beta))
# OR
# beta |> matrix(nrow=1, ncol=length(beta)) |> fromTL()
# OR
# fromTL(matrix(beta, nrow=1, ncol=length(beta)))
# kriging LC variables
err = data$setLocators(nm_var, ELoc_Z(), cleanSameLocator = TRUE)
krigopt = KrigOpt()
krigopt$setMatLC(fromTL(matrix(beta, nrow = 1, ncol = length(beta))))## [1] 0err = kriging(data, grd, mod_Z_SOS, neigh,
flag_est = TRUE, flag_std = TRUE, krigopt=krigopt,
namconv = NamingConvention("dSK"))
# statistics
vn = c(
paste(c("SK-LC.LC", "dSK.LC", "LC.SK"), "estim", sep = "."),
paste(c("SK-LC.LC", "dSK.LC"), "stdev", sep = ".")
)
tab = dbStatisticsMono(grd, names = vn, opers = opers, flagIso = FALSE)$toTL()
knitr::kable(tab, caption = "Statistics on the grid for Simple Kriging of LC", digits = 4)| Number | Minimum | Maximum | Mean | St. Dev. | |
|---|---|---|---|---|---|
| SK-LC.LC.estim | 10000 | -5.2554 | -4.4973 | -5.0028 | 0.0814 |
| dSK.LC.estim | 10000 | -5.2554 | -4.4973 | -5.0028 | 0.0814 |
| LC.SK.estim | 10000 | -5.2554 | -4.4973 | -5.0028 | 0.0814 |
| SK-LC.LC.stdev | 10000 | 5.8940 | 5.8968 | 5.8962 | 0.0005 |
| dSK.LC.stdev | 10000 | 5.8940 | 5.8968 | 5.8962 | 0.0005 |
# plot direct LC estimate vs. estimate using the LC model
p1 = plot.init() +
plot.correlation(grd,
namex = paste("SK-LC", "LC", "estim", sep = "."),
namey = paste("dSK", "LC", "estim", sep = "."),
asPoint = TRUE, flagDiag = TRUE, diagColor = "red", flagSameAxes = TRUE) +
plot.decoration(xlab = "Estimate using LC model", ylab = "Direct estimate of the LC",
title = "Simple coK: estimated value")
p2 = plot.init() +
plot.correlation(grd,
namex = paste("SK-LC", "LC", "stdev", sep = "."),
namey = paste("dSK", "LC", "stdev", sep = "."),
asPoint = TRUE, flagDiag = TRUE, diagColor = "red", flagSameAxes = TRUE) +
plot.decoration(xlab = "Estimate using LC model", ylab = "Direct estimate of the LC",
title = "Simple coK: Std. of kriging error")
ggarrange(p1, p2, ncol = 2, nrow = 1)
Kriging with heterotopic data
Creating heterotopic variables
z1 = data["Z1"]
z1[sample(c(TRUE, FALSE), size = length(z1), prob = c(0.2, 0.8), replace = TRUE)] <- NaN
data["Z1_miss"] <- z1
z2 = data["Z2"]
z2[sample(c(TRUE, FALSE), size = length(z2), prob = c(0.5, 0.5), replace = TRUE)] <- NaN
data["Z2_miss"] <- z2
vn = c("Z1", "Z1_miss", "Z2", "Z2_miss")
tab = dbStatisticsMono(data, names = vn, opers = opers, flagIso = FALSE)$toTL()
rownames(tab) <- vn
knitr::kable(tab, caption = "Statistics on the heterotopic data", digits = 4)| Number | Minimum | Maximum | Mean | St. Dev. | |
|---|---|---|---|---|---|
| Z1 | 50 | 0.6299 | 1.5488 | 1.0274 | 0.2347 |
| Z1_miss | 40 | 0.6299 | 1.5421 | 1.0248 | 0.2219 |
| Z2 | 50 | 0.2570 | 10.3052 | 2.3650 | 1.9771 |
| Z2_miss | 26 | 0.5143 | 7.3356 | 2.3495 | 1.9224 |
Ordinary Kriging of a linear combination with heterotopic data
err = grd$deleteColumns("dOK_He*")
if(flag.verbose) {
print(paste0(">>> Direct computation of linear combination kriging with heterotopic data"))
}
# kriging LC variables
err = data$setLocators(paste(nm_var, "miss", sep = "_"), ELoc_Z(), cleanSameLocator = TRUE)
krigopt = KrigOpt()
krigopt$setMatLC(fromTL(matrix(beta, nrow = 1, ncol = length(beta))))## [1] 0err = kriging(data, grd, mod_Z_IRF0, neigh, krigopt=krigopt,
flag_est = TRUE, flag_std = TRUE,
namconv = NamingConvention("dOK_He"))
# statistics
vn = c(
paste(c("dOK_He.LC", "dOK.LC", "LC.OK"), "estim", sep = "."),
paste(c("dOK_He.LC", "dOK.LC"), "stdev", sep = ".")
)
tab = dbStatisticsMono(grd, names = vn, opers = opers, flagIso = FALSE)$toTL()
knitr::kable(tab, caption = "Statistics on the grid for Ordinary Kriging with heterotopic data",
digits = 4)| Number | Minimum | Maximum | Mean | St. Dev. | |
|---|---|---|---|---|---|
| dOK_He.LC.estim | 10000 | -8.2712 | -2.5472 | -5.3995 | 0.8466 |
| dOK.LC.estim | 10000 | -8.5591 | -3.7712 | -6.0292 | 0.9327 |
| LC.OK.estim | 10000 | -8.5591 | -3.7712 | -6.0292 | 0.9327 |
| dOK_He.LC.stdev | 10000 | 6.0274 | 6.4982 | 6.1254 | 0.0810 |
| dOK.LC.stdev | 10000 | 5.9915 | 6.3041 | 6.0301 | 0.0417 |
# plot direct LC estimate: Heterotopic data vs. Homotopic data
p1 = plot.init() +
plot.correlation(grd,
namex = paste("dOK", "LC", "estim", sep = "."),
namey = paste("dOK_He", "LC", "estim", sep = "."),
asPoint = TRUE, flagDiag = TRUE, diagColor = "red", flagSameAxes = TRUE) +
plot.decoration(xlab = "Homotopic data", ylab = "Heterotopic data",
title = "Ordinary coK of LC: estimated value")
p2 = plot.init() +
plot.correlation(grd,
namex = paste("dOK", "LC", "stdev", sep = "."),
namey = paste("dOK_He", "LC", "stdev", sep = "."),
asPoint = TRUE, flagDiag = TRUE, diagColor = "red", flagSameAxes = TRUE) +
plot.decoration(xlab = "Homotopic data", ylab = "Heterotopic data",
title = "Ordinary coK of LC: Std. of kriging error")
ggarrange(p1, p2, ncol = 2, nrow = 1)
# estimate maps: Heterotopic data vs. Homotopic data
p1 = plot.init(asp=1) +
plot.raster(grd, name = paste("dOK", "LC", "estim", sep = "."),
flagLegend = TRUE, legendName = "Z*", palette = "Spectral",
limits = c(-25,0)) +
plot.decoration(title = "Estimate: homotopic data")
p2 = plot.init(asp=1) +
plot.raster(grd, name = paste("dOK_He", "LC", "estim", sep = "."),
flagLegend = TRUE, legendName = "Z*", palette = "Spectral",
limits = c(-25, 0)) +
plot.decoration(title = "Estimate: heterotopic data")
ggarrange(p1, p2, ncol = 2, nrow = 1)
# Std. maps: Heterotopic data vs. Homotopic data
lim = range(c(grd[paste("dOK", "LC", "stdev", sep = ".")], grd[paste("dOK_He", "LC", "stdev",
sep = ".")]))
p1 = plot.init(asp=1) +
plot.raster(grd, name = paste("dOK", "LC", "stdev", sep = "."),
flagLegend = TRUE, legendName = "Std.", palette = "Spectral", limits = lim) +
plot.decoration(title = "IRF-0: homotopic data")
p2 = plot.init(asp=1) +
plot.raster(grd, name = paste("dOK_He", "LC", "stdev", sep = "."),
flagLegend = TRUE, legendName = "Std.", palette = "Spectral", limits = lim) +
plot.decoration(title = "IRF-0: heterotopic data")
ggarrange(p1, p2, ncol = 2, nrow = 1)
Simple Kriging of a linear combination with heterotopic data
err = grd$deleteColumns("dSK_He*")
if(flag.verbose) {
print(paste0(">>> Direct computation of linear combination simple cokriging with heterotopic data"))
}
# kriging LC variables
err = data$setLocators(paste(nm_var, "miss", sep = "_"), ELoc_Z(), cleanSameLocator = TRUE)
krigopt = KrigOpt()
krigopt$setMatLC(fromTL(matrix(beta, nrow = 1, ncol = length(beta))))## [1] 0err = kriging(data, grd, mod_Z_SOS, neigh, krigopt=krigopt,
flag_est = TRUE, flag_std = TRUE,
namconv = NamingConvention("dSK_He"))
# statistics
vn = c(
paste(c("dSK_He.LC", "dSK.LC", "LC.SK"), "estim", sep = "."),
paste(c("dSK_He.LC", "dSK.LC"), "stdev", sep = ".")
)
tab = dbStatisticsMono(grd, names = vn, opers = opers, flagIso = FALSE)$toTL()
knitr::kable(tab, caption = "Statistics on the grid for simple coKriging with heterotopic data",
digits = 4)| Number | Minimum | Maximum | Mean | St. Dev. | |
|---|---|---|---|---|---|
| dSK_He.LC.estim | 10000 | -5.2200 | -4.6393 | -4.9946 | 0.0576 |
| dSK.LC.estim | 10000 | -5.2554 | -4.4973 | -5.0028 | 0.0814 |
| LC.SK.estim | 10000 | -5.2554 | -4.4973 | -5.0028 | 0.0814 |
| dSK_He.LC.stdev | 10000 | 5.8940 | 5.8968 | 5.8964 | 0.0004 |
| dSK.LC.stdev | 10000 | 5.8940 | 5.8968 | 5.8962 | 0.0005 |
# plot direct LC estimate: Heterotopic data vs. Homotopic data
p1 = plot.init() +
plot.correlation(grd,
namex = paste("dSK", "LC", "estim", sep = "."),
namey = paste("dSK_He", "LC", "estim", sep = "."),
asPoint = TRUE, flagDiag = TRUE, diagColor = "red", flagSameAxes = TRUE) +
plot.decoration(xlab = "Homotopic data", ylab = "Heterotopic data",
title = "Simple coK of LC: estimated value")
p2 = plot.init() +
plot.correlation(grd,
namex = paste("dSK", "LC", "stdev", sep = "."),
namey = paste("dSK_He", "LC", "stdev", sep = "."),
asPoint = TRUE, flagDiag = TRUE, diagColor = "red", flagSameAxes = TRUE) +
plot.decoration(xlab = "Homotopic data", ylab = "Heterotopic data",
title = "Simple coK of LC: Std. of kriging error")
ggarrange(p1, p2, ncol = 2, nrow = 1)
# estimate maps: Heterotopic data vs. Homotopic data
p1 = plot.init(asp=1) +
plot.raster(grd, name = paste("dSK", "LC", "estim", sep = "."),
flagLegend = TRUE, legendName = "Z*", palette = "Spectral",
limits = c(-25,0)) +
plot.decoration(title = "SOS estimate: homotopic data")
p2 = plot.init(asp=1) +
plot.raster(grd, name = paste("dSK_He", "LC", "estim", sep = "."),
flagLegend = TRUE, legendName = "Z*", palette = "Spectral",
limits = c(-25, 0)) +
plot.decoration(title = "SOS estimate: heterotopic data")
ggarrange(p1, p2, ncol = 2, nrow = 1)
# Std. maps: Heterotopic data vs. Homotopic data
lim = range(c(grd[paste("dSK", "LC", "stdev", sep = ".")], grd[paste("dSK_He", "LC", "stdev",
sep = ".")]))
p1 = plot.init(asp=1) +
plot.raster(grd, name = paste("dSK", "LC", "stdev", sep = "."),
flagLegend = TRUE, legendName = "Std.", palette = "Spectral", limits = lim) +
plot.decoration(title = "SOS: homotopic data")
p2 = plot.init(asp=1) +
plot.raster(grd, name = paste("dSK_He", "LC", "stdev", sep = "."),
flagLegend = TRUE, legendName = "Std.", palette = "Spectral", limits = lim) +
plot.decoration(title = "SOS: heterotopic data")
ggarrange(p1, p2, ncol = 2, nrow = 1)
QC OK vs. SK
v = "Z1"
# cross plot: SOS vs. IRF0
p1 = plot.init() +
plot.correlation(grd,
namex = paste("SK", v, "estim", sep = "."),
namey = paste("OK", v, "estim", sep = "."),
asPoint = TRUE, flagDiag = TRUE, diagColor = "red", flagSameAxes = TRUE) +
plot.decoration(xlab = "SOS model (SK)", ylab = "IRF-0 model (OK)",
title = paste0(v, " Estimation"))
p2 = plot.init() +
plot.correlation(grd,
namex = paste("SK", v, "stdev", sep = "."),
namey = paste("OK", v, "stdev", sep = "."),
asPoint = TRUE, flagDiag = TRUE, diagColor = "red", flagSameAxes = TRUE) +
plot.decoration(xlab = "SOS model (SK)", ylab = "IRF-0 model (OK)",
title = "Std. of estimation error")
ggarrange(p1, p2, ncol = 2, nrow = 1)
v = "Z2"
# cross plot: SOS vs. IRF0
p1 = plot.init() +
plot.correlation(grd,
namex = paste("SK", v, "estim", sep = "."),
namey = paste("OK", v, "estim", sep = "."),
asPoint = TRUE, flagDiag = TRUE, diagColor = "red", flagSameAxes = TRUE) +
plot.decoration(xlab = "SOS model (SK)", ylab = "IRF-0 model (OK)",
title = paste0(v, " Estimation"))
p2 = plot.init() +
plot.correlation(grd,
namex = paste("SK", v, "stdev", sep = "."),
namey = paste("OK", v, "stdev", sep = "."),
asPoint = TRUE, flagDiag = TRUE, diagColor = "red", flagSameAxes = TRUE) +
plot.decoration(xlab = "SOS model (SK)", ylab = "IRF-0 model (OK)",
title = "Std. of estimation error")
ggarrange(p1, p2, ncol = 2, nrow = 1)
# direct kriging of LC
v = "LC"
# cross plot: SOS vs. IRF0
p1 = plot.init() +
plot.correlation(grd,
namex = paste("dSK", v, "estim", sep = "."),
namey = paste("dOK", v, "estim", sep = "."),
asPoint = TRUE, flagDiag = TRUE, diagColor = "red", flagSameAxes = TRUE) +
plot.decoration(xlab = "SOS model (SK)", ylab = "IRF-0 model (OK)",
title = "LC estimation (direct)")
p2 = plot.init() +
plot.correlation(grd,
namex = paste("dSK", v, "stdev", sep = "."),
namey = paste("dOK", v, "stdev", sep = "."),
asPoint = TRUE, flagDiag = TRUE, diagColor = "red", flagSameAxes = TRUE) +
plot.decoration(xlab = "SOS model (SK)", ylab = "IRF-0 model (OK)",
title = "Std. of LC estimation error (direct)")
ggarrange(p1, p2, ncol = 2, nrow = 1)