I've been trying to calculate the number of pixels of three raster files that satisfy certain condition with R (in Windows 10 and R version 3.6.3). The idea is to use these three images (which are pre-processed) and count the amount of pixels that satisfy the condition in the if
statement of the following code.
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")
As it is possible to observe, I have three types of pre-processed raster (*.tif
) files that start with crNDVI
, crNBR
and dNBR
, all those files are listed in each "listado
" variables. The idea is to iterate over each type of files and obtain the values of each raster file per day. The code does what it should (I get the .csv
file with the burned area (pixels that satisfy the condition), however, it takes forever to compute because each raster file corresponds to an specific day and thus the result in burned_area
is a value per day. Therefore, the only possibility for that to happen is when i=j=k
, but applying the latter code the iteration is over each possibility in the for
loop. The code runs from i=1, then j=1 and then k=1,2,3... until the last file of listado3[k]
, and when the length of listado3[k]
finishes, it passes to the following j+1
index and looping all over again in k
to each (unnecesary) iteration.
Can anyone help me in order to do the same in a more efficient way? Is there any possibility to force an "only i=j=k
" for the whole nested for
loops? Id be very grateful for any advice. Thanks in advance, Jorge.
Note: all the raster files have the same extent, nrow
, ncol
and ncell
.
What is the need for the nested loops? It appears if you use one loop then i=j=k will always be true.