在 foreach 处理栅格堆栈图不正确之后

After foreach processing raster stack plots incorrectly

我正在通过未标记的预测函数处理大型数据堆栈 (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", format = "GTiff", overwrite = TRUE)
writeRaster(paroutLower, "Plower_PEFA.tif", format = "GTiff", overwrite = TRUE)
writeRaster(paroutUpper, "Pupper_PEFA.tif", format = "GTiff", 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", "Pse_PEFA.tif", "Plower_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 使用少得多。取出模板并意识到我可以在循环后处理光栅格式(duh)。

     #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 = paroutPred)
paroutSE<- raster(predictions[[2]], template = paroutSE)
paroutLower<- raster(predictions[[3]], template = paroutLower)
paroutUpper<- raster(predictions[[4]], template = paroutUpper)