Raster calculation loop - R

Question

I have a raster stack of 108 elements

The stack 'images' is created using 4 bands (RGB+NIR) of 27 satellite images. The sequence of the rasters in the stack is:

Position [[1]]: blue band image 1

Position [[2]]: green band image 1

Position [[3]]: red band image 1

Position [[4]]: NIR band image 1

Position [[5]]: blue band image 2 

...

Position [[108]]: NIR band image 27

I want to derive two basic vegetation indices for each image (27 NDVI and 27 NDWI2 indices): NDVI - (NIR-R)/(NIR+R) and NDWI2 (G-NIR)/(G+NIR)

So far I have doing it manually. For example:

NDVI_1-> (images[[4]] - images[[3]]) / (images[[4]] + images[[3]])

NDVI_2-> (images[[8]] - images[[7]]) / (images[[8]] + images[[7]])

and

NDWI2_1-> (images[[2]] - images[[4]]) / (images[[2]] + images[[4]])


NDWI2_2-> (images[[6]] - images[[8]]) / (images[[6]] + images[[8]])

Since my raster stack as a pattern, I am looking for an automatic way to derive the indices. Something like a loop that calculates the indeces, save it with the appropriate name and then moves to the next 'position' to repeat the process.

Any help on that?

Answer 1

Just use a for loop with a little trick to iterate.

First, as always is recommended, a reproducible example:

library(raster)
data('lsat',package='RStoolbox')

names(lsat)[1:4] <- c('blue','green','red','NIR')

images <- list()

for(i in 1:27){
  images[[i]] <- lsat[[1:4]]
}

images <- stack(images)

images
## class       : RasterStack 
## dimensions  : 310, 287, 88970, 108  (nrow, ncol, ncell, nlayers)
## resolution  : 30, 30  (x, y)
## extent      : 619395, 628005, -419505, -410205  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=22 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## names       : blue.1, green.1, red.1, NIR.1, blue.2, green.2, red.2, NIR.2, blue.3, green.3, red.3, NIR.3, blue.4, green.4, red.4, ... 
## min values  :     54,      18,    11,     4,     54,      18,    11,     4,     54,      18,    11,     4,     54,      18,    11, ... 
## max values  :    185,      87,    92,   127,    185,      87,    92,   127,    185,      87,    92,   127,    185,      87,    92, ... 

Create two lists to store NDVI and NDWI indices:

ndvi_list <- list()
ndwi_list <- list()

Loop through images with a condition to skip every 4 layers per iteration:

for(i in 1:27){
  ndvi_list[[i]] <- overlay(images[[((i-1)*4+3)]],images[[((i-1)*4+4)]],fun=function(x,y) (y-x)/(y+x))
  ndwi_list[[i]] <- overlay(images[[((i-1)*4+4)]],images[[((i-1)*4+2)]],fun=function(x,y){(y-x)/(y+x)})
}

Check results (need to be the same, I used one scene in this example):

plot(stack(ndvi_list[[1]],ndwi_list[[1]],ndvi_list[[27]],ndwi_list[[27]]))