重新缩放 R 中的植被指数

Question

我在将重新缩放方法集成到计算栅格植被指数的函数中时遇到了一些问题。我尝试使用 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 中的公式，但运行变成了另一个错误。这次在运行之后 VARI 变量。我希望尽可能将重新缩放方法保持为 "bland"，以便我可以将其集成到其他指数中，例如三角绿度指数 (TGI)。有什么建议么？

方法一：

# 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')

方法一错误：

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

方法二：

# 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)

方法二错误：

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

Answer 1

你可以这样做

示例数据

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

函数。注意 !is.finite。这是为了捕获 (b1 + b3 - b2) == 0 的情况。发生这种情况时，最大值变为 Inf，结果不好。

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

v <- overlay(d, fun=varifun)

现在计算最小值和最大值。您必须 而不是 在上面的函数中这样做，因为这将对数据集的块进行操作，因此每个块可能会获得不同的最小值和最大值。大概您想要整个数据集的这些。

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

现在合并

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)

顺便说一下，您从方法 1 获得的消息不是错误。这些消息源于仅计算 NA 的最小值或最大值。

重新缩放 R 中的植被指数

Rescaling vegetation index in R

r

raster

geospatial

rescale