Function for hyperspectral rasters with 3 variables unexpected error

Question

I am writing a script to extract vegetation indices from hyperspectral datacubes. One of the indices is the Triangular Vegetation Index for which I wrote the following function:

    TVI2 <- function(h,m,l) 
{
  (1.5*(2.5*(h-m)-1.3*(h-l))/sqrt((2*h+1)^2-(6*h-5*sqrt(m)-0.5)))
}

The variables stand for different wavelengths in the forumula. I have a rasterbrick object containing all 8 rasters and 101 bands per raster. To call a specific raster and wavelength I use:

raster(r[[1]], layer=1)

And this works just fine for calculating an index with just 2 variables, like the NDVI function:

ndvi1 <- NDVI(raster(r[[1]], layer=41),(raster(r[[1]], layer=71)))

However if I want to use this same set up for the TVI2 function:

tvi1 <- TVI2(raster(r[[1]], layer=71),(raster(r[[1]], layer=45),(raster(r[[1]], layer=21))))

I keep getting the following error:

Error: unexpected ',' in "tvi1 <- TVI2(raster(r[[1]], layer=71),(raster(r[[1]], layer=45),"

I cannot seem to figure out why. If I don' t use the raster call as the last variable, but just a number, the function runs fine.

Any ideas on how to solve this?

Answer 1

You need to pass your function to "overlay" along with the appropriate spectral bands in your raster brick. This example would pull the 71st, 45th and 21st bands of a brick object.

tvi1 <- overlay(r[[71]], r[[45]], r[[21]], fun = TVI2)

You may need to re-write your function so error conditions are accounted for.

I am not clear on your data structure, it seems odd and you do not provide enough information to understand how you formatted your data. Since syntax is often at the heart of R issues this is a problem. You should not need to coerce each band using raster. The bands in a raster stack/brick are already raster objects. If "r" is just a list of raster names then you should create an actual raster stack or brick. Are you storing a series of raster bricks in a list? If so, why?

Answer 2

So you have a list of 8 RasterBrick objects

r <- lapply(1:8, function(x) brick(system.file("external/rlogo.grd", package="raster")) )

And a function

TVI2 <- function(h,m,l) {
  (1.5*(2.5*(h-m)-1.3*(h-l))/sqrt((2*h+1)^2-(6*h-5*sqrt(m)-0.5)))
}

Now you can do

x <- TVI2(raster(r[[1]], layer=1), raster(r[[1]], layer=2), raster(r[[1]], layer=3))

Which is equivalent to

x <- TVI2(raster(r[[1]][[1]]), raster(r[[1]][[2]]), raster(r[[1]][[3]]))

This did not work for you because your parenthesis did not match. Or, following Jeffrey Evans:

y <- overlay(r[[1]][[1:3]], fun=TVI2)

with the indices you use (this won't work with the example data)

overlay(r[[1]][[c(71, 45, 21)]], fun = TVI2)

calling TVI2 with overlay should be more efficient than calling it directly with RasterLayer arguments