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:

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.1", "Y1")
err = grd$setName("Simu.2", "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$getSampleNumber(), 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 = ggDefaultGeographic() + 
  plot.grid(grd, nameRaster = "Z1", palette = "Spectral") +
  plot.point(data, colour = "red") +
  plot.decoration(title="Z1")
p2 = ggDefault() +
  plot.hist(grd, name='Z1',  bins = 25, bg = "orange", color= "gray") +
  plot.decoration(xlab = "Raw variable", title="Histogram of the initial variable Z1")
p3 = ggDefaultGeographic() + 
  plot.grid(grd, nameRaster = "Z2", palette = "Spectral") +
  plot.point(data, colour = "red") +
  plot.decoration(title="Z2")
p4 = ggDefault() +
  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$getSampleNumber(), 
                    ncol = dim(matLC)[2]) %*% t(matLC)[,1]
grd["LC"]  = matrix(grd$getColumns(colnames(matLC)),  nrow = grd$getSampleNumber(),  
                    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)
Statistics on the reference data set
Number Minimum Maximum Mean St. Dev.
grid-Z1 10000 0.4528 2.2903 1.0419 0.2619
grid-Z2 10000 0.1181 15.8685 1.9705 1.7018
grid-LC 10000 -47.0365 1.5622 -4.8696 5.2569
samples-Z1 50 0.7668 1.6995 1.1311 0.2404
samples-Z2 50 0.2648 4.6042 1.6170 1.1103
samples-LC 50 -12.9096 0.7581 -3.7199 3.4745

Global parameters for spatial modelling

# Variogram parameters (omnidirectional)
varioParam = VarioParam_createOmniDirection(npas=10, dpas=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$getVariableNumber()
  stopifnot(nv == dim(A)[1])
  ns = model$getCovaNumber()
  model_LC = model$clone()
  for (k in 1:ns) {
    CC = A %*% model$getSillValues(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$getSampleNumber(),  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)
Statistics on the grid for Ordinary Kriging with the initial model
Number Minimum Maximum Mean St. Dev.
Z1 10000 0.4528 2.2903 1.0419 0.2619
OK.Z1.estim 10000 0.7978 1.6791 1.1037 0.1053
OK.Z1.stdev 10000 0.0443 0.2663 0.2206 0.0273
Z2 10000 0.1181 15.8685 1.9705 1.7018
OK.Z2.estim 10000 0.2454 3.4234 1.7027 0.6338
OK.Z2.stdev 10000 0.7099 1.2544 0.8919 0.0659 
# plot of estimated values
p1 = ggDefaultGeographic() + 
    plot.grid(grd, nameRaster =  paste("OK", "Z1", "estim", sep = "."),
            palette = "Spectral", flagLegendRaster = TRUE, legendNameRaster = "estim") +
    plot.decoration(xlab = "", ylab = "", title = paste0("Ordinary coK for Z1"))
p2 = ggDefaultGeographic() + 
    plot.grid(grd, nameRaster =  paste("OK", "Z1", "stdev", sep = "."),
            palette = "Spectral", flagLegendRaster = TRUE, legendNameRaster = "stdev") +
    plot.decoration(xlab = "", ylab = "", title = paste0("Std. of Ordinary coK for Z1"))
p3 = ggDefaultGeographic() + 
    plot.grid(grd, nameRaster =  paste("OK", "Z2", "estim", sep = "."),
            palette = "Spectral", flagLegendRaster = TRUE, legendNameRaster = "estim") +
    plot.decoration(xlab = "", ylab = "", title = paste0("Ordinary coK for Z2"))
p4 = ggDefaultGeographic() + 
    plot.grid(grd, nameRaster =  paste("OK", "Z2", "stdev", sep = "."),
            palette = "Spectral", flagLegendRaster = TRUE, legendNameRaster = "stdev") +
    plot.decoration(xlab = "", ylab = "", title = paste0("Std. of Ordinary coK for Z2"))
ggarrange(p1, p2, p3, p4, ncol = 2, nrow = 2)
## Warning: `guide_colourbar()` cannot be used for fill_ggnewscale_1.
## ℹ Use one of colour, color, or fill instead.
## `guide_colourbar()` cannot be used for fill_ggnewscale_1.
## ℹ Use one of colour, color, or fill instead.
## `guide_colourbar()` cannot be used for fill_ggnewscale_1.
## ℹ Use one of colour, color, or fill instead.
## `guide_colourbar()` cannot be used for fill_ggnewscale_1.
## ℹ Use one of colour, color, or fill instead.

# plot Z vs. Z*
p1 = ggDefault() + 
  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 = ggDefault() + 
  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)
Statistics on the grid for Ordinary Kriging with the LC model
Number Minimum Maximum Mean St. Dev.
OK.Z1.estim 10000 0.7978 1.6791 1.1037 0.1053
OK-LC.Z1.estim 10000 0.7978 1.6791 1.1037 0.1053
OK-LC.LC.estim 10000 -9.3628 0.7857 -4.0043 1.9728
LC.OK.estim 10000 -9.3628 0.7857 -4.0043 1.9728
OK.Z1.stdev 10000 0.0443 0.2663 0.2206 0.0273
OK-LC.Z1.stdev 10000 0.0443 0.2663 0.2206 0.0273
OK-LC.LC.stdev 10000 2.1581 3.9132 2.7951 0.2130 
# plot Z vs. Z*
p1 = ggDefault() + 
  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 = ggDefault() + 
  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 MatrixRectangular:
# MatrixRectangular_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)
err = kriging(data, grd, mod_Z_IRF0, neigh, 
              matLC = fromTL(matrix(beta, nrow=1, ncol=length(beta))),
              flag_est = TRUE, flag_std = TRUE,
              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)
Statistics on the grid for Ordinary Kriging of LC
Number Minimum Maximum Mean St. Dev.
OK-LC.LC.estim 10000 -9.3628 0.7857 -4.0043 1.9728
dOK.LC.estim 10000 -9.3628 0.7857 -4.0043 1.9728
LC.OK.estim 10000 -9.3628 0.7857 -4.0043 1.9728
OK-LC.LC.stdev 10000 2.1581 3.9132 2.7951 0.2130
dOK.LC.stdev 10000 2.1581 3.9132 2.7951 0.2130 
# plot direct LC estimate vs. estimate using the LC model
p1 = ggDefault() + 
  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 = ggDefault() + 
  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$getSampleNumber(),  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)
