GIMMS和MODIS NDVI时间序列对比分析

前面我们进行了GIMMS和MODIS NDVI的月度数据合成，接下来进行一下时间序列的对比分析。

数据准备

• 栅格数据裁剪，使用terra包中trim和mask函数
• MODIS NDVI数据需要除以10000，将值域转换为[-1,1]
• seq和as.Date生成时间序列，作为栅格波段名称
• writeRaster输出栅格
• animate可以输出动画，使用Screen to GIF录制GIF图，n设置循环次数

library(terra)
cngimmslst = list.files(path = "D:/R/CMR/ChinaMonthMax", pattern = ".tif$")
cngimms = paste0("D:/R/CMR/ChinaMonthMax/", cngimmslst)
cngimms_rast = rast(cngimms)
#获取京津冀范围
JJJ = vect("./SHP/JJJ.shp")
gimms_mon_msk = trim(mask(cngimms_rast, JJJ))
#时间
gimms_temp = seq(as.Date("1982-01-01"), as.Date("2015-12-01"), "month")
names(gimms_mon_msk) = gimms_temp
writeRaster(gimms_mon_msk, filename = "./MVC/gimms_mvc.tif", names=gimms_mon_msk@ptr[["names"]])

#MODIS数据读取
modlst = list.files(path = "./GEEMVC1k", pattern = ".tif$")
mod = paste0("./GEEMVC1k/", modlst)
mod_rast= rast(mod)
#时间
mod_temp = seq(as.Date("2001-01-01"), as.Date("2020-12-01"), "month")
names(mod_rast) = mod_temp
#京津冀范围裁剪
mod_JJJ = trim(mask(mod_rast, JJJ))
#MODIS NDVI有系数
mod_JJJ = mod_JJJ/10000
writeRaster(mod_JJJ, filename = "./MVC/mod_mvc.tif", names=mod_JJJ@ptr[["names"]])
#NDVI数据查看
plot(mod_JJJ[[1:12]])
animate(mod_JJJ[[1:12]], n=99)

GIMMS NDVI 1982年时间序列

MODIS NDVI 2001年时间序列

MODIS NDVI2001年变化

栅格统计生成时间序列并制图

栅格统计

• global进行像元统计，在这里计算了整个像元的均值

#获取重合时间2001-2015的MODIS和GIMMS数据
t1 = seq(as.Date("1982-01-01"), as.Date("2001-01-01"), "month")
length(t1)   #获取2001年1月的波段号
print(length(t1)+12*15)   #获取2015年12月波段号
#获取2001-2015的GIMMS
gimms_rs = gimms_mon_msk[[229:409]]
#获取2001-2015年MODIS
mod_rs = mod_JJJ[[1:180]]

library(tidyverse)
#MODIS最大值时间序列
MODIS_Mon_Series = global(mod_rs, "mean", cores=4, na.rm=TRUE)
MODIS_Mon_Series = as.data.frame(MODIS_Mon_Series, make.names=TRUE)
colnames(MODIS_Mon_Series)= "MODIS"  #修改列名
MODIS_Mon_Series = rownames_to_column(MODIS_Mon_Series, "temp")

# 月度NDVI时间序列
#月最大值的时间序列
GIMMS_Mon_Series <- global(gimms_rs, "mean", cores=4, na.rm=TRUE)
colnames(GIMMS_Mon_Series) = "GIMMS"  #修改列名
GIMMS_Mon_Series = rownames_to_column(GIMMS_Mon_Series, "temp")

时间序列与制图

• 由于用到了两种NDVI数据，为了对比，需要把俩NDVI列合起来
• full_join将两组NDVI合并
• 关于R的几种数据链接形式可以参考这个文章：https://www.jianshu.com/p/1f4c7bfed3d4
• 里面的一些动图出自于https://github.com/gadenbuie/tidyexplain#relational-data
• ggplot2制图

#时间序列合并
Mon_series = full_join(GIMMS_Mon_Series, MODIS_Mon_Series, by="temp")

#数据转换
MonthlyNDVIdf = gather(Mon_series, key = "var", value = "NDVI", - temp)
MonthlyNDVIdf$temp = as.Date(MonthlyNDVIdf$temp)

#绘制月度变化图
NDVI.formula <- y~x
p1 <- ggplot(data = MonthlyNDVIdf, aes(temp, NDVI, color=var))+
  labs(x="Month", y="NDVI")+
  theme_bw()+
  theme(title = element_text(size = 22),
        axis.title = element_text(size = 20),    #调整标题大小
        axis.text.x = element_text(size = 18),    #x轴标签大小
        axis.text.y = element_text(size = 18), 
        legend.position = "bottom",
        legend.text = element_text(size = 18),
        legend.title = element_blank())+     #y轴标签大小
  geom_point()+
  geom_line()
p1

GIMMS NDVI和MODIS NDVI时间序列对比

参考文献

点击阅读原文查看视频课程学习