R语言NDVI变异系数计算与结果分析

变异系数是衡量观测序列值变异程度的一个统计量，可以很好地反映空间数据在时间序列上变化的差异程度，评价数据时间序列的稳定性。计算公式为：

式中： 为标准差， 为算术平均值。CV越大，数据波动越大。

一般情况下，将区域NDVI CV划分五个等级，来表征区域植被覆盖稳定情况：

CV	稳定性
<0.05	低波动
0.05≤CV<0.10	相对低波动
0.10≤CV< 0.15	中等波动
0.15  ≤CV<0.20	相对高波动
≥0.20	高波动

接下来介绍几种不同的代码，还有我踩的坑。

R语言变异系数计算

栅格转数据框逐像元计算

栅格的本质就是矩阵，在R语言里面可以把栅格转为数据框或者tibble，然后进行数值计算，再将计算结果转回栅格。下面是代码：

#合成栅格砖
flnames <- list.files(path = './monmean/', pattern = '.tif$')
fl <- paste0("./monmean/", flnames)
Raster_MonMean <- brick(lapply(fl, raster))
writeRaster(Raster_MonMean, "./月均值多波段分析/MonMeanstack.tif", format="GTiff")

#月度均值CV
raster('./月均值多波段分析/MonMeanstack.tif', band = 2) -> test

test[] %>% 
  as_tibble()

#每个波段读取为数据框
lapply(1:408, function(x) {
  raster('./月均值多波段分析/MonMeanstack.tif', 
         band = x) %>% 
    `[` %>% 
    as_tibble() %>% 
    mutate(id = as.numeric(row.names(.)))
}) %>% 
  bind_rows() -> df

df
df %>% 
  arrange(id) %>% 
  nest(-id) -> nested_df

nested_df$data[1][[1]] %>% 
  pull(value) %>% 
  cv()  #变异系数计算

nested_df %>% 
  mutate(cv = map_dbl(data, function(x) cv(x$value))) %>% 
  dplyr::select(-data) %>% 
  pull(cv) -> cv
cv %>% 
  raster(vals = ., nrow = 295, ncol = 1366,
         crs = "+proj=longlat +datum=WGS84 +no_defs",
         ext = extent(test)) -> cv_raster

plot(cv_raster)

# 输出
cv_raster %>% 
  writeRaster("./月均值多波段分析/MonMeanCV.tif", format="GTiff", overwrite = TRUE)

使用Terra包计算

最近发现一个新的R包，Terra,支持并行运算，比Raster包更加强大，也可以用来计算变异系数。

library(terra)
library(raster)
GIMMS <- rast("./Yearmaxdata/YearMaxstack.tif")  #读取GIMMS数据
sr <- app(GIMMS, cv)
plot(sr)
writeRaster(sr, filename = "./Yearmaxdata/YearMaxCV_terra.tif", filetype="GTiff")

计算结果分析

上面的代码计算完变异系数后，CV值域从2.2到583

变异系数计算结果

再看看文献中写的，0.2以上就是高波动变化了，而我计算的是2.2起步！这是什么情况呢？数据有问题？还是计算方法的问题呢？

植被稳定性空间格局

CV就是标准差除均值，计算原理很简单，为了检查我到底是数据问题还是程序算法问题，我使用Terra包分步计算了CV。

Terra包分步计算CV

#前面数据和包均已导入，继续上面的代码的计算
stdrast <- app(GIMMS, sd, cores=4)  #计算标准差
plot(stdrast)
meanrast <- app(GIMMS, mean)  #计算均值
plot(meanrast)
cvrast <- stdrast/meanrast  #CV计算，简单粗暴，标准差均值直接相除
writeRaster(cvrast, filename = "./Yearmaxdata/YearMaxCV_terra2.tif")

标准差和均值都[0,1]之间

那么，为什么前面两种方法求出来的CV那么大呢？

罪魁祸首就是这个cv函数，这个cv函数是raster包中的一个函数。

罪魁祸首就是这个cv函数

我再仔细看看帮助

原因就是，这个CV是百分比！

难怪我算的CV比文献里大那么多，原来是乘了100的

写代码还是要仔细看帮助啊，完全理解输出结果。

参考文献

1. Jiang W, Yuan L, Wang W, etal. Spatio-temporal analysis of vegetation variation in the Yellow River Basin[J]. Ecological Indicators, 2015, 51: 117–126.
2. 秦格霞, 芦倩, 孟治元, 等. 1982-2015年中国北方草地NDVI时空动态及其对气候变化的响应[J]. 水土保持研究, 2021, 28(01): 101-108+117.
3. https://github.com/rspatial/terra
4. https://rspatial.org/terra/rs/index.html
5. https://www.rdocumentation.org/packages/raster/versions/3.4-5