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

Question

我正在通过未标记的预测函数处理大型数据堆栈 (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)

Answer 1

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

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

只有 catch 这样做会占用大量 RAM。每个栅格模板的大小为 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)

找到了更好的方法。更少的内存使用。取出模板并意识到我可以在循环之后处理光栅格式（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)