多个嵌套 for 循环仅在 R 中的 i=j=k 索引下迭代
Multiple nested for loops to iterate only under i=j=k indexes in R
我一直在尝试用 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
开头,所有这些文件都列在每个“listado
”变量中。这个想法是迭代每种类型的文件并每天获取每个光栅文件的值。代码做了它应该做的事情(我得到了 .csv
文件,其中包含燃烧区域(满足条件的像素),但是,它需要永远计算,因为每个光栅文件对应于特定的一天,因此结果 burned_area
是每天的值。因此,发生这种情况的唯一可能性是 i=j=k
,但应用后一种代码,迭代遍历了 for
循环中的每种可能性。代码运行从i=1,然后j=1,再到k=1,2,3...,直到listado3[k]
的最后一个文件,当listado3[k]
的长度结束后,传递给下面的j+1
索引并在 k
中再次循环到每个(不必要的)迭代。
任何人都可以帮助我以更有效的方式做同样的事情吗?是否有可能为整个嵌套 for
循环强制使用 "only i=j=k
"?我非常感谢任何建议。提前致谢,豪尔赫。
注意:所有光栅文件的范围相同,nrow
、ncol
和ncell
。
嵌套循环需要什么?看起来如果你使用一个循环那么 i=j=k 将永远是真的。
这是一个最小的独立示例。
library(raster)
set.seed(0)
r1 <- raster(nrow=10, ncol=10, vals=1:100) / 100
r2 <- raster(nrow=10, ncol=10, vals=c(1:50, 50:1)) / 100
r3 <- raster(nrow=10, ncol=10, vals=c(rep(1, 10), 11:100)) / 100
x <- ((r1==r2) & (r2==r3))
s <- stack(r1, r2, r3)
s <- mask(s, x, maskvalue=1)
我一直在尝试用 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
开头,所有这些文件都列在每个“listado
”变量中。这个想法是迭代每种类型的文件并每天获取每个光栅文件的值。代码做了它应该做的事情(我得到了 .csv
文件,其中包含燃烧区域(满足条件的像素),但是,它需要永远计算,因为每个光栅文件对应于特定的一天,因此结果 burned_area
是每天的值。因此,发生这种情况的唯一可能性是 i=j=k
,但应用后一种代码,迭代遍历了 for
循环中的每种可能性。代码运行从i=1,然后j=1,再到k=1,2,3...,直到listado3[k]
的最后一个文件,当listado3[k]
的长度结束后,传递给下面的j+1
索引并在 k
中再次循环到每个(不必要的)迭代。
任何人都可以帮助我以更有效的方式做同样的事情吗?是否有可能为整个嵌套 for
循环强制使用 "only i=j=k
"?我非常感谢任何建议。提前致谢,豪尔赫。
注意:所有光栅文件的范围相同,nrow
、ncol
和ncell
。
嵌套循环需要什么?看起来如果你使用一个循环那么 i=j=k 将永远是真的。
这是一个最小的独立示例。
library(raster)
set.seed(0)
r1 <- raster(nrow=10, ncol=10, vals=1:100) / 100
r2 <- raster(nrow=10, ncol=10, vals=c(1:50, 50:1)) / 100
r3 <- raster(nrow=10, ncol=10, vals=c(rep(1, 10), 11:100)) / 100
x <- ((r1==r2) & (r2==r3))
s <- stack(r1, r2, r3)
s <- mask(s, x, maskvalue=1)