I am using gstat package in R to generate sequential gaussian simulations. My pc have 4 cores and I tried to parallelize the krige() function using the parallel package following the script provided by Guzmán to answer the question How to achieve parallel Kriging in R to speed up the process?.

The resulting simulations are, however, wrong, it looks a geometry problem, but i can't find how to tackle the problem. Next i will provide an example (4 cores). You will see that after running the code, the maps from parallelization have vertical lines, and are different from the ones using only one core at the time.


# create a regular grid
nx=100 # number of columns
ny=100 # number of rows
srgr <- expand.grid(1:ny, nx:1)
names(srgr) <- c('x','y')

# generate a spatial process (unconditional simulation)
g<-gstat(formula=z~x+y, locations=~x+y, dummy=T, beta=15, model=vgm(psill=3, range=10, nugget=0,model='Exp'), nmax=20)
sim <- predict(g, newdata=srgr, nsim=1)

# generate sample data (Poisson process)  

# assign raster values to sample data
dfpp <-raster::extract(r,df,df=TRUE)
smp<-smp[complete.cases(smp), ]

# fit variogram to sample data
vs <- variogram(sim1~1, data=smp)
m <- fit.variogram(vs, vgm("Exp"))
plot(vs, model = m)

# generate 2 conditional simulations with one core processor
one <- krige(formula = sim1~1, locations = smp, newdata = srgr, model = m,nmax=12,nsim=2)

# plot simulation 1 and 2: statistics (min, max) are ok, simulations are also ok.
spplot(one["sim1"], main = "conditional simulation")
spplot(one["sim2"], main = "conditional simulation")

# generate 2 conditional with parallel processing
parts <- split(x = 1:length(srgr), f = 1:no_cores)
clusterExport(cl = cl, varlist = c("smp", "srgr", "parts","m"), envir = .GlobalEnv)
clusterEvalQ(cl = cl, expr = c(library('sp'), library('gstat')))
par <- parLapply(cl = cl, X = 1:no_cores, fun = function(x) krige(formula=sim1~1, locations=smp, model=m, newdata=srgr[parts[[x]],],  nmax=12, nsim=2))

# merge all parts    
mergep <- maptools::spRbind(par[[1]], par[[2]])
mergep <- maptools::spRbind(mergep, par[[3]])
mergep <- maptools::spRbind(mergep, par[[4]])

# create SpatialPixelsDataFrame from mergep
mergep <- SpatialPixelsDataFrame(points = mergep, data = mergep@data)

# plot mergep: statistics (min, max) are ok, but simulated maps show "vertical lines". i don't understand why.
spplot(mergep[1], main = "conditional simulation")
spplot(mergep[2], main = "conditional simulation")
  • It's not really a geometry problem. I think the main reason is that using nsim in the krige function will follow a single random path through the data. So running it 4 times using each core will follow four different random paths. The result will have a vertical lines pattern (that is due how the splitting was done previously, try split(x = 1:length(srgr), f = cut(1:length(srgr), no_cores)) and you will notice horizontal lines following the splitted blocks). I tried using set.seed(somevalue) to have a reproducible example but each time krige runs follow a different path. – Guzmán May 27 '18 at 23:34
  • Thank you @Guzmán for your comments. I only could not understand your last sentence ("I tried using...path."). – maluicr May 28 '18 at 16:48
  • I mean, a random process take place when using nsim. So I wanted to generate the same random process for each process in each core. I thought that if I set the same seed e.g. set.seed(43) before running the krige function I would keep the same random process and achieve the correct result. But it didn't work. – Guzmán May 29 '18 at 14:00

Your Answer

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

Browse other questions tagged or ask your own question.