裁剪和屏蔽 rasterstack 的并行处理速度较慢

Question

我正在解决一个问题，我需要将栅格堆栈裁剪并遮盖到一系列 shapefile 中。这是一个可重现的例子：

# Example data ------------------------------------------------------------

#create example raster stack
r1 = raster(nrows=1000,ncol=1000,xmn=60,xmx=90,ymn=0,ymx=25)
rr = lapply(1:10, function(i) setValues(r1,runif(ncell(r1))))
rrstack=stack()
for (i in 1:length(rr)){
  stacknext=rr[[i]]
  rrstack=stack(rrstack,stacknext)
}


#create example shapefile list

lats=runif(10,min=0,max=25)
lons=runif(10,min=60,max=90)
exnames=paste0("city_",letters[1:10])
coords=data.frame(names=exnames,lats=lats,lons=lons)
coords_sf = st_as_sf(coords,coords=c("lons","lats"),crs=4326,dim ="XY")
circle1=st_buffer(coords_sf, 1E3)
circle100=st_buffer(coords_sf,1E5)
circle500=st_buffer(coords_sf,5E5)
circlist=list(circle1=circle1,circle100=circle100,circle500=circle500)

circlist_reproj=lapply(circlist,function(x) st_transform(x,crs(rrstack[[1]])))

我开发了一系列嵌套的 for 循环来裁剪栅格堆栈并将其屏蔽到 shapefile：

start <- proc.time()
citlist <- vector(mode='list',length=nrow(circlist_reproj[[1]]))
dellist <- vector(mode='list',length=length(circlist_reproj))
mystack <- stack()
for(k in 1:nrow(circlist_reproj[[1]])) {
  for(j in 1:length(circlist_reproj)) { 
    for (i in 1:nlayers(rrstack)){
      maskraster <- raster::mask(rrstack[[i]],circlist_reproj[[j]][k,])
      maskraster <- raster::crop(maskraster,circlist_reproj[[j]][k,])
      mystack <- stack(mystack,maskraster)
    }
    dellist[[j]] <- mystack
    mystack <- stack()
  }
  citlist[[k]] <- dellist
  dellist <- vector(mode='list',length=length(circlist_reproj))
}
basetime <- proc.time()-start

对于我的实际数据集，这需要几天才能完成。我开发了一个并行进程来加速上面的代码

cores <- detectCores()
cl <- makeCluster(cores[1]-2) #not to overload your computer
registerDoParallel(cl)
parstart=proc.time()
citlist <- vector(mode='list',length=nrow(circlist_reproj[[1]]))
dellist <- vector(mode='list',length=length(circlist_reproj))
for(k in 1:nrow(circlist_reproj[[1]])) {
  for(j in 1:length(circlist_reproj)) { 
    parrasterstack <- foreach(i=1:nlayers(rrstack),.packages=c('raster','sf')) %dopar% {
    maskraster <- raster::mask(rrstack[[i]],circlist_reproj[[j]][k,])
    raster::crop(maskraster,circlist_reproj[[j]][k,])
    }
  parrasterstack <- stack(parrasterstack)
  dellist[[j]] <- parrasterstack
  parrasterstack <- NULL
  }
citlist[[k]] <- dellist
dellist <- vector(mode='list',length=length(circlist_reproj))
}

stopCluster(cl)

运行两种处理方式，我发现顺序处理比并行处理快

> basetime #sequential process
   user  system elapsed 
 19.858   4.124  24.283 

> partime #parallel process
   user  system elapsed 
 41.415  10.153 132.262

我在另一台具有更多内核的机器上重复了代码，得到了更大的差异。为什么并行处理各方面都比较慢，有没有办法加快处理速度？

Answer 1

通过你的顺序方法，我得到了

proc.time()-start
#   user  system elapsed 
#  18.86    2.57   21.44

剪掉一个不必要的循环，并在 before 掩码下使用裁剪使速度提高大约 8 倍

start <- proc.time()
citlist <- vector(mode='list',length=nrow(circlist_reproj[[1]]))

for(k in 1:nrow(circlist_reproj[[1]])) {
    for(j in 1:length(circlist_reproj)) { 
        r <- crop(rrstack, circlist_reproj[[j]][k,])
        citlist[[k]][[j]] <- mask(r, circlist_reproj[[j]][k,])
    }
}
proc.time()-start
#   user  system elapsed 
#   2.41    0.30    2.70

并且 terra 比 raster

快 3 倍

library(terra)
rrstack <- rast(rrstack)
circlist_reproj <- lapply(circlist_reproj, vect)

start <- proc.time()
citlist <- vector(mode='list',length=nrow(circlist_reproj[[1]]))

for(k in 1:nrow(circlist_reproj[[1]])) {
    for(j in 1:length(circlist_reproj)) { 
        r <- crop(rrstack, circlist_reproj[[j]][k,])
        citlist[[k]][[j]] <- mask(r, circlist_reproj[[j]][k,])
    }
}
proc.time()-start
#   user  system elapsed 
#   0.70    0.19    0.89

所以速度提高了 24 倍。

另外，为了安全起见，循环真的应该这样写：

start <- proc.time()
citlist <- vector(mode='list',length=length(circlist_reproj))
for(i in 1:length(circlist_reproj)) { 
    for(j in 1:nrow(circlist_reproj[[i]])) {
        r <- crop(rrstack, circlist_reproj[[i]][j,])
        citlist[[i]][[j]] <- mask(r, circlist_reproj[[i]][j,])
    }
}
proc.time()-start

使用 terra，您可以一步完成裁剪和蒙版，使用它可能会快 35 倍。

start <- proc.time()
citlist <- vector(mode='list',length=length(circlist_reproj))
for(i in 1:length(circlist_reproj)) { 
    for(j in 1:nrow(circlist_reproj[[i]])) {
        citlist[[i]][[j]] <- crop(rrstack, circlist_reproj[[i]][j,], mask=TRUE)
    }
}
proc.time()-start
# user  system elapsed 
# 0.46    0.11    0.56

您现在可以再次查看并行化是否有帮助（如果并行化仍然有必要）（对于此类问题很少值得这样做）。也许对 outer-most 循环这样做，以限制开销。但请注意，虽然您可以使用 raster 做到这一点，但您不能使用 terra 做到这一点（最终将在“引擎盖下”提供）。

最后，你没有提供上下文，但你所做的看起来很奇怪。您想改用 extract 吗？你可以这样做，例如：

e <- lapply(circlist_reproj, function(v) extract(rrstack, v))

这不是最快的方法，但可能是最清晰的代码。

裁剪和屏蔽 rasterstack 的并行处理速度较慢

parallel processing slower for cropping and masking rasterstack

parallel-processing

foreach

r

raster