# R/turningLayers.R In gstat: Spatial and Spatio-Temporal Geostatistical Modelling, Prediction and Simulation

#### Documented in krigeSTSimTB

```#
# library(gstat)
# # separable model: spatial and temporal sill will be ignored
# # and kept constant at 1-nugget respectively. A joint sill is used.
# separableModel <- vgmST("separable",
#                         space=vgm(0.9,"Exp", 147, 0.1),
#                         time =vgm(0.9,"Exp", 3.5, 0.1),
#                         sill=40)
#
# covSTMat <- ceWrapSpaceTimeOnTorusCalcCovRow1(c(250,20), c(1,16), separableModel)
# dim(covSTMat)
#
# image(matrix(ceSim(covSTMat, 1, c(20, 16)),
#              nrow=20, ncol=16), asp=1, ylab="time", xlab="space")
#
# # product sum model: spatial and temporal nugget will be ignored and kept
# # constant at 0. Only a joint nugget is used.
# prodSumModel <- vgmST("productSum",
#                       space=vgm(39, "Sph", 343, 0),
#                       time= vgm(36, "Exp",   3, 0),
#                       k=15)
#
# covSTMat <- ceWrapSpaceTimeOnTorusCalcCovRow1(c(250,20), c(1,16), prodSumModel)
# dim(covSTMat)
#
# image(matrix(ceSim(covSTMat, 1, c(20, 16)),
#              nrow=20, ncol=16), asp=1, ylab="time", xlab="space")
#
# # sum metric model: spatial, temporal and joint nugget will be estimated
# sumMetricModel <- vgmST("sumMetric",
#                         space=vgm( 6.9, "Lin", 200, 3.0),
#                         time =vgm(10.3, "Lin",  15, 3.6),
#                         joint=vgm(37.2, "Exp",  84,11.7),
#                         stAni=77.7)
#
# covSTMat <- ceWrapSpaceTimeOnTorusCalcCovRow1(c(250,20), c(1,16), sumMetricModel)
# dim(covSTMat)
#
# image(matrix(ceSim(covSTMat, 1, c(20, 16)),
#              nrow=20, ncol=16), asp=1, ylab="time", xlab="space")
#
# # simplified sumMetric model, only a overall nugget is fitted. The spatial,
# # temporal and jont nuggets are set to 0.
# simpleSumMetricModel <- vgmST("simpleSumMetric",
#                               space=vgm(20,"Lin", 150, 0),
#                               time =vgm(20,"Lin", 10,  0),
#                               joint=vgm(20,"Exp", 150, 0),
#                               nugget=1, stAni=15)
#
# covSTMat <- ceWrapSpaceTimeOnTorusCalcCovRow1(c(250,20), c(1,16), simpleSumMetricModel)
# dim(covSTMat)
#
# image(matrix(ceSim(covSTMat, 1, c(20, 16)),
#              nrow=20, ncol=16), asp=1, ylab="time", xlab="space")
#
# # metric model
# metricModel <- vgmST("metric",
#                      joint=vgm(60, "Exp", 150, 10),
#                      stAni=60)
#
# covSTMat <- ceWrapSpaceTimeOnTorusCalcCovRow1(c(250,20), c(1,16), metricModel, turningLayers = FALSE)
# dim(covSTMat)
#
# image(matrix(ceSim(covSTMat, 1, c(20, 16)),
#              nrow=20, ncol=16), asp=1, ylab="time", xlab="space")

## turning bands:
# for cov-functions valid in 3-dim

# adjust the covariance matrix (Schlahter, Chap. 2, eq 2.25)
# simulate several layers (up-right hyper-planes in the 3D+T-cube) with random directions on a unit hemi-sphere
# orthogonally project target points onto rndm lyrs and average the values for all rndm lyrs

# -> one simulation of a 3d+T Gaussian random field, for 2D use a hyperplane in 3D

# the tunring bands operator (numerical, 3-dimensional case)
# @input:
#  model target covariance function
#  dist_grid: data.frame wtih columns spacelag and timelag for which the covariances are calculated
#
# @value: data.frame with columns spacelag, timelag and gamma where the latter contains the covariances for the 1-dim spatial and temporal layers

# model <- separableModel
# dist_grid <- as.data.frame(cbind("spacelag" = rep(1:150*1., each=4),
#                                  "timelag" = rep(1:4, 150)))

# CAVE: distgrid must ahve the correct spatial and temporal metrics

tbOperator <- function(model, dist_grid) {
r <- dist_grid\$spacelag

derFun <- function(r) r * variogramSurface(model, dist_grid = data.frame(spacelag=r, timelag=dist_grid\$timelag),
covariance = TRUE)\$gamma

cbind(dist_grid, "gamma" = diag(attr(numericDeriv(quote(derFun(r)), "r"), "gradient")))
}

randomDirections <- function(n) {
u <- runif(n, 0, 2*pi)
v <- runif(n, 0, 1)
s <- sqrt(1-v^2)

cbind(s*cos(u), s*sin(u), v)
}

# library(rgl)
# plot3d(randomDirections(100), aspect = c(1,1,0.5))

# nLyrs <- 500
# coordMat <- matrix(c(1,5,5,
#                      1,7,7,
#                      1,5,7,
#                      1,7,5), byrow=T, ncol=3)
#
# coordMat <- matrix(runif(3*1001,1,5), byrow=T, ncol=3)
#
# # random directions in 3D
# rndDir <- randomDirections(nLyrs)
#
# # how much does each direction contribute to the point
# cntrbtn <- rndDir %*% t(coordMat)
# # -> each column corresponds to one location, each row i to the "contribution" of the i-th rndDir
#
# cl_cntrbtn <- ceiling(cntrbtn)+10 # centre for grid; mind spatial index
# fl_cntrbtn <- floor(cntrbtn) +10
#
# covRow1 <- ceWrapSpaceTimeOnTorusCalcCovRow1(c(250,20), c(1,16), metricModel, turningLayers = TRUE)
#
# sTime <- Sys.time()
# simLyrs <- ceSim(covRow1, nLyrs, c(20,16))
# simLyrs <- lapply(1:nLyrs, function(col) matrix(simLyrs[,col], nrow=20, ncol=16))
#
# lambda <- cl_cntrbtn - cntrbtn -10
#
# cntrbSngSimLyr <- function(lyrId) {
#   simLyrs[[lyrId]][cl_cntrbtn[lyrId,],] * (1-lambda)[lyrId,] + simLyrs[[lyrId]][fl_cntrbtn[lyrId,],] * lambda[lyrId,]
# }
#
# # reduce to the one realisation based on the combination of nLyrs turning bands
# simTs <- Reduce('+', lapply(1:nLyrs, cntrbSngSimLyr))/sqrt(nLyrs)
# eTime <- Sys.time()
#
# eTime - sTime # 0.7 secs/simulation of 100 random points
#
# image(simTs)
#
# dim(simTs)
# # 1001 locations at 16 time stamps
#
# plot(simTs[500,])

# computes the covariance matrixes and weights once, applied to series of
# variables/simulations where each variable/simulation is stored in one column of
# the multiVarMatrix copied from krigeST to avoid repeted calls to krige with
# multiple, identical inversions of the weights matrix

# TODO: add functionality for temporal sandwich-wise processing: i.e. use the
# +/- nmaxTime time slices to predict one time slice

krigeSTMultiple <- function(formula, from, to, modelList, multiVarMatrix, nmaxTime=Inf) {
lst = extractFormula(formula, from, to)

separate <- length(from) > 1 && length(to) > 1 &&
inherits(from, "STF") && inherits(to, "STF")

X = lst\$X
x0 = lst\$x0

V = covfn.ST(from, model = modelList, separate=separate)
v0 = covfn.ST(from, to, modelList)

if (modelList\$stModel == "separable" & separate)
skwts <- STsolve(V, v0, X) # use Kronecker trick
else
skwts <- CHsolve(V, cbind(v0, X))

npts = length(to)
ViX = skwts[,-(1:npts)]
skwts = skwts[,1:npts]

idPredFun <- function(sim) {
sim <- matrix(sim, ncol = 1)
beta = solve(t(X) %*% ViX, t(ViX) %*% sim)
x0 %*% beta + t(skwts) %*% (sim - X %*% beta)
}

apply(multiVarMatrix, 2, idPredFun)
}

# unconcitional
## STFDF
# library(sp)
# library(spacetime)
# library(zoo)
# library(xts)
# data(meuse, package = "sp")
# coordinates(meuse) <- ~x+y
# proj4string(meuse) <- CRS("+init=epsg:28992")

# krigeST <- function(formula, data, newdata, modelList, y, beta, nmax=Inf, stAni=NULL,
#                     computeVar = FALSE, fullCovariance = FALSE,
#                     bufferNmax=2, progress=TRUE)

##
krigeSTSimTB <- function(formula, data, newdata, modelList, nsim,
progress=TRUE, nLyrs=500, tGrid=NULL, sGrid=NULL, ceExt=2,
nmax=Inf) {
stopifnot(zoo::is.regular(newdata@time))

condSim <- TRUE
if (missing(data)) {
condSim <- FALSE
message("[No data provided: performing unconditional simulation.]")
} else {
message("[Performing conditional simulation.]")
}

pb <- txtProgressBar(0,nsim,style=3)

# ST-simulation grid
if (is.null(tGrid)) {
tDis <- diff(c(index(newdata@time[1]), newdata@endTime[1]))
if (!is.null(attr(modelList, "temporal unit"))) {
units(tDis) <- attr(modelList, "temporal unit")
} else {
message("[The spatio-temporal variogram model does not carry a time unit attribute: krigeST cannot check whether the temporal distance metrics coincide.]")
}

tGrid <- c(as.numeric(tDis), length(newdata@time))
attr(tGrid, "units") <- c("", units(tDis))

debug_time_unit(units(tDis))
}

if (is.null(sGrid)) {
if (gridded(newdata@sp)) {
# SpatialPixels/SpatialGrid:
# based on GridTopology: use minimal cellsize of both directions; take enough to cover the diagonal
sDis <- min(newdata@sp@grid@cellsize)
sDim <- ceiling(sqrt(sum((newdata@sp@grid@cellsize * newdata@sp@grid@cells.dim)^2))/sDis)
sGrid <-  c(sDis, sDim)
} else {
# treat (as) SpatialPoints:
# Average area per location --assuming-a-regular-squared-outline-taking-the-sqrt--> length/location --take-later-twice-as-many-->
bboxExt <- apply(newdata@sp@bbox, 1, diff)
sDis <- sqrt(prod(bboxExt)/length(newdata@sp))
sDim <- ceiling(sqrt(sum((bboxExt/sDis)^2)))
sGrid <- c(sDis, sDim)
}
}

# random directions in 3D
rndDir <- randomDirections(nLyrs)

# coordinates
coordMat <- coordinates(newdata@sp)

# coordinates (embedded in 3D) shifted + scaled to the grid index of the spatial gridding sGrid
coordMat[,1] <- (coordMat[,1] - newdata@sp@bbox[1,1])/sGrid[1]
coordMat[,2] <- (coordMat[,2] - newdata@sp@bbox[2,1])/sGrid[1]

if (ncol(coordMat) == 2) {
coordMat <- cbind(coordMat,1)
} else {
coordMat[,3] <- (coordMat[,3] - newdata@sp@bbox[3,1])/sGrid[1]
}

# how much does each direction contribute to the point
cntrbtn <- rndDir %*% t(coordMat)
# -> each column corresponds to one location, each row i to the "contribution" of the i-th rndDir

# ceiling and floor + shifted to avoid negative indices
cl_cntrbtn <- ceiling(cntrbtn) + sGrid[2]
fl_cntrbtn <- floor(cntrbtn) + sGrid[2]

covRow1 <- ceWrapSpaceTimeOnTorusCalcCovRow1(c(sGrid[1], 2*sGrid[2]), tGrid, modelList, turningLayers = TRUE, ext=ceExt)

origDim <- c(2*sGrid[2], tGrid[2])
sims <- list()

##
for (i in 1:nsim) {
setTxtProgressBar(pb, i)
simLyrs <- ceSim(covRow1, nLyrs, origDim)
simLyrs <- lapply(1:nLyrs, function(col) matrix(simLyrs[,col],
nrow=origDim[1],
ncol=origDim[2]))
lambda <- cl_cntrbtn - cntrbtn - sGrid[2]

cntrbSngSimLyr <- function(lyrId) {
simLyrs[[lyrId]][cl_cntrbtn[lyrId,],] * (1-lambda)[lyrId,] + simLyrs[[lyrId]][fl_cntrbtn[lyrId,],] * lambda[lyrId,]
}

# reduce to the one realisation based on the combination of nLyrs turning bands
sims[[paste0("sim",i)]] <- Reduce('+', lapply(1:nLyrs, cntrbSngSimLyr))/sqrt(nLyrs)
}
close(pb)

sims <- do.call(cbind, lapply(sims, as.numeric))

# bind simulations to newdata geometry
if (!condSim) {
if ("data" %in% slotNames(newdata))
newdata@data <- cbind(newdata@data, sims)
else
return(newdata)
}

# function call ends here if no data has been provided -> unconditional case
varName <- all.vars(formula[[2]])

## conditioning
# interpolate the observations to the simulation grid
obsMeanField <- krigeST(formula=formula, data=data, newdata=newdata, modelList=modelList)

# interpolate to observation locations from the simulated grids for each simulation
simMeanObsLoc <- krigeSTMultiple(as.formula(paste0("var1.pred ~", formula[[3]])),
obsMeanField, data, modelList, sims)

# interpolate from kriged mean sim at observed locations back to the grid for mean surface of the simulations
simMeanFields <- krigeSTMultiple(as.formula(paste0(varName, "~", formula[[3]])),
data, newdata, modelList, simMeanObsLoc)

# add up the mean field and the corrected data
sims <- obsMeanField@data\$var1.pred + sims - simMeanFields

# bind simulations to newdata geometry
if ("data" %in% slotNames(newdata)) {
newdata@data <- cbind(newdata@data, sims)
return(newdata)
}

}

#
# sTime <- Sys.time()
# krigedSim <- krigeSTUncSimTB(stf, metricModel, 100)
# Sys.time() - sTime
#
# # 27 secs for 100 simulated ST fields of 155 locations and 21 time steps: 325500 values
#
# # plot one simulation along time
# stplot(krigedSim[,1:12])
#
# # plot one simulation along time
# stplot(krigedSim[1:12,,"sim1"], mode="ts")
#
# # plot the ten simulations of the first day
# spplot(krigedSim[,1], paste0("sim",1:10), as.table=TRUE)
```

## Try the gstat package in your browser

Any scripts or data that you put into this service are public.

gstat documentation built on April 6, 2023, 5:21 p.m.