3D Kriging with Conditional Gaussian Simulation in R

3d_kriging_cgs.R

library(gstat)
library(lattice)

# 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 empty grid to krige
range1D <- seq(from = 0, to = 1, length = 20)
grid3D <- expand.grid(x = range1D, y = range1D, z = range1D)
gridded(grid3D) = ~x+y+z

# Perform CGS with 10 realizations; maxdist & nmax important for speed of calculation.
res3D <- krige(formula = v ~ 1, data3D, grid3D, model = vgm(1, "Exp", .2),nsim=10,maxdist=10,nmax=9)

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

Modified version 2025: add sp() package, show ordinary kriging, finer x/y grid

library(gstat)
library(lattice)
library(sp)

# 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 empty grid to krige
range1D <- seq(from = 0, to = 1, length = 51)
rangez <- seq(from = 0, to = 1, length = 25)
grid3D <- expand.grid(x = range1D, y = range1D, z = rangez)
gridded(grid3D) = ~x+y+z

# Perform CGS with 10 realizations; maxdist & nmax important for speed of calculation.
res3D <- krige(formula = v ~ 1, data3D, grid3D,
               model = vgm(1, "Exp", .2),nsim=10,maxdist=10,nmax=9)

## ordinary kriging
res3Dx <- krige(formula = v ~ 1, data3D, grid3D,
               model = vgm(1, "Exp", .2),maxdist=10,nmax=9)

# Plot Results
levelplot(var1.pred ~ x + y | z, as.data.frame(res3Dx),
          useRaster = TRUE)