Statistics on the grid for Simple Kriging with the initial model
Number Minimum Maximum Mean St. Dev.
Z1 10000 0.4528 2.2903 1.0419 0.2619
SK.Z1.estim 10000 0.7961 1.6773 1.0570 0.1068
SK.Z1.stdev 10000 0.0443 0.2466 0.2177 0.0257
Z2 10000 0.1181 15.8685 1.9705 1.7018
SK.Z2.estim 10000 0.2331 3.4294 1.7523 0.6352
SK.Z2.stdev 10000 0.7099 1.2486 0.8910 0.0651 
# plot of estimated values
p1 = ggDefaultGeographic() + 
    plot.grid(grd, nameRaster =  paste("SK", "Z1", "estim", sep = "."),
            palette = "Spectral", flagLegendRaster = TRUE, legendNameRaster = "estim") +
    plot.decoration(xlab = "", ylab = "", title = paste0("Simple coK for Z1"))
p2 = ggDefaultGeographic() + 
    plot.grid(grd, nameRaster =  paste("SK", "Z1", "stdev", sep = "."),
            palette = "Spectral", flagLegendRaster = TRUE, legendNameRaster = "stdev") +
    plot.decoration(xlab = "", ylab = "", title = paste0("Std. of Simple coK for Z1"))
p3 = ggDefaultGeographic() + 
    plot.grid(grd, nameRaster =  paste("SK", "Z2", "estim", sep = "."),
            palette = "Spectral", flagLegendRaster = TRUE, legendNameRaster = "estim") +
    plot.decoration(xlab = "", ylab = "", title = paste0("Simple coK for Z2"))
p4 = ggDefaultGeographic() + 
    plot.grid(grd, nameRaster =  paste("SK", "Z2", "stdev", sep = "."),
            palette = "Spectral", flagLegendRaster = TRUE, legendNameRaster = "stdev") +
    plot.decoration(xlab = "", ylab = "", title = paste0("Std. of Simple coK for Z2"))
ggarrange(p1, p2, p3, p4, ncol = 2, nrow = 2)
## Warning: `guide_colourbar()` cannot be used for fill_ggnewscale_1.
## ℹ Use one of colour, color, or fill instead.
## `guide_colourbar()` cannot be used for fill_ggnewscale_1.
## ℹ Use one of colour, color, or fill instead.
## `guide_colourbar()` cannot be used for fill_ggnewscale_1.
## ℹ Use one of colour, color, or fill instead.
## `guide_colourbar()` cannot be used for fill_ggnewscale_1.
## ℹ Use one of colour, color, or fill instead.

