1

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?

2 Answers 2

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?

2
  • Thankyou for your reply, however the code does not pull the bands of the brick object, but rather the 71th, 45th and 21th raster of the brick. The data structure of r is a rasterbrick containing 8 rasters, each with 101 layers: files <- list.files(path=".", pattern = ".bsq$") r <- lapply(files, brick). To call a certain raster and layer I use: r1l41 <- raster(r[[1]], layer=41) which then contains the first raster of the brick and the 41th band. Is there another, better way, to call these bands and use them in a function like the TVI example? Many thanks!
    – cf2
    Nov 3, 2015 at 19:00
  • You do not need to use an lapply to create the brick. What you are ending up with is a list object containing raster bricks. If you would like to keep that data structure there is no problem but you need to understand the object structure that you are working with. In my example just use an index for the raster brick, contained in the list object, and an index for the band ie., to access the first brick object and the specific bands use: overlay(r[[1]][[71]], r[[1]][[45]], r[[1]][[21]], fun = TVI2) Nov 3, 2015 at 20:44
0

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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