微信公众号搜"智元新知"关注
微信扫一扫可直接关注哦!

匹配栅格在 foreach 循环中不起作用 我最初为堆叠对齐栅格所采取的步骤 - 仅供参考从这里开始重现错误

如何解决匹配栅格在 foreach 循环中不起作用 我最初为堆叠对齐栅格所采取的步骤 - 仅供参考从这里开始重现错误

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

我的代码一个可重现的简短示例如下。经过大量故障排除后,我发现我需要使用 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),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,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,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,tmp2) # works

基本原理

事实证明,您的坏栅格缺少值:

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

如果我们人为地去除NA中的sNoJoy,并在NA中随机写入1个sampleJoy,那么状态翻转:

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

### run same code as above

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

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

版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。