Rescaling vegetation index in R

I am having some trouble integrating a rescaling method within a function that is calculating a vegetation index for a raster. I tried using a formula from this solution. The code would run, but I would get two warning messages and my image would be blank. I checked the rasters min and max values and they read "-Inf" "Inf" respectively. I also tried another way using the RPMG library from this post, but ran into another error. This time after running the VARI variable. I want to keep the rescaling method as "bland" as possible so I can integrate it into other indices such as the Triangular Greeness Index (TGI). Any suggestions?

Method 1:

# Visable Atmospherically Resistant Index
VARI.Overlay <- function(b1, b2, b3){
  VARI.Calc <- (b1 - b3) / (b1 + b3 -b2)
  VARI.Scale <- ((VARI.Calc - min(VARI.Calc)) / (max(VARI.Calc) - min(VARI.Calc)) - 0.5 ) * 2
  return(VARI.Scale)
}
VARI <- overlay(img[[1]], img[[2]], img[[3]], fun = VARI.Overlay)

image(VARI, main = 'VARI')

Method 1 Error:

1: In min(x) : no non-missing arguments to min; returning Inf
2: In max(x) : no non-missing arguments to max; returning -Inf

Method 2:

# Visable Atmospherically Resistant Index
VARI.Overlay <- function(b1, b2, b3){
  VARI.Calc <- (b1 - b3) / (b1 + b3 -b2)
  VARI.min <- min(VARI.Calc)
  VARI.max <- max(VARI.Calc)
  VARI.Scale <- RESCALE(VARI.Calc, -1, 1, VARI.min, VARI.max)
  return(VARI.Scale)
}
VARI <- overlay(img[[1]], img[[2]], img[[3]], fun = VARI.Overlay)

Method 2 Error:

Error in (function (x, fun, filename = "", recycle = TRUE, forcefun = FALSE,  : 
  cannot use this formula, probably because it is not vectorized

Solution

You can do something like this

Example data

library(raster)
d <- brick(system.file("external/rlogo.grd", package="raster"))

Function. Note the !is.finite. This is to catch the cases where (b1 + b3 - b2) == 0. When that happens, the max value becomes Inf and the results are no good.

varifun <- function(b1, b2, b3){
  x <- (b1 - b3) / (b1 + b3 -b2)
  x[!is.finite(x)] <- NA
  x
}

v <- overlay(d, fun=varifun)

Now compute the min and max values. You must not do that in the function above, because that will operate on chunks of the dataset, and therefore each chunk may get different min and max values. Presumably you want these for the entire dataset.

vmn <- cellStats(v, "min", na.rm=TRUE)
vmx <- cellStats(v, "max", na.rm=TRUE)

And now combine

vari <- 2 * ((v - vmn) / (vmx - vmn) - 0.5)
vari
#class      : RasterLayer 
#dimensions : 77, 101, 7777  (nrow, ncol, ncell)
#resolution : 1, 1  (x, y)
#extent     : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
#crs        : +proj=merc +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
#source     : memory
#names      : layer 
#values     : -1, 1  (min, max)

By the way, the messages you get from Method 1 are not errors. These messages stem from computing min or max values for NA only.