# How to achieve parallel Kriging in R to speed up the process?

I have tried the `gstat` package in `R` to perform bathymetry predictions based on Kriging interpolation methods. However, the process is very time consuming due the large amount of sampling data I'm using. Also, the whole process in `R` is running only in one core of the four cores my PC have.

1. I wonder if it is possible to parallelize the Kriging (`krige`) function of the `gstat` package or
2. Any `R` package to perform parallel Kriging

Once you've got a variogram fitted then kriging is trivially parallelizable. Kriging predictions are independent of each other.

So, divide your prediction points (grid) into N sets, where N is your number of cores, and do the predictions for each of the point sets on a separate core. Merge the predictions afterwards.

You can use any of the paralleling tricks for R, the appropriate one depends on your architecture. Got a cluster? Spread over the cluster. Read the Parallel Processing Task View for more precise details.

I don't know of any extras to gstat that do this already, but I assume you've done your research and at least scanned CRAN and not found anything.

• Nice! I didn't realize that dividing the prediction points grid would to the trick! Thanks! I will try the `paralell` package! – Guzmán Apr 22 '17 at 13:28

This is a simple and reproducible example with `meuse` dataset of `gstat` package based on @Spacedman answer and using the `parallel` package:

``````# Load libraries
library('sp')
library('gstat')
library('parallel')

# Load example data: meuse dataset
data("meuse") # data
data("meuse.grid") # 40m x 40m prediction grid for meuse dataset

# Create SpatialPointsDataFrame from meuse
coordinates(meuse) = ~x+y

# Create SpatialPixelsDataFrame from meuse.grid
gridded(meuse.grid) = ~x+y

# Create Spheric Variogram model with predefined parameters (psill, range and nugget)
m <- vgm(psill = 0.59, model = "Sph", range = 874, nugget = 0.04)

# Ordinary kriging (one core processor)
system.time(x <- krige(formula = log(zinc)~1, locations = meuse, newdata = meuse.grid, model = m))

# Plot x
spplot(x["var1.pred"], main = "ordinary kriging predictions")

# Ordinary kriging (parallel processing)

# Calculate the number of cores
no_cores <- detectCores() - 1

# Initiate cluster (after loading all the necessary object to R environment: meuse, meuse.grid, m)
cl <- makeCluster(no_cores)

parts <- split(x = 1:length(meuse.grid), f = 1:no_cores)

clusterExport(cl = cl, varlist = c("meuse", "meuse.grid", "m", "parts"), envir = .GlobalEnv)
clusterEvalQ(cl = cl, expr = c(library('sp'), library('gstat')))

system.time(parallelX <- parLapply(cl = cl, X = 1:no_cores, fun = function(x) krige(formula = log(zinc)~1, locations = meuse, newdata = meuse.grid[parts[[x]],], model = m)))

stopCluster(cl)

# Merge all the predictions
mergeParallelX <- maptools::spRbind(parallelX[], parallelX[])
mergeParallelX <- maptools::spRbind(mergeParallelX, parallelX[])
mergeParallelX <- maptools::spRbind(mergeParallelX, parallelX[])
mergeParallelX <- maptools::spRbind(mergeParallelX, parallelX[])
mergeParallelX <- maptools::spRbind(mergeParallelX, parallelX[])
mergeParallelX <- maptools::spRbind(mergeParallelX, parallelX[])

# Create SpatialPixelsDataFrame from mergeParallelX
mergeParallelX <- SpatialPixelsDataFrame(points = mergeParallelX, data = mergeParallelX@data)

# Plot mergeParallelX
spplot(mergeParallelX["var1.pred"], main = "ordinary kriging predictions")
`````` Note:

• I tested the code above with my own data and in a total of 6364 data locations, 3429 newdata prediction positions and 7 core processors the amount of time of the parallel process of Kriging was half the non-parallel process!