GIMMS NDVI数据ENVI裁剪和R语言时间序列处理分析
GIMMS NDVI数据R语言时间序列处理分析
拿到的GIMMS NDVI数据是ENVI格式（不是NASA的原始数据），要求做时间序列分析，查看变化情况，那么该怎么做呢？
ENVI裁剪数据
R语言裁剪数据效率不高，推荐使用ENVI进行数据的预处理，先把数据按照需要进行一些裁剪，方便后续分析。
在ENVI中，裁剪数据之前需要先划定感兴趣区（Region of Interest, ROI），在这我们先使用ROI工具对研究区范围进行划定。ROI可以手动划定，也可以在ROI工具中导入SHP文件。
在这里我直接画了一个矩形，需要注意的是，画完ROI后，需要接受一下，让区域颜色变成不透明的，这样就生成了一个ROI
如下所示，ROI区域变成不透明的颜色后，这个ROI就绘制完成了，这个时候就可以使用裁剪工具（Subset Data from ROIs）进行裁减了。
接下来选择裁剪范围，要注意的是，Mask pixels outside of ROI这里要选择Yes，否则的话裁剪范围外还会有像元。然后选择保存位置，即可完成裁剪。
R语言读取ENVI数据
R语言的raster包可以读取ENVI数据的dat文件，两种常用的函数，一个是brick，一个是stack
#在这里说明一下代码中用到的R包
library(raster)
library(tidyverse)
library(tibble)
library(stringr)
#读取数据，两种方法
NDVI_br <- brick("GIMMS_NDVI_CMR_1981_2015.dat")
NDVI_stack <- stack("GIMMS_NDVI_CMR_1981_2015.dat")
读取后查看数据，发现stack是分层的，每一个层里面都有独立的名称等信息，而brick就像它的名字，好像一个砖块，更像一个整体。
仔细查看一下stack，每一层里面有单独的数据，数据中有单独的名称。
R语言时间序列计算
计算思路
我采用下面文献5中提到的方法进行计算：
对15d数据采用平均值方法，计算得到月值数据，月值数据采用最大值合成（Maximum Value Composition,MVC）计算得到1982-2015年年值数据。
具体数据处理流程如下：
逐stack layer计算NDVI均值 计算相同月的NDVI均值，得月值数据 对月值数据求最大值，计算年值数据
几个要点：
NDVI计算结果转时间序列 月度均值计算 年度求最大值
逐图层计算均值
逐stack layer计算NDVI均值，代码如下：
avgNDVI <- cellStats(NDVI_stack, mean) #逐图层计算均值
计算结果为一个向量（vector），带有图层的名称信息，也就是时间信息，但是我们需要进行一些格式转换，让R能够从名称中提取出时间。
时间序列转换
要想把时间信息提出来，我们得弄一个数据框，让数据框中NDVI值和时间对起来。
avgNDVIdf <- as.data.frame(avgNDVI) #将计算结果转数据框
avgNDVIdf <- rownames_to_column(avgNDVIdf, "Date") #将行名转为一列
avgNDVIdf$Date <- str_remove(avgNDVIdf$Date, "X") #去除X
avgNDVIdf$Date2 <- as.Date(avgNDVIdf$Date, format = '%Y.%m.%d')
format内部代码解释：
%y：两位数字表示的年份（00-99），不带世纪，例如，数值是18，格式%y，表示2018年 %Y：四位数字表示的年份（0000-9999） %m：两位数字的月份，取值范围是01-12，或1-12 %d：月份中的天，取值范围是01-31 %e：月份中的天，取值范围是1-31 %b：缩写的月份（Jan、Feb、Mar等） %B：英语月份全名（January、February 、March等） %a：缩写的星期名（Mon、Tue、Wed、Thur、Fri、Sat、Sun） %A：星期全名
运行上述代码后，我们实现了NDVI和时间的对应，下图的数据框中：
Date列为日期（字符串型），R语言不能识别 avgNDVI为计算的各个图层的NDVI均值 Date2为日期，R语言可以识别的形式
月度年度NDVI计算
时间序列生成后，我们要根据时间序列计算每月均值和每年的最大值了。
最开始我只是计算了月度均值和年度最大值，数据框中加上了年份和月份，但是由于最后制图会出问题，所以后来又调整代码，将月度和年度调整为了完整的月度和年度的时间序列。
以下是完整的NDVI年度和月度值计算和时间序列生成代码：
#月均NDVI和年度每月最大NDVI计算
avgNDVIdf$Month <- months(avgNDVIdf$Date2) #提取月信息
#函数强制转换月份汉语为数字
numMonth<-function(x)
c(一月="01",二月="02",三月="03",四月="04",五月="05",六月="06",七月="07",八月="08",九月="09",十月=10,十一月=11,十二月=12)[tolower(x)]
avgNDVIdf$Month <- numMonth(avgNDVIdf$Month)
avgNDVIdf$Year <- format(avgNDVIdf$Date2, format="%Y") #提取年份信息
Mon_Mean_NDVI <- aggregate(avgNDVI ~ Month+Year, avgNDVIdf, mean) #计算月度均值
Mon_Mean_NDVI$Date <- as.Date(paste(Mon_Mean_NDVI$Year, Mon_Mean_NDVI$Month,"15", sep = "-"), format = "%Y-%m-%d") #月度15号
Year_Max_NDVI <- aggregate(avgNDVI ~ Year, Mon_Mean_NDVI, max)
Year_Max_NDVI$Date <- as.Date(paste(Year_Max_NDVI$Year, "07", "01", sep = "-"), format = "%Y-%m-%d") #每年7月1日
R语言时间序列图
推荐使用ggplot2进行NDVI时间变化曲线的绘制，代码示例如下：
p2 <- ggplot(data = Year_Max_NDVI, aes(Date, avgNDVI))+
labs(x="Year", y="NDVI")+
theme_bw()+
geom_smooth()+
geom_point()
p2
R语言栅格动图
R语言可以使用raster包中的anmimate函数生成动图，但是我不太清楚如何输出，我是用Screen to GIF软件截得屏
animate(NDVI_br)
参考文献
https://www.neonscience.org/dc-raster-time-series-r https://www.neonscience.org/julian-day-conversion-r https://www.r-graph-gallery.com/279-plotting-time-series-with-ggplot2.html https://www.neonscience.org/dc-ndvi-calc-raster-time-series 王涛,赵元真,王慧,曹亚楠,彭静,曹亚楠.基于GIMMS NDVI的青藏高原植被指数时空变化及其气温降水响应[J].冰川冻土,2020,42(02):641-652. https://www.cnblogs.com/ljhdo/p/4804113.html https://www.neonscience.org/dc-convert-date-time-POSIX-r https://stackoverflow.com/questions/16652199/compute-monthly-averages-from-daily-data