更改栅格值的并行化

Question

尊敬的 stackoverflow 社区，

我阅读了一堆关于如何并行化栅格计算的答案。但是，我没有找到适合我的问题的答案。 我的问题是关于用数据框替换栅格中的值的并行化。这似乎很容易，我找到了一种方法来做到这一点。 但是，我有一个超过 500 万个像元的栅格和一个超过 340 000 行的数据框。因此，我的 for 循环需要很长时间。 因此，我想知道 foreach 循环来做我的计算。我的主要问题是最后结合我的结果。

我提供了一个简单的示例（来自：http://j.p.rossi.free.fr/rpackages/ecpaysage/TD_Analyse_quanti_paysage.html#(29)），以便您了解我想要做什么：

library(ECPaysage)
library(raster)
r <- raster(system.file("extdata/r.tif", package="ECPaysage")) 
plot(r)

这是我们需要一些补丁信息的土地利用栅格。为此，我们使用 SDMTools 包：

library(SDMTools)
res(r)

matbi <- r
matbi[] <- 0 # let's create the same raster with 0 values
w <- which(r[]==1) 
matbi[w] <- 1 # this creates a binary raster for the specific landuse == 1 

plot(matbi, axes=F, box=F, main="landuse 1 : build") 

我们要从这个二进制栅格中创建补丁信息。因此，我们创建了一个补丁栅格：

matpatch <- ConnCompLabel(matbi)
plot(matpatch, axes=F,box=F, main="patch ID landuse 1 :build") 

然后，我们从这个补丁栅格中提取补丁级别的信息：

patch <- PatchStat(mat=matpatch, cellsize = res(r)[1], latlon = FALSE)
names(patch)

dim(patch)

patch$patchID
dim(patch) 

现在，我们有了栅格中包含的补丁数量（每个像元都有一个补丁 ID）。 我们要在栅格中提取指定补丁 ID 的每个像元的补丁值。

shapeindex <- matpatch

w <- which(shapeindex[]==0) # cells that do not belong to any patch
shapeindex[w] <- NA

因此，要为每个补丁 ID 单元格提取相应的补丁索引，需要执行一个 for 循环：

system.time(for(p in 2:(dim(patch)[1])) {
  w <- which(shapeindex[]==patch$patchID[p])
  shapeindex[w] <- patch$shape.index[p]
})
plot(shape, axes=F,main="shape index - build", box=F, col=rev(terrain.colors(dim(patch)[1]-1))) 

这非常有效，但是对于我的数据（超过 500 万个单元格的栅格和超过 340 000 行的补丁数据框）需要太多时间（这不利于再现性）。 我的一个想法是用 foreach 并行化这个循环。但是，我的主要问题是最后合并我的栅格。如果我对数据框的行进行并行化，我将有 340 000 个矩阵要组合。它可能需要与 for 循环相同的时间。 因此，我想知道是否有一种方法（或多种方法，如函数）来加速处理和/或以简单的方式组合栅格。 感谢您提供的所有帮助。 最好的， 阿德里安娜

Answer 1

感谢您的回复@ulfelder 我按照你说的尝试了一些东西。我创建了一个函数，使用 lapply 并结合它们之间的栅格：

如果我在循环之前使用前面的代码，它会给出：

shape <- matpatch
w <- which(shape[]==0)
shape[w] <- NA
class(shape)
shape
plot(shape)

# create my function for one iteration of the for loop
RastVal <- function(monr,vec,i){
  require(raster)
  w <- which(monr[]==vec$patchID[i])
  x <- which(monr[]!=vec$patchID[i])
  monr[w] <- vec$shape.index[i]
  monr[x] <- 0
  return(monr)
}

# create a list of rasters containing the patch information
RastL <- lapply(X=c(2:(dim(patch)[1])),
                   FUN=function(x)RastVal(shape,patch,x))

# Combining the rasters
new.rast <- RastL[[1]]
for(i in 2:length(RastL)){
  x <- RastL[[i]][] != 0
  new.rast[x] <- RastL[[i]][x]
}
plot(new.rast) # should be the same as the raster shape in the previous code

但是，我不确定这是否更有效...也许，如果栅格变大，两种方式之间的时间差往往会减小...

system.time(for(p in 2:(dim(patch)[1])) {
  w <- which(shape[]==patch$patchID[p])
  shape[w] <- patch$shape.index[p]
})

system.time(RastL <- lapply(X=c(2:(dim(patch)[1])),
                   FUN=function(x)RastVal(shape,patch,x)))

new.rast <- RastL[[1]]
system.time(for(i in 2:length(RastL)){
  x <- RastL[[i]][] != 0
  new.rast[x] <- RastL[[i]][x]
})

此外，创建一个包含大约 340 000 个栅格的列表可能会导致内存问题... 再次感谢！

最好， 阿德里安娜