r - 裁剪和 mask 光栅堆栈的并行处理速度较慢

Question

我正在解决一个问题，我需要将光栅堆栈裁剪并屏蔽为一系列形状文件。这是一个可重现的示例:

# 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 循环来裁剪并将光栅堆栈屏蔽到形状文件:

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 

剪切不必要的循环，并在之前使用裁剪蒙版可以使速度提高约 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，您可以一步完成裁剪和 mask ，使用它可能会快 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 

现在，您可以再次查看并行化是否有帮助(如果仍然有必要的话)(对于此类问题来说，并行化很少值得)。也许对最外层循环执行此操作，以限制开销。但请注意，虽然您可以使用raster 做到这一点，但您不能使用terra 做到这一点(最终将在“幕后”可用)。

最后，您没有提供任何上下文，但您所做的事情看起来很奇怪。您想使用 extract 来代替吗？例如，您可以这样做:

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

这不是最快的方法，但也许是最清晰的代码。