【发布时间】:2020-07-12 11:52:57
【问题描述】:
我一直在尝试计算满足 R 特定条件的三个光栅文件的像素数(在 Windows 10 和 R 版本 3.6.3 中)。思路是利用这三张图片(经过预处理),统计下面代码if语句中满足条件的像素数量。
library(rgdal)
library(raster)
mydir<- file.path("D:/files")
setwd(mydir)
listado1 <- list.files(pattern="crNDVI*")
listado2 <- list.files(pattern="crNBR*")
listado3 <- list.files(pattern="dNBR*")
df <- data.frame(NULL)
for (i in 1:length(listado1)) {
r1 <- raster(listado1[i])
v1 <- values(r1)
for (j in 1:length(listado2)) {
r2 <- raster(listado2[j])
v2 <- values(r2)
for (k in 1:length(listado3)) {
r3 <- raster(listado3[k])
v3 <- values(r3)
if ((i==j) & (j==k)){
df2 <- data.frame(NULL)
df2 <- data.frame(v1, v2, v3)
burn <- subset(df2, df2$v1 >= 0.3 & df2$v1 <= 1.0 & df2$v2 >= 0.2 & df2$v2 <= 1.3 & df2$v3 > 0.1 & df2$v3 <= 1.3)
burned_area <- nrow(burn)*(xres(r1)*yres(r1))*10^(-4) # Hectares
name <- paste(substr(listado1[i], 8, 15), sep="")
df <- rbind(df, data.frame(name,burned_area))
print(paste("iteracion:",i,j,k," ok",sep = " "))
}
else {
print(paste("iteracion:",i,j,k," passed",sep = " "))
invisible()}
}
}
}
write.csv(df,file = "burned_area.csv")
可以观察到,我有三种类型的预处理光栅 (*.tif) 文件,它们以 crNDVI、crNBR 和 dNBR 开头,所有这些文件都列在每个“@987654327 @”变量。这个想法是迭代每种类型的文件并每天获取每个光栅文件的值。代码做了它应该做的(我得到了带有烧毁区域(满足条件的像素)的.csv文件,但是,它需要永远计算,因为每个光栅文件对应于特定的一天,因此结果是burned_area是每天的一个值。因此,发生这种情况的唯一可能性是i=j=k,但应用后一个代码,迭代会超过for循环中的每个可能性。代码从i = 1运行,然后j = 1 然后 k=1,2,3... 直到listado3[k] 的最后一个文件,当listado3[k] 的长度结束时,它传递到下面的j+1 索引并在k 中再次循环每次(不必要的)迭代。
谁能帮助我以更有效的方式做同样的事情?是否有可能为整个嵌套的for 循环强制“仅i=j=k”?非常感谢您的任何建议。在此先感谢豪尔赫。
注意:所有栅格文件的范围相同,nrow、ncol 和 ncell。
【问题讨论】:
-
嵌套 for 循环适用于每种组合。如果你想并行迭代,那就是
i = 1, j = 1, k = 1,然后是i =2, j = 2, k = 2,那么这只是一个循环。for(i in 1:length(listado1)) {j <- i; k <- i, ...}(或者,你甚至不需要j和k,只要在任何地方使用i。 -
哦,格雷戈尔,我刚读到这个......如果这是一个答案,我会给你支票。感谢您的热心帮助。我做到了(正如其他用户回答的那样)并且工作得很好。
标签: r gis raster rgdal r-modis