其他
R语言NDVI变异系数计算与结果分析
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")
那么，为什么前面两种方法求出来的CV那么大呢？
罪魁祸首就是这个cv函数，这个cv函数是raster包中的一个函数。
我再仔细看看帮助
原因就是，这个CV是百分比！
难怪我算的CV比文献里大那么多，原来是乘了100的
写代码还是要仔细看帮助啊，完全理解输出结果。
参考文献
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. 秦格霞, 芦倩, 孟治元, 等. 1982-2015年中国北方草地NDVI时空动态及其对气候变化的响应[J]. 水土保持研究, 2021, 28(01): 101-108+117. https://github.com/rspatial/terra https://rspatial.org/terra/rs/index.html https://www.rdocumentation.org/packages/raster/versions/3.4-5