# plot Z vs. Z*
p1 = ggDefault() + 
  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 = ggDefault() + 
  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)
Statistics on the grid for Simple Kriging with the LC model
Number Minimum Maximum Mean St. Dev.
SK.Z1.estim 10000 0.7961 1.6773 1.0570 0.1068
SK-LC.Z1.estim 10000 0.7961 1.6773 1.0570 0.1068
SK-LC.LC.estim 10000 -9.3904 0.8032 -4.2000 1.9793
LC.SK.estim 10000 -9.3904 0.8032 -4.2000 1.9793
SK.Z1.stdev 10000 0.0443 0.2466 0.2177 0.0257
SK-LC.Z1.stdev 10000 0.0443 0.2466 0.2177 0.0257
SK-LC.LC.stdev 10000 2.1580 3.8860 2.7909 0.2096 
# plot Z vs. Z*
p1 = ggDefault() + 
  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 = ggDefault() + 
  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 MatrixRectangular:
# MatrixRectangular_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)
err = kriging(data, grd, mod_Z_SOS, neigh, 
              matLC = fromTL(matrix(beta, nrow=1, ncol=length(beta))),
              flag_est = TRUE, flag_std = TRUE,
              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)
Statistics on the grid for Simple Kriging of LC
Number Minimum Maximum Mean St. Dev.
SK-LC.LC.estim 10000 -9.3904 0.8032 -4.2000 1.9793
dSK.LC.estim 10000 -9.3904 0.8032 -4.2000 1.9793
LC.SK.estim 10000 -9.3904 0.8032 -4.2000 1.9793
SK-LC.LC.stdev 10000 2.1580 3.8860 2.7909 0.2096
dSK.LC.stdev 10000 2.1580 3.8860 2.7909 0.2096 
# plot direct LC estimate vs. estimate using the LC model
p1 = ggDefault() + 
  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 = ggDefault() + 
  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)
Statistics on the heterotopic data
Number Minimum Maximum Mean St. Dev.
Z1 50 0.7668 1.6995 1.1311 0.2404
Z1_miss 40 0.7668 1.6995 1.1455 0.2535
Z2 50 0.2648 4.6042 1.6170 1.1103
Z2_miss 26 0.2648 4.1471 1.6766 1.0977

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)
err = kriging(data, grd, mod_Z_IRF0, neigh, 
              matLC = fromTL(matrix(beta, nrow=1, ncol=length(beta))),
              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)
Statistics on the grid for Ordinary Kriging with heterotopic data
Number Minimum Maximum Mean St. Dev.
dOK_He.LC.estim 10000 -10.1558 0.8856 -4.1178 1.7557
dOK.LC.estim 10000 -9.3628 0.7857 -4.0043 1.9728
LC.OK.estim 10000 -9.3628 0.7857 -4.0043 1.9728
dOK_He.LC.stdev 10000 2.2084 3.9516 2.9518 0.2158
dOK.LC.stdev 10000 2.1581 3.9132 2.7951 0.2130 
# plot direct LC estimate: Heterotopic data vs. Homotopic data
p1 = ggDefault() + 
  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 = ggDefault() + 
  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 = ggDefaultGeographic() + 
  plot.grid(grd, nameRaster =  paste("dOK", "LC", "estim", sep = "."), 
            flagLegendRaster = TRUE, legendNameRaster = "Z*", palette = "Spectral", 
            limits = c(-25,0)) +
  plot.decoration(title = "Estimate: homotopic data")
p2 = ggDefaultGeographic() + 
  plot.grid(grd, nameRaster =  paste("dOK_He", "LC", "estim", sep = "."),
            flagLegendRaster = TRUE, legendNameRaster = "Z*", palette = "Spectral", 
            limits = c(-25, 0)) +
  plot.decoration(title = "Estimate: heterotopic data")
ggarrange(p1, p2, ncol = 2, nrow = 1)
## Warning: `guide_colourbar()` cannot be used for fill_ggnewscale_1.
## ℹ Use one of colour, color, or fill instead.
## `guide_colourbar()` cannot be used for fill_ggnewscale_1.
## ℹ Use one of colour, color, or fill instead.

# Std. maps: Heterotopic data vs. Homotopic data
lim = range(c(grd[paste("dOK", "LC", "stdev", sep = ".")], grd[paste("dOK_He", "LC", "stdev", 
                                                                     sep = ".")]))
