怎么查看表用了那个序列_GIMMS NDVI数据R语言时间序列处理分析

GIMMS NDVI数据R语言时间序列处理分析

拿到的GIMMS NDVI数据是ENVI格式(不是NASA的原始数据)，要求做时间序列分析，查看变化情况，那么该怎么做呢？

R语言读取ENVI数据

R语言的raster包可以读取ENVI数据的dat文件，两种常用的函数，一个是brick，一个是stack

#在这里说明一下代码中用到的R包
library(raster)
library(tidyverse)
library(tibble)
library(stringr)
#读取数据，两种方法
NDVI_br "GIMMS_NDVI_CMR_1981_2015.dat")
NDVI_stack "GIMMS_NDVI_CMR_1981_2015.dat")

读取后查看数据，发现stack是分层的，每一个层里面都有独立的名称等信息，而brick就像它的名字，好像一个砖块，更像一个整体。

stack和brick对比

仔细查看一下stack，每一层里面有单独的数据，数据中有单独的名称。

stack详情

查看层数据文件名

R语言时间序列计算

计算思路

我采用下面文献5中提到的方法进行计算：

对15d数据采用平均值方法，计算得到月值数据，月值数据采用最大值合成(Maximum Value Composition,MVC)计算得到1982-2015年年值数据。

具体数据处理流程如下：

1. 逐stack layer计算NDVI均值
2. 计算相同月的NDVI均值，得月值数据
3. 对月值数据求最大值，计算年值数据

几个要点：

• NDVI计算结果转时间序列
• 月度均值计算
• 年度求最大值

逐图层计算均值

逐stack layer计算NDVI均值，代码如下：

avgNDVI #逐图层计算均值

计算结果为一个向量(vector)，带有图层的名称信息，也就是时间信息，但是我们需要进行一些格式转换，让R能够从名称中提取出时间。

avgNDVI

时间序列转换

要想把时间信息提出来，我们得弄一个数据框，让数据框中NDVI值和时间对起来。

avgNDVIdf #将计算结果转数据框
avgNDVIdf "Date")   #将行名转为一列
avgNDVIdf$Date $Date, "X")   #去除X
avgNDVIdf$Date2 $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 $Date2)  #提取月信息
#函数强制转换月份汉语为数字
numMonthfunction(x) 
  c(一月="01",二月="02",三月="03",四月="04",五月="05",六月="06",七月="07",八月="08",九月="09",十月=10,十一月=11,十二月=12)[tolower(x)]
avgNDVIdf$Month $Month)
avgNDVIdf$Year $Date2, format="%Y")  #提取年份信息
Mon_Mean_NDVI #计算月度均值
Mon_Mean_NDVI$Date $Year, Mon_Mean_NDVI$Month,"15", sep = "-"), format = "%Y-%m-%d")  #月度15号
Year_Max_NDVI Year_Max_NDVI$Date $Year, "07", "01", sep = "-"), format = "%Y-%m-%d")   #每年7月1日

数据计算结果1

修改时间序列代码后，添加Date列

R语言时间序列图

推荐使用ggplot2进行NDVI时间变化曲线的绘制，代码示例如下：

p2   labs(x="Year", y="NDVI")+
  theme_bw()+
  geom_smooth()+
  geom_point()
p2

ggplot2绘制时间序列变化

R语言栅格动图

R语言可以使用raster包中的anmimate函数生成动图，但是我不太清楚如何输出，我是用Screen to GIF软件截得屏

animate(NDVI_br)

NDVI时间序列图

参考文献

1. https://www.neonscience.org/dc-raster-time-series-r
2. https://www.neonscience.org/julian-day-conversion-r
3. https://www.r-graph-gallery.com/279-plotting-time-series-with-ggplot2.html
4. https://www.neonscience.org/dc-ndvi-calc-raster-time-series
5. 王涛,赵元真,王慧,曹亚楠,彭静,曹亚楠.基于GIMMS NDVI的青藏高原植被指数时空变化及其气温降水响应[J].冰川冻土,2020,42(02):641-652.
6. https://www.cnblogs.com/ljhdo/p/4804113.html
7. https://www.neonscience.org/dc-convert-date-time-POSIX-r
8. https://stackoverflow.com/questions/16652199/compute-monthly-averages-from-daily-data