Raster calculation loop - R

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 []: blue band image 1

Position []: green band image 1

Position []: red band image 1

Position []: NIR band image 1

Position []: blue band image 2

...

Position []: 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[] - images[]) / (images[] + images[])

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

and

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

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

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?

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[],ndwi_list[],ndvi_list[],ndwi_list[])) • Hi @aldo_tapia, thanks for your reply. I was struggling with the trick to iterate properly. My approach was to define a function and then use lapply but I guess this would end up in the same result. What do you mean by saying that position [] and [] have to be the same? – GCGM Oct 31 '17 at 13:35
• Because I copied lsat 27 times. The loop works because the first and last item in ndvi_list and ndwi_list are good calculated... As are copies of 1 scene, the results, in this case, are the same... It's only a test to try loop effectiveness – aldo_tapia Oct 31 '17 at 13:39