p1 = ggDefaultGeographic() + 
  plot.grid(grd, nameRaster =  paste("dOK", "LC", "stdev", sep = "."), 
            flagLegendRaster = TRUE, legendNameRaster = "Std.", palette = "Spectral", limits = lim) +
  plot.decoration(title = "IRF-0: homotopic data")
p2 = ggDefaultGeographic() + 
  plot.grid(grd, nameRaster =  paste("dOK_He", "LC", "stdev", sep = "."),
            flagLegendRaster = TRUE, legendNameRaster = "Std.", palette = "Spectral", limits = lim) +
  plot.decoration(title = "IRF-0: heterotopic data")
ggarrange(p1, p2, ncol = 2, nrow = 1)
## Warning: `guide_colourbar()` cannot be used for fill_ggnewscale_1.
## ℹ Use one of colour, color, or fill instead.
## `guide_colourbar()` cannot be used for fill_ggnewscale_1.
## ℹ Use one of colour, color, or fill instead.

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)
err = kriging(data, grd, mod_Z_SOS, neigh, 
              matLC = fromTL(matrix(beta, nrow=1, ncol=length(beta))),
              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)
Statistics on the grid for simple coKriging with heterotopic data
Number Minimum Maximum Mean St. Dev.
dSK_He.LC.estim 10000 -10.2113 1.0058 -4.1182 1.7736
dSK.LC.estim 10000 -9.3904 0.8032 -4.2000 1.9793
LC.SK.estim 10000 -9.3904 0.8032 -4.2000 1.9793
dSK_He.LC.stdev 10000 2.2068 3.9388 2.9489 0.2150
dSK.LC.stdev 10000 2.1580 3.8860 2.7909 0.2096 
# plot direct LC estimate: Heterotopic data vs. Homotopic data
p1 = ggDefault() + 
  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 = ggDefault() + 
  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 = ggDefaultGeographic() + 
  plot.grid(grd, nameRaster =  paste("dSK", "LC", "estim", sep = "."), 
            flagLegendRaster = TRUE, legendNameRaster = "Z*", palette = "Spectral", 
            limits = c(-25,0)) +
  plot.decoration(title = "SOS estimate: homotopic data")
p2 = ggDefaultGeographic() + 
  plot.grid(grd, nameRaster =  paste("dSK_He", "LC", "estim", sep = "."),
            flagLegendRaster = TRUE, legendNameRaster = "Z*", palette = "Spectral", 
            limits = c(-25, 0)) +
  plot.decoration(title = "SOS estimate: heterotopic data")
ggarrange(p1, p2, ncol = 2, nrow = 1)
## Warning: `guide_colourbar()` cannot be used for fill_ggnewscale_1.
## ℹ Use one of colour, color, or fill instead.
## `guide_colourbar()` cannot be used for fill_ggnewscale_1.
## ℹ Use one of colour, color, or fill instead.

# Std. maps: Heterotopic data vs. Homotopic data
lim = range(c(grd[paste("dSK", "LC", "stdev", sep = ".")], grd[paste("dSK_He", "LC", "stdev", 
                                                                     sep = ".")]))
p1 = ggDefaultGeographic() + 
  plot.grid(grd, nameRaster =  paste("dSK", "LC", "stdev", sep = "."), 
            flagLegendRaster = TRUE, legendNameRaster = "Std.", palette = "Spectral", limits = lim) +
  plot.decoration(title = "SOS: homotopic data")
p2 = ggDefaultGeographic() + 
  plot.grid(grd, nameRaster =  paste("dSK_He", "LC", "stdev", sep = "."),
            flagLegendRaster = TRUE, legendNameRaster = "Std.", palette = "Spectral", limits = lim) +
  plot.decoration(title = "SOS: heterotopic data")
ggarrange(p1, p2, ncol = 2, nrow = 1)
## Warning: `guide_colourbar()` cannot be used for fill_ggnewscale_1.
## ℹ Use one of colour, color, or fill instead.
## `guide_colourbar()` cannot be used for fill_ggnewscale_1.
## ℹ Use one of colour, color, or fill instead.

QC OK vs. SK

v = "Z1"

# cross plot: SOS vs. IRF0
p1 = ggDefault() + 
  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 = ggDefault() + 
  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 = ggDefault() + 
  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 = ggDefault() + 
  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 = ggDefault() + 
  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 = ggDefault() + 
  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)