匹配栅格在 foreach 循环中不起作用
Matching rasters do not work in foreach loop
我正在研究涵盖整个怀俄明州的栖息地占用预测。某些站点协变量栅格在预测中起作用，而其他具有匹配分辨率、范围等的栅格则不起作用。
下面是我的代码的一个简短的可重现示例。经过广泛的故障排除后，我发现我需要使用 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
回答
出现错误是因为您在 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中的NA，在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
因此，我会尝试弄清楚为什么您有 NA 开头。
我正在研究涵盖整个怀俄明州的栖息地占用预测。某些站点协变量栅格在预测中起作用，而其他具有匹配分辨率、范围等的栅格则不起作用。
下面是我的代码的一个简短的可重现示例。经过广泛的故障排除后，我发现我需要使用 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
回答
出现错误是因为您在 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中的NA，在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
因此，我会尝试弄清楚为什么您有 NA 开头。