拿到的GIMMS NDVI數據是ENVI格式（不是NASA的原始數據），要求做時間序列分析，查看變化情況，那麼該怎麼做呢？
ENVI裁剪數據R語言裁剪數據效率不高，推薦使用ENVI進行數據的預處理，先把數據按照需要進行一些裁剪，方便後續分析。
在ENVI中，裁剪數據之前需要先劃定感興趣區（Region of Interest, ROI），在這我們先使用ROI工具對研究區範圍進行劃定。ROI可以手動劃定，也可以在ROI工具中導入SHP文件。
新建ROI在這裡我直接畫了一個矩形，需要注意的是，畫完ROI後，需要接受一下，讓區域顏色變成不透明的，這樣就生成了一個ROI
接受ROI如下所示，ROI區域變成不透明的顏色後，這個ROI就繪製完成了，這個時候就可以使用裁剪工具（Subset Data from ROIs）進行裁減了。
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和brick對比仔細查看一下stack，每一層裡面有單獨的數據，數據中有單獨的名稱。
stack詳情查看層數據文件名R語言時間序列計算計算思路我採用下面文獻5中提到的方法進行計算：
對15d數據採用平均值方法，計算得到月值數據，月值數據採用最大值合成（Maximum Value Composition,MVC）計算得到1982-2015年年值數據。
具體數據處理流程如下：
幾個要點：
逐圖層計算均值逐stack layer計算NDVI均值，代碼如下：
avgNDVI <- cellStats(NDVI_stack, mean) #逐圖層計算均值計算結果為一個向量（vector），帶有圖層的名稱信息，也就是時間信息，但是我們需要進行一些格式轉換，讓R能夠從名稱中提取出時間。
avgNDVI時間序列轉換要想把時間信息提出來，我們得弄一個數據框，讓數據框中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年%m：兩位數字的月份，取值範圍是01-12，或1-12%B：英語月份全名（January、February 、March等）%a：縮寫的星期名（Mon、Tue、Wed、Thur、Fri、Sat、Sun）運行上述代碼後，我們實現了NDVI和時間的對應，下圖的數據框中：
帶時間序列的數據框月度年度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日
數據計算結果1修改時間序列代碼後，添加Date列R語言時間序列圖推薦使用ggplot2進行NDVI時間變化曲線的繪製，代碼示例如下：
p2 <- ggplot(data = Year_Max_NDVI, aes(Date, avgNDVI))+
labs(x="Year", y="NDVI")+
theme_bw()+
geom_smooth()+
geom_point()
p2
ggplot2繪製時間序列變化R語言柵格動圖R語言可以使用raster包中的anmimate函數生成動圖，但是我不太清楚如何輸出，我是用Screen to GIF軟體截得屏
animate(NDVI_br)
NDVI時間序列圖參考文獻https://www.neonscience.org/dc-raster-time-series-rhttps://www.neonscience.org/julian-day-conversion-rhttps://www.r-graph-gallery.com/279-plotting-time-series-with-ggplot2.htmlhttps://www.neonscience.org/dc-ndvi-calc-raster-time-series王濤,趙元真,王慧,曹亞楠,彭靜,曹亞楠.基於GIMMS NDVI的青藏高原植被指數時空變化及其氣溫降水響應[J].冰川凍土,2020,42(02):641-652.https://www.cnblogs.com/ljhdo/p/4804113.htmlhttps://www.neonscience.org/dc-convert-date-time-POSIX-rhttps://stackoverflow.com/questions/16652199/compute-monthly-averages-from-daily-data