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

Question

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

Answer 1

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

More info and help here: Parallelizing and clustering in R

# 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[[1]], parallelX[[2]])
mergeParallelX <- maptools::spRbind(mergeParallelX, parallelX[[3]])
mergeParallelX <- maptools::spRbind(mergeParallelX, parallelX[[4]])
mergeParallelX <- maptools::spRbind(mergeParallelX, parallelX[[5]])
mergeParallelX <- maptools::spRbind(mergeParallelX, parallelX[[6]])
mergeParallelX <- maptools::spRbind(mergeParallelX, parallelX[[7]])

# 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!

Answer 2

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.