利用raster包进行栅格数据处理（3）

      ,欢迎也在微信上查看相关内容。

有关栅数据处理的更多信息，还可查看第1部分、第2部分。本次主要利用Peneda-Geres国家公园(PGNP，位于葡萄牙西北部)的5年2012年至2016年)增强植被指数(EVI)来进行一些分析。该数据自MODIS的MOD13Q1数据产品（版本-006），其空间分辨率为250m，时间分辨率为16天(更准确地说，该产品是由每个16天周期的最大日值复合生成的)。这意味着每年总共有23次观测。该数据是从EarthData平台下载的，随后使用MODIS重投影工具重新投影到WGS 1984 - UTM 29N坐标参考系统(CRS)。每年有23张影像，总共有115个图层。MODIS数据产品，如MOD13Q1，使用由年份(4位)和儒略日(3位，即当天位于这一年的第几天)组成的7位长日期来确定图像合成的参考日期。例如，日期码2012001对应的日期为2012-01-01(YYYY-mm-dd)和2012161 为 2012-06-09。通常，这些日期被嵌入到图像文件或元数据中，首先生成它们，然后对它们进行处理，这样我们就可以为每一层获得一个日期对象。然后我们可以在rts函数中使用这些正确格式化的日期来创建一个raster时间序列。

  RasterBrick是一个典型的由多层(或多波段)文件创建的多层raster对象，与RasterStacks类似，但使用RasterBrick时处理时间更短。然而，它只能指向单个文件，而RasterStacks可以指向多个文件。

    除了raster包，我们还将使用rts，它提供了操作和处理raster时间序列数据的类和方法(例如，卫星图像的时间序列)。raster时间序列对象是通过组合raster包中的RasterStack或RasterBrick对象和POSIXct、POSIXt、Date、timeDate类的一组日期来创建的。然后，rts中的时间信息由xts对象处理。

rts函数用于构建raster时间序列(RasterBrickTS或RasterStackTS)，即包含两个槽的S4对象:

    Slot [raster]: RasterStack或RasterBrick对象。

   Slot[time]:一个xts对象，为栅格对象中的每一层设置日期。

下载地址：

https://raw.githubusercontent.com/joaofgoncalves/R_exercises_raster_tutorial/master/data/MODIS_EVI_TS_PGNP_MultiBand.zip

#---时间序列------

library(raster)

setwd("F:/test")

rst <- brick("data-raw/MOD13Q1.2012_2016.PGNP_250m_EVI_16days.tif")

names(rst) <- paste("EVI",1:nlayers(rst),sep="_")

padIt <- function(x)

  {if(x<10)

    paste("00",x,sep="")

  else if(x<100 & x>=10)

    paste("0",x, sep="")

else

  return(as.character(x))}

padWithZeros <- function(x) sapply(x, FUN = padIt)

# Generate the MODIS-like dates for each layer

julDay <- padWithZeros(rep(seq(from = 1,to = 365,by = 16),5))

yrs <- as.character(rep(2012:2016, each = 23))

MODISYrJday <- paste(yrs, julDay, sep="")

# Extract the year

MOD.getYear<-function(x){

  as.integer(sapply(x,FUN=function(x) substr(x,1,4)))}

# Extract de julian day

MOD.getDOY<-function(x) {

  as.integer(sapply(x,FUN=function(x) substr(x,5,7)))}

# Process the MODIS-like date into YYYY-mm-dd format as Date object

MOD.getDate<-function(x){

  as.Date(sapply(x,FUN=function(x) as.character(as.Date(MOD.getDOY(x)-1,origin=paste(MOD.getYear(x),"01-01",sep="-")))))}

MODdates <- MOD.getDate(MODISYrJday)

#install.packages("rts")

library(rts)

rstTS <- rts(rst, MODdates)

#Sub-setting Raster Time-series

# Subset a specific period

rstTSsubset1 <- subset(rstTS,"2013-05-15/2014-08-25")

# Subset the whole year of 2012

rstTSsubset2 <- subset(rstTS,"2012")

# Subset years from 2013 to 2014

rstTSsubset3 <- subset(rstTS,"2013/2014")

# Subset all years from (and including) 2014 to the series end

rstTSsubset4 <- subset(rstTS,"2014/")

# Subset all to the end of 2014

rstTSsubset5 <- subset(rstTS,"/2014")

# Subset all to the end of July 2014

rstTSsubset6 <- subset(rstTS,"/2014-06")

# Subset a specific month

rstTSsubset7 <- subset(rstTS,"2016-05")

# Plot the May 2016 data

plot(rstTSsubset7)

# 计算每个季度的平均值

# Mean

rstTS_quarterlyMN <- apply.quarterly(rstTS, FUN = mean, na.rm=TRUE)

# Standard-deviation

rstTS_quarterlySD <- apply.quarterly(rstTS, FUN = sd, na.rm=TRUE)

# Apply function to each year

# Mean

rstTS_yearlyMN <- apply.yearly(rstTS, FUN = mean, na.rm=TRUE)

# Standard-deviation

rstTS_yearlySD <- apply.yearly(rstTS, FUN = sd, na.rm=TRUE)

 

# Plot the time-series for annual EVI

plot(rstTS_yearlyMN)

如果希望raster时间序列对象使用raster包提供的“更简单”/更通用的函数，可使用calc和stackApply函数。这些函数之间的区别是，calc逐像素应用定义的函数，而stack对RasterStack或RasterBrick的子集应用函数。对于这个函数，输出结果图层名以index_加序号。所使用的函数应该返回单个值，并且输出rasterc中的层数等于索引中的惟一值的数量。另一方面，calc允许定义多个计算类型如mean和sum，而stackApply只能指定一个计算方式。

# Apply function to each year

# Mean

rstTS_yearlyMN <- apply.yearly(rstTS, FUN = mean, na.rm=TRUE)

# Standard-deviation

rstTS_yearlySD <- apply.yearly(rstTS, FUN = sd, na.rm=TRUE)

 

# Plot the time-series for annual EVI

plot(rstTS_yearlyMN)

rstMean <- calc(rst, fun = mean)

rstQuantiles <- calc(rst, fun = function(x,...) {

  as.numeric(quantile(x,probs=c(0.05, 0.5, 0.95),...))},

  na.rm=TRUE)

#每年的均值，计算速度比apply.yearly函数快

rstYrMean <- stackApply(rst, fun=mean,

                        indices = rep(1:5,each=23))

# stackApply with RasterBrick的时间

system.time({rstYrMean <- stackApply(rst, fun=mean, indices = rep(1:5,each=23))})

# apply.yearly with RasterBrickTS的时间

system.time({rstTS_yearlyMN <- apply.yearly(rstTS, FUN = mean, na.rm=TRUE)})

plot(rstYrMean)