匹配栅格在 foreach 循环中不起作用

Question

我正在进行一项涵盖整个怀俄明州的栖息地占用率预测。某些站点协变量栅格在预测中起作用，而其他具有匹配分辨率、范围等的栅格则不起作用。

我的代码的一个简短的可复制示例如下。经过广泛的故障排除后，我发现我需要使用 5 个栅格中的 3 个栅格，导致该脚本失败，所有栅格都出现相同的错误。 我假设我的栅格以某种方式损坏（？）但想看看是否有人对可能发生的事情有其他想法。

数据位于this link。数据是未标记的对象（保存为 .rds）和 2 个非常小的剪辑：1. 有效的栅格，以及 2. 无效的栅格之一

我最初对齐栅格以进行堆叠的步骤 - 仅供参考

#ndvi <- raster(paste(getwd(), "./Original_rasters/ndvi_summer.TIF", sep = ""))

#precip <- raster(paste(getwd(), "./Original_rasters/bioclim15.TIF", sep = ""))

#temp <- projectExtent(original, original)
#res(temp) <- 220

#sNoJoy <- resample(ndvi, temp)
#sampleJoy <- resample(precip, temp)

从这里开始重现错误

library(raster)
library(unmarked)
sampleData <- readRDS("./sampleData.RDS")

# Formula referencing raster that does not work
fmTest <- occu(~ Day + Min_TempC + AvailTN_prop ~ ndvi.summer, 
               sampleData)
fmTest

# Formula referencing raster that DOES work
#fmTest <- occu(~ Day + Min_TempC + AvailTN_prop ~ bc15_220, 
#               sampleData)

# Formula for both variables that fails also
#fmTest <- occu(~ Day + Min_TempC + AvailTN_prop ~ ndvi.summer + bc15_220,
#               sampleData)

# Load and name rasters so formulas can find them
sampleJoy <- raster(paste(getwd(), "./sampleGoodRas.TIF", sep = ""))
names(sampleJoy) <- "bc15_220"
sNoJoy <- raster(paste(getwd(), "./sampleBadRas.TIF", sep = ""))
names(sNoJoy) <- "ndvi.summer"

compareRaster(sampleJoy, sNoJoy) # Returns "[1] TRUE"
 
# pm <- stack(sampleJoy)
pm <- stack(sNoJoy)
# pm <- stack(sNoJoy, sampleJoy)

###########################Combine Function#####################################

comb <- function(x, ...) {
  mapply("rbind", x, ..., SIMPLIFY = F)
}

# This combine function allows foreach to return a list containing multiple 
# matrices making it easy to insert results into raster templates

############################ Foreach Code #####################################
#Assemble cluster for parallel processing.  Code works in Windows or other O/S#

ifelse(Sys.info()["sysname"] != "Windows", 
       c(require(doMC), nc <- detectCores()-1, registerDoMC(nc)),
       c(require(doParallel), nc <- detectCores()-1, cl <- makeCluster(nc), 
         registerDoParallel(cl)))

# Foreach loop returning predicted values, SE, LCI, and UCI
pred <- foreach(i = 1:nrow(pm), .combine = comb, .multicombine = T,
                .maxcombine = 90,
                .packages = c("unmarked", "raster")) %dopar% {
  
  # make raster into a data.frame row by row for prediction
  tmp <- as.data.frame(pm[i,], xy = T)
  
  # Predict the new data
  pred <- predict(fmTest, "state", tmp)
  
  # Make a list of 4 matrices to retrieve them from the loop      
  list(Predicted = pred$Predicted,
       SE = pred$SE,
       lower = pred$lower,
       upper = pred$upper)
}

# Close the cluster
stopCluster(cl)

## Using sampleJoy produces a list of 4 matrices which are easily coerced into
## raster format: Prediction, SE, lower, and upper, as it should.

## Using sNoJoy produces:  
# Error in { : 
# task 1 failed - "Matrices must have same number of rows in 
# cbind2(.Call(dense_to_Csparse, x), y)" 

## Rasters are the same extent, same origin, same resolution, etc.
# 

### Not Working
# > pm
# class      : RasterStack 
# dimensions : 15, 2675, 40125, 1  (nrow, ncol, ncell, nlayers)
# resolution : 220, 220  (x, y)
# extent     : 201539.7, 790039.7, 647050.2, 650350.2  (xmin, xmax, ymin, ymax)
# crs        : +proj=lcc +lat_0=41 +lon_0=-107.5 +lat_1=41 +lat_2=45 +x_0=500000 +y_0=200000 +datum=NAD83 +units=m +no_defs 
# names      : ndvi.summer 
# min values :  0.09507491 
# max values :   0.8002191 
# 

### Working Raster
# > pm
# class      : RasterStack 
# dimensions : 15, 2675, 40125, 1  (nrow, ncol, ncell, nlayers)
# resolution : 220, 220  (x, y)
# extent     : 201539.7, 790039.7, 647050.2, 650350.2  (xmin, xmax, ymin, ymax)
# crs        : +proj=lcc +lat_0=41 +lon_0=-107.5 +lat_1=41 +lat_2=45 +x_0=500000 +y_0=200000 +datum=NAD83 +units=m +no_defs 
# names      : bc15_220 
# min values :       14 
# max values :       66

Answer 1

回答

出现错误是因为您在sNoJoy 中有缺失。如果那些没有丢失，它会工作得很好。

问题重写

您的问题与您的并行代码无关。归结为：

fmTest <- occu(~ Day + Min_TempC + AvailTN_prop ~ ndvi.summer, sampleData)
fmTest2 <- occu(~ Day + Min_TempC + AvailTN_prop ~ bc15_220, sampleData)

pm <- stack(sNoJoy)
pm2 <- stack(sampleJoy)

tmp <- as.data.frame(pm[1,], xy = T)
tmp2 <- as.data.frame(pm2[1,], xy = T)

pred <- predict(fmTest, "state", tmp) # fails
pred2 <- predict(fmTest2, "state", tmp2) # works

基本原理

事实证明，您的错误栅格缺少值：

table(is.na(sNoJoy[]))
#FALSE  TRUE 
#35998  4127 

如果我们人为去掉sNoJoy中的NAs，在sampleJoy中随机写入1个NA，那么状态翻转：

sNoJoy[is.na(sNoJoy[])] <- 1
sampleJoy[10] <- NA

### run same code as above

pred <- predict(fmTest, "state", tmp) # now works
pred2 <- predict(fmTest2, "state", tmp2) # now fails

因此，我会尝试弄清楚您为什么要从 NAs 开始。