3

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?

5

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]]))

plot

  • 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 [[1]] and [[27]] 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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.