在 foreach 错误处理栅格堆栈图后

编程问答 2022-05-26

问题描述

我正在通过未标记的预测函数处理大型数据堆栈 (pm)。处理后，我使用 paroutPred、SE、Lower 和 Upper 作为模板来粘贴 Predicted、SE、lower 和 Upper 数据并堆叠这些栅格以进行绘图。

我当前的代码如下。 foreach 循环似乎运行良好，我正在获取必要的变量。毕竟说了又做了， compareRaster(rs,pm) 结果为TRUE。只是数值不同。

执行 plot(rs) 后，四个光栅绘制，但它们都被提升，如下所示：What is drawn should be occupancy probability for a species across a map of Wyoming。

我还没弄清楚出了什么问题。我在 for 循环中按顺序运行（处理时间为 3 天）以确认栅格输出是错误的。

有人对我的问题有任何见解吗？非常感谢所有帮助。

paroutPred<- pm[[1]]
paroutSE<- pm[[1]]
paroutLower<- pm[[1]]
paroutUpper<- pm[[1]]
paroutPred[]<- NA
paroutSE[]<- NA
paroutLower[]<- NA
paroutUpper[]<- NA

library(doSNow)

nc<- detectCores()-1
cl<- makeCluster(nc);cl
registerDoSNow(cl)

comb<- function(...){
  mapply('rbind',...,SIMPLIFY = F)
}
predictions<- 
  foreach(i = 1:nrow(pm),.combine = 'comb',.multicombine = T,.maxcombine = 200,.packages = c("unmarked","raster"),.verbose = T
          )%dopar%{ 
            test<- cellFromrow(pm,row=i)
            # make into a data.frame for prediction
            tmp<- data.frame(pm[test])
            # test which are na
            na<- sapply(1:nrow(tmp),FUN = function(x){any(is.na(tmp[x,]))})
            # deal with writing the data back to raster
            if(length(which(na)) != nrow(tmp)){
              # Predict the new data
              pred<- predict(fmBest,"state",tmp)
              }
          list(test,na,pred)  
         }
stopCluster(cl)

predlist<- list(predictions[[3]][["Predicted"]],predictions[[3]][["SE"]],predictions[[3]][["upper"]],predictions[[3]][["lower"]])
pred<- do.call(cbind,lapply(predlist,data.frame))
names(pred)<- c("Predicted","SE","upper","lower")
test<- predictions[[1]]
na<- data.frame(predictions[[2]])

paroutPred[test[!na]]<- pred$Predicted
names(paroutPred)<- "Predicted"
# Save prediction
paroutSE[test[!na]]<- pred$SE
names(paroutSE)<- "SE"
# Save prediction
paroutLower[test[!na]]<- pred$lower
names(paroutLower)<- "lower"
# Save prediction
paroutUpper[test[!na]]<- pred$upper
names(paroutUpper)<- "upper"

writeraster(paroutPred,"Ppred_PEFA.tif",format = "GTiff",overwrite = TRUE)
writeraster(paroutSE,"Pse_PEFA.tif",overwrite = TRUE)
writeraster(paroutLower,"Plower_PEFA.tif",overwrite = TRUE)
writeraster(paroutUpper,"Pupper_PEFA.tif",overwrite = TRUE)

Ppred.PEFA <- raster(paste(getwd(),"/Ppred_PEFA.tif",sep=""))
Pse.PEFA <- raster(paste(getwd(),"/Pse_PEFA.tif",sep=""))
Plower.PEFA <- raster(paste(getwd(),"/Plower_PEFA.tif",sep=""))
Pupper.PEFA <- raster(paste(getwd(),"/Pupper_PEFA.tif",sep=""))

rs<- stack(c("Ppred_PEFA.tif","Pupper_PEFA.tif"))
plot(rs)

解决方法

我终于搞定了。我使用了下面的代码。原来我一直在努力做事情，导出变量然后尝试将它们操作成光栅形式。一旦我合并了光栅模板并使用快速光栅函数将它们作为输出，BOOM！

从 foreach 循环中正确绘制的栅格。

唯一的捕获就是这样做会占用大量内存。由于每个栅格模板的大小为 1.2Gb，并且所有四个栅格模板都导出到每个集群，因此加起来很快。在具有 128 Gb RAM 的机器上，我将内存最大化以在 Windows 的 9 个内核上运行。

# Loop through the rows and columns,subset the data,make a prediction
paroutPred<- pm[[1]]
paroutSE<- pm[[1]]
paroutLower<- pm[[1]]
paroutUpper<- pm[[1]]
paroutPred[]<- NA
paroutSE[]<- NA
paroutLower[]<- NA
paroutUpper[]<- NA

library(doSNOW)
#Assemble cluster for parallel processing
nc<- detectCores()-8
cl<- makeCluster(nc);cl
registerDoSNOW(cl)
#Combine function to use in foreach loop w/multiple returns
comb<- function(...){
  mapply('rbind',...,SIMPLIFY = F)
}
#foreach loop for predictions
system.time(
predictions<- 
  foreach(i = 1:nrow(pm),.combine = 'comb',.multicombine = T,.maxcombine = 200,.packages = c("unmarked","raster"),.verbose = T
          )%dopar%{ 
            
            test<- cellFromRow(pm,i)
            
            # make into a data.frame for prediction
            tmp<- data.frame(pm[test])
            
            # test which are na
            na<- sapply(1:nrow(tmp),FUN = function(i){any(is.na(tmp[i,]))})
            
            # deal with writing the data back to raster
            if(length(which(na)) != nrow(tmp)){
              #   # Predict the new data
              pred<- predict(fmBest,"state",tmp)
            }
 #List of raster format non-raster objects to return  
      list(       
      paroutPred[test[!na]]<- pred$Predicted,paroutSE[test[!na]]<- pred$SE,paroutLower[test[!na]]<- pred$lower,paroutUpper[test[!na]]<- pred$upper)
          })
#Close the cluster
stopCluster(cl)
#Return data to raster form
paroutPred<- raster(predictions[[1]],template = paroutPred)
paroutSE<- raster(predictions[[2]],template = paroutSE)
paroutLower<- raster(predictions[[3]],template = paroutLower)
paroutUpper<- raster(predictions[[4]],template = paroutUpper)

找到了更好的方法。 RAM使用少得多。取出模板并意识到我可以在循环后处理光栅格式（废话）。

     #predict function of foreach loop here
}
#List of non-raster prediction results to return  
      list(pred$Predicted,pred$SE,pred$lower,pred$upper)
          }
#Return data to raster form
paroutPred<- raster(predictions[[1]],template = paroutUpper)
parallel.foreachrrr-raster