3d_kriging_cgs_cosimulation.R

library(gstat)
library(sp)
library(lattice)
library(data.table)
library(ggplot2)

# Create Data Points (Random)
n <- 50
data3D <- data.frame(x = runif(n), y = runif(n), z = runif(n), v = rnorm(n))
coordinates(data3D) = ~x+y+z

# Create Pseudo Wells
n <- 5
data3D <- data.table(merge(data.frame(x = runif(n),y = runif(n)),data.frame(z = seq(0,1,by=0.05))))
data3D[,avg_phie:=rnorm(1,mean = 0.30,sd = 0.02),by=c("x","y")]
data3D[,avg_swe:=rnorm(1,mean = 0.27,sd = 0.04),by=c("x","y")]
data3D[,avg_h:=runif(1)*(1-0.2)+0.2,by=c("x","y")]
data3D[,phie:=rnorm(.N,mean=avg_phie,sd=0.02),by=c("x","y")]
data3D[,swe:=rnorm(.N,mean=avg_swe,sd=0.04),by=c("x","y")]
data3D[z>avg_h,phie:=0,]
data3D <- data.frame(subset(data3D,select=c(x,y,z,phie,swe)))
coordinates(data3D) = ~x+y+z

# Build a gstat structure for co-kriging
model <- gstat(NULL, id = "phie", form = phie ~ 1, data=data3D, beta = 0, nmax=10)
model <- gstat(model, id = "swe", form = swe ~ 1, data=data3D, beta = 1, nmax=10)
model <- gstat(model, id = "phie", model = vgm(1,"Exp",0.2),fill.all=TRUE)
v.cross <- variogram(model)
str(v.cross)
plot(v.cross, pl=T)
model <- fit.lmc(v.cross, model)
plot(variogram(model), model = model$model)

# Create empty grid to krige
range1D <- seq(from = 0, to = 1, length = 21)
grid3D <- expand.grid(x = range1D, y = range1D, z = range1D)
gridded(grid3D) = ~x+y+z

res3D <- predict(model, newdata = grid3D,nsim=100)

# Plot Results
levelplot(phie.sim1 ~ x + y | z, as.data.frame(res3D))

# HPV Assessment:

resData <- data.table(res3D@data)

temp <- data.table(melt(resData))
temp[,variable:=as.character(variable)]
temp[,cell:=seq(1,.N,by=1),by=c("variable")]
temp[,sim:=strsplit(variable,"\\.")[[1]][2],by=c("variable")]
temp[,variable:=strsplit(variable,"\\.")[[1]][1],by=c("variable")]

temp <- data.table(dcast(temp,cell+sim~variable))

temp[,hpv:=phie*(1-swe),]

result <- temp[,list(HPV=sum(hpv)),by=c("sim"),]

ggplot(result,aes(x=HPV))+geom_histogram()