how to loop ndvi using R in sentinel-2

Question

I have a two different stacks: one for band 4 (RED) and other with band 8 (NIR). I'm trying to loop the calculation of NDVI in an area of interest using R.

Well, this is the loop that I'm trying to use: First, I'm trying to crop the stacks using a shapefile, then calculate the NDVI of the area of interest and save this to the disk.

for(i in 1){

  b4_crop <- crop(stack_b4[i],AOI_Repro)

  b8_crop <- crop(stack_b8[i],AOI_Repro)

  ndvi <- (b4_crop[i]-b8_crop[i])/(b4_crop[i]-b8_crop[i])

  writeRaster(ndvi,output_name[i])

} 

And this the message I get from it:

Error in (function (classes, fdef, mtable) : unable to find an inherited method for function ‘crop’ for signature ‘"matrix"’

Answer 1

You are calling with single brackets, use double. Raster stack is a kind of list:

Reproducible example (different sensor, same procedure)

library(raster)
library(RStoolbox)

data(lsat)

output_name <- paste0('raster',1:10,'.tif')

stack_b4 <-list()
stack_b8 <-list()

for(i in 1:10){
  stack_b4[[i]] <- lsat[[3]]
  stack_b8[[i]] <- lsat[[4]]
}

stack_b4 <- stack(stack_b4)
stack_b8 <- stack(stack_b8)

AOI_Repro <- as(extent(lsat,50,150,50,150),'SpatialPolygons')

For your procedure:

for (i in 1:length(stack_b4)) {

  b4_crop <- crop(stack_b4[[i]],AOI_Repro)

  b8_crop <- crop(stack_b8[[i]],AOI_Repro)

  ndvi <- (b4_crop-b8_crop)/(b4_crop-b8_crop)

  writeRaster(ndvi,output_name[i])

} 

A better approach (IMO):

for (i in 1:length(stack_b4)) {

  b4_crop <- crop(stack_b4[[i]],AOI_Repro)

  b8_crop <- crop(stack_b8[[i]],AOI_Repro)

  ndvi <- overlay(x = b4_crop, y = b8_crop, fun = function(x,y) (y-x)/(y+x))

  writeRaster(ndvi,output_name[i])

} 

A faster approach:

NDVI <- function(x,y){(y-x)/(y+x)}

beginCluster(12) # your availables cores

for (i in 1:length(stack_b4)) {

  b4_crop <- crop(stack_b4[[i]],AOI_Repro)

  b8_crop <- crop(stack_b8[[i]],AOI_Repro)

  ndvi <- clusterR(stack(b4_crop,b8_crop), calc, args = list(fun = NDVI))

  writeRaster(ndvi,output_name[i])

} 

endCluster()

Less lines:

stack_b4c <- crop(stack_b4, AOI_Repro)
stack_b8c <- crop(stack_b8, AOI_Repro)
mapply(FUN = function(x,filename) writeRaster(x, filename),x = mapply(function(x,y) (y-x)/(y+x), as.list(stack_b4c),  as.list(stack_b8c)), filename=as.list(output_name))

Answer 2

First, what does this for( i in 1) loop do? i in 1 does not make any sense. Something like for(i in 1:length(stack_b4)) would be more appropriate here.

Second, you need double brackets and since you are creating a new object of index i, you do not need to call the iterator again. You could, however, crop your stacks first and not do this within the for loop.

for(i in for(i in 1:length(stack_b4))){
  b4_crop <- crop(stack_b4[[i]],AOI_Repro)
  b8_crop <- crop(stack_b8[[i]],AOI_Repro)
  ndvi <- (b4_crop-b8_crop)/(b4_crop-b8_crop)
  writeRaster(ndvi, output_name[i])
}