1

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.

Jorge V
  • 73
  • 7
  • 1
    Nested for loops are for every combination. If you want to iterate in parallel, that is `i = 1, j = 1, k = 1`, then `i =2, j = 2, k = 2`, then that's only one loop. `for(i in 1:length(listado1)) {j <- i; k <- i, ...}` (Or, you don't even need `j` and `k`, just use `i` everywhere. – Gregor Thomas Mar 31 '20 at 20:05
  • Oh Gregor, I just read this... if this would have been an answer I would have given you the check. Thank you for your kind help. I did it (as other user answered) and worked perfectly. – Jorge V Apr 01 '20 at 01:24

2 Answers2

1

What is the need for the nested loops? It appears if you use one loop then i=j=k will always be true.

DPek
  • 180
  • 2
  • 15
0

Here is a minimal self-contained example.

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)
Nimantha
  • 6,405
  • 6
  • 28
  • 69
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • Hey Robert. Thank you for your advice. Nevertheless, I do not find it necessary because I believe that would make the question more complex. The question is related to nested loops in R, not regarding to raster analysis. That is why I posted in this forum and not in GIS. – Jorge V Apr 01 '20 at 01:03
  • Then it is even more important that you rework your question. In that case you should not include raster data and create some simple example data. It is much to difficult to understand what you want to do or what your question really is. – Robert Hijmans Apr 01 '20 at 17:26
  • Hey Robert, thanks again. My problem was solved and thus it was understood by those who tried to help. Have a nice week. – Jorge V Apr 01 '20 at 17:30