4

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
3

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
2

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")

krigingPlot

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!

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.