Rescaling vegetation index in R

Question

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

score 2 · Accepted Answer · answered May 27 '20 at 00:41

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.

Thank you for catching the division by zero. That makes sense. I am little confused still about why you cannot have the `vmn` and `vmx` within the function. I understand what you said, but why would it operate in chunks? — Binx, May 27 '20 at 02:05
The raster methods operate in chunks to allow processing of large datasets that cannot be loaded into memory. In this example, there is probably only 1 chunk, so it would not mater. — Robert Hijmans, May 27 '20 at 04:34

Rescaling vegetation index in R

1 Answers1