1

我在 R 中使用 gstat 包来生成顺序高斯模拟。我的电脑有 4 个内核,我尝试使用并行包按照Guzmán提供的脚本来并行化 krige() 函数,以回答问题如何在 R 中实现并行克里金法以加快进程?.

然而,由此产生的模拟与当时仅使用一个内核(无并行化)的模拟不同。它看起来是一个几何问题,但我不知道如何解决它。

接下来我将提供一个生成 2 个模拟的示例(使用 4 个内核)。您会看到,运行代码后,从并行化导出的模拟地图显示了一些伪影(如垂直线),并且与当时仅使用一个内核的地图不同。

代码需要库gstatspraster和。如果任何一行不起作用,请先运行。parallelspatstatlibrary()install.packages()

library(gstat)
library(sp)
library(raster)
library(parallel)
library(spatstat)

# 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')
gridded(srgr)<-~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)
r<-raster(sim)

# generate sample data (Poisson process)  
int<-0.02
rpp<-rpoispp(int,win=owin(c(0,nx),c(0,ny)))
df<-as.data.frame(rpp)
coordinates(df)<-~x+y 

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

# 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
no_cores<-detectCores()
cl<-makeCluster(no_cores)
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))
stopCluster(cl)

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

1 回答 1

0

我已经尝试了您的代码,我认为问题在于您拆分工作的方式:

parts <- split(x = 1:length(srgr), f = 1:no_cores)

在我的双核机器上,这意味着 srgr 中的所有奇数索引由一个进程处理,而所有偶数索引由另一个进程处理。这可能是您看到的垂直伪影的来源。

更好的方法应该是将数据拆分成连续的块,如下所示:

parts <- parallel::splitIndices(length(srgr), no_cores)

将此拆分与您的其余代码一起使用,我得到的结果看起来与顺序的结果相当。至少在我未经训练的眼睛里......


原始答案,这只是一个次要的影响。set.seed使用顺序和clusterSetRNGStream并行处理修复种子仍然可能有意义。

从我所读到的关于克里金的内容中,它要求您绘制随机数。这些随机数将与并行处理不同。vignette("parallel")有关详细信息,请参阅平行插图 ( ) 的第 6 部分。

于 2018-05-25T13:05:32.190 回答