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?

Improve this question

1 Answer 1

Reset to default
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

Improve this answer
2
  • 1
    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, 2017 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, 2017 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.