0.1 R Markdown
1 load required libraries and global functions
2 load all original data
3 Transfer all data to data frame
4 Transfer monthly climate data to annual
5 Transfer monthly climate data to season (preseason, growseason,posseason)
6 Staistic all variable for each vegetation type and plot
- 6.1 sta mean value of Phe. for each grid
- 6.2 sta mean time series value of Phe. for each vegetation type
7 Grid cor analysis of climate and Phe.
- 7.1 Mean annual relationship
- 7.2 Seasonal relationship

0.1 R Markdown

This is Script for processing SGS of Zagunao watershed for three main vegetation types.

1 load required libraries and global functions

# check, load and install libraries
# load common functions
#source("http://146.118.107.77/Functions.R")
#source("/mnt/Funcs_AU_PRE.R")

if ("dWaSSI" %in% installed.packages()) {library(dWaSSI)} else{install_github("ln1267/dWaSSI") }

librs<-c("ggplot2","reshape2","pryr","caTools","raster","parallel","RColorBrewer","abind")

f_lib_check(librs)

2 load all original data

The spatial resolution of data is 250m, Phe. data temporal resolution is annual. While the temporal resolution of Climate data is monthly.

DEM<-read.ENVI("ZGNPhen/ZGNDEM")
VEG<-read.ENVI("ZGNPhen/ZGNforests")
#plot(raster(VEG))
SGS<-read.ENVI("ZGNPhen/Phenology/zgnSGS")
GSL<-read.ENVI("ZGNPhen/Phenology/zgnGSL")
END<-read.ENVI("ZGNPhen/Phenology/zgnEND")

nrows<-nrow(DEM)
ncols<-ncol(DEM)

TEM_mon<-array(NA,dim=c(nrows,ncols,192))
PRE_mon<-array(NA,dim=c(nrows,ncols,192))

TEM_mon[,,1:132]<-read.ENVI("ZGNPhen/TEM")
TEM_mon[1:293,1:255,133:192]<-read.ENVI("ZGNPhen/TEM_ZGN_11-15") 
PRE_mon[,,1:132]<-read.ENVI("ZGNPhen/PRE")
PRE_mon[1:293,1:255,133:192]<-read.ENVI("ZGNPhen/PRE_ZGN_11-15")

3 Transfer all data to data frame

Vegetation info : 1- Planted; 2 - Mixed; 3- Natural
Fill value : 0

Data_ZGN<-data.frame(ID=c(1:75264),Year=rep(c(2000:2014),each=75264),DEM=as.vector(DEM),VEG=as.vector(VEG),SGS=as.vector(SGS),GSL=as.vector(GSL),END=as.vector(END))

summary(Data_ZGN)

##        ID             Year           DEM            VEG        
##  Min.   :    1   Min.   :2000   Min.   :   0   Min.   :0.0000  
##  1st Qu.:18817   1st Qu.:2003   1st Qu.:3420   1st Qu.:0.0000  
##  Median :37632   Median :2007   Median :3981   Median :0.0000  
##  Mean   :37632   Mean   :2007   Mean   :3737   Mean   :0.5911  
##  3rd Qu.:56448   3rd Qu.:2011   3rd Qu.:4395   3rd Qu.:1.0000  
##  Max.   :75264   Max.   :2014   Max.   :5714   Max.   :3.0000  
##       SGS             GSL             END       
##  Min.   :  0.0   Min.   :  0.0   Min.   :  0.0  
##  1st Qu.:  0.0   1st Qu.:  0.0   1st Qu.:  0.0  
##  Median :146.6   Median :123.9   Median :285.9  
##  Mean   :107.6   Mean   :104.2   Mean   :206.3  
##  3rd Qu.:167.3   3rd Qu.:151.2   3rd Qu.:299.8  
##  Max.   :210.0   Max.   :300.0   Max.   :350.0

# set fill value to NA
for (i in c(3:7)){

  Data_ZGN[[i]][Data_ZGN[[i]]==0]<-NA
}

summary(Data_ZGN)

##        ID             Year           DEM             VEG        
##  Min.   :    1   Min.   :2000   Min.   :1700    Min.   :1.0     
##  1st Qu.:18817   1st Qu.:2003   1st Qu.:3518    1st Qu.:2.0     
##  Median :37632   Median :2007   Median :4017    Median :2.0     
##  Mean   :37632   Mean   :2007   Mean   :3925    Mean   :2.1     
##  3rd Qu.:56448   3rd Qu.:2011   3rd Qu.:4411    3rd Qu.:3.0     
##  Max.   :75264   Max.   :2014   Max.   :5714    Max.   :3.0     
##                                 NA's   :54015   NA's   :817830  
##       SGS              GSL              END        
##  Min.   : 90.0    Min.   : 50.0    Min.   :210.0   
##  1st Qu.:146.0    1st Qu.:118.4    1st Qu.:283.9   
##  Median :159.8    Median :135.6    Median :294.7   
##  Mean   :159.4    Mean   :140.5    Mean   :295.0   
##  3rd Qu.:173.2    3rd Qu.:162.0    3rd Qu.:305.7   
##  Max.   :210.0    Max.   :300.0    Max.   :350.0   
##  NA's   :367005   NA's   :291675   NA's   :339675

4 Transfer monthly climate data to annual

# Function for transfering monthly grid (array) to annual, by selecting "method"= "mean" or "sum"
f_mon2ann_grid<-function(da,method="total") {
  da_ann<-array(0,c(nrow(da),ncol(da),dim(da)[3]/12))
  a<-1
  b<-1
    for (y in 1:(dim(da)[3]/12)){
      linshi<-matrix(0,ncol=ncol(da),nrow=nrow(da))
      for (m in 1:12){
        
        linshi<-linshi+da[,,a]
        #print(range(linshi,na.rm=T))
        a<-a+1
      }
     
      da_ann[,,b]<-linshi
      
      b<-b+1
      print(y)
    }
    range(da_ann,na.rm=T)
    
  if (method=="Mean" |method=="MEAN" |method=="mean"){
    da_ann<-da_ann/12
  }
    da_ann
}

PRE_ann<-f_mon2ann_grid(PRE_mon)

## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12
## [1] 13
## [1] 14
## [1] 15
## [1] 16

summary(PRE_ann)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##     0.0   977.5  1071.0  1034.0  1158.0  1609.0    2745

TEM_ann<-f_mon2ann_grid(TEM_mon,method = "mean")

## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12
## [1] 13
## [1] 14
## [1] 15
## [1] 16

summary(TEM_ann)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
## -7.4130 -0.6821  1.1490  1.7270  3.8680 13.1000    2745

5 Transfer monthly climate data to season (preseason, growseason,posseason)

Data info for each season

Season	months
preseason	2,3,4
growseason	5,6,7,8,9
postseason	8,9

# Function for transfering monthly grid (array) to annual, by selecting "method"= "mean" or "sum"
f_mon2season_grid<-function(index_season,da,method="total") {
  da_season<-array(0,c(nrow(da),ncol(da),dim(da)[3]/12))
  
    for (y in 1:(dim(da)[3]/12)){
      linshi<-matrix(0,ncol=ncol(da),nrow=nrow(da))
      
      for (m in index_season){
        index_sea<- m+(y-1)*12
        linshi<-linshi+da[,,index_sea]
        #print(range(linshi,na.rm=T))
      }
     
      da_season[,,y]<-linshi
     # print(y)
    }
    range(da_season,na.rm=T)
    
  if (method=="Mean" |method=="MEAN" |method=="mean"){
    da_season<-da_season/length(index_season)
  }
    da_season
}

index_season_ZGN<-list("preseason"=c(2,3,4),"growseason"=c(5:9),"postseason"=c(8,9))


PRE_season<-lapply(index_season_ZGN,f_mon2season_grid,da=PRE_mon)
print("the summary info of seasonal (PRE-,Grow-,Post-) rainfall")

## [1] "the summary info of seasonal (PRE-,Grow-,Post-) rainfall"

summary(PRE_season$preseason)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##     0.0   103.5   121.4   117.2   136.8   224.2    2745

summary(PRE_season$growseason)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##     0.0   754.4   840.8   812.4   920.8  1390.0    2745

summary(PRE_season$postseason)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##     0.0   287.4   328.8   323.7   379.4   567.8    2745

TEM_season<-lapply(index_season_ZGN,f_mon2season_grid,da=TEM_mon,method = "mean")
print("the summary info of seasonal (PRE-,Grow-,Post-) temperature")

## [1] "the summary info of seasonal (PRE-,Grow-,Post-) temperature"

summary(TEM_season$preseason)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
## -11.870  -3.691  -1.358  -1.136   1.006  11.190    2745

summary(TEM_season$growseason)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##  -3.593   3.992   6.505   6.863   9.518  20.000    2745

summary(TEM_season$postseason)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##  -3.523   4.292   6.817   7.124   9.761  20.080    2745

6 Staistic all variable for each vegetation type and plot

6.1 sta mean value of Phe. for each grid

# sta mean value of Phe. for each grid
a<-melt(Data_ZGN,id=c(1,2))
Data_ZGN_mean<-dcast(a,ID~variable,mean,na.rm=T)
Data_ZGN_mean$VEG<-factor(Data_ZGN_mean$VEG,levels=c(1,2,3),labels = c("Planted","Mixed","Natural"))

#plot(raster(matrix(Data_ZGN_mean$DEM,nrow = 294)))
Data_ZGN_mean_NONA<-na.omit(Data_ZGN_mean)

ggplot(data=Data_ZGN_mean_NONA,aes(x=DEM,y=SGS,color=VEG))+ geom_point(na.rm = T,alpha=0.05)+facet_wrap(~VEG,nrow = 3) #+geom_smooth(method=lm,se=FALSE) #+ggtitle("SGS")

ggplot(data=Data_ZGN_mean_NONA,aes(x=DEM,y=GSL,color=VEG))+ geom_point(na.rm = T,alpha=0.05)+facet_wrap(~VEG,nrow = 3) #+geom_smooth(method=lm,se=FALSE) #+ggtitle("SGS")

ggplot(data=Data_ZGN_mean_NONA,aes(x=DEM,y=END,color=VEG))+ geom_point(na.rm = T,alpha=0.05)+facet_wrap(~VEG,nrow = 3) #+stat_smooth(method=lm,se=FALSE) #+ggtitle("SGS")

6.2 sta mean time series value of Phe. for each vegetation type

a<-melt(Data_ZGN,id=c(1,2,4))
Data_ZGN_veg_mean<-dcast(a,Year+VEG~variable,mean,na.rm=T)
Data_ZGN_veg_mean$VEG<-factor(Data_ZGN_veg_mean$VEG,levels=c(1,2,3),labels = c("Planted","Mixed","Natural"))

Data_ZGN_veg_mean_NONA<-na.omit(Data_ZGN_veg_mean)

ggplot(data=Data_ZGN_veg_mean_NONA,aes(x=Year,y=SGS,color=VEG))+ geom_line()

ggplot(data=Data_ZGN_veg_mean_NONA,aes(x=Year,y=GSL,color=VEG))+ geom_line()

ggplot(data=Data_ZGN_veg_mean_NONA,aes(x=Year,y=END,color=VEG))+ geom_line()

7 Grid cor analysis of climate and Phe.

7.1 Mean annual relationship

# Function for calculating the correlationship
f_cor<-function(da) {
    y_e<-length(da)
    x_e<-y_e/2
    y_s<-x_e+1
    x<-da[1:x_e];y<-da[y_s:y_e]

    res <- try(cor.test(x,y), silent=TRUE) 
    
    if (class(res)=="try-error") {
     res <- setNames(c(NA, NA), c("estimate","p.value"))
            
    }else{
      .res<-unclass(res)[c("estimate","p.value")]
      res<-unlist(.res)
    }
    res 
}

# function for calculating cor of two array data in gird scale
f_cor_var<-function(da1,da2,Data_ZGN_mean){
  da_sta<-abind(da1,da2)
  dim(da_sta)<-c(dim(da_sta)[1]*dim(da_sta)[2],dim(da_sta)[3])
  da_sta<-t(da_sta)
  #plot(raster(matrix(da_sta[,1],nrow = 294)))
  
  a<-parCapply(cl,da_sta,f_cor)
  Cor_result<-matrix(a,nrow=2)
  op <- par(no.readonly = TRUE) # the whole list of settable par's.

  par(mfrow=c(1,2),mar=c(2.5, 3, 2, 5))
  # plot R2
    ind<-seq(-1,1,by=0.2)
    arg <- list(at=ind, labels=as.character(ind))
    breakpoints<-ind
  #we will select the first 4 colors in the Set1 palette
    n<-length(breakpoints)-1
    color<-brewer.pal(n=n,name="RdYlGn")
  
    .a<-matrix(Cor_result[1,],ncol=ncols,nrow=nrows)
    .a[VEG==0]<-NA
      plot(raster(.a,xmn=102.5809, xmx=103.2294, ymn=31.1837, ymx=31.9289),xlab=expression(paste("Latitude (",degree,")")),ylab=expression(paste("Longitude (",degree,")")),breaks=breakpoints,col=color, axis.arg=arg,legend=T,las=1,cex.lab=1.2,cex.axis=0.8,font.lab=2,horizontal=F,legend.width=1.5,legend.shrink = 0.9,main="R2")  
  #text(120,-40,"R all month",cex=1.1,font=2)
  
  # plot P
  arg <- list(at=c(0,0.01,0.05,0.1,1), labels=c("0","0.01","0.05","0.1","1")) #these are the class names
    breakpoints<-c(0,0.01,0.05,0.1,1)
  color=c("red","orange","skyblue","darkolivegreen1") #and color representation
    .a<-matrix(Cor_result[2,],ncol=ncols,nrow=nrows)
    .a[VEG==0]<-NA
    plot(raster(.a,xmn=102.5809, xmx=103.2294, ymn=31.1837, ymx=31.9289),xlab=expression(paste("Latitude (",degree,")")),ylab=expression(paste("Longitude (",degree,")")),breaks=breakpoints,col=color, axis.arg=arg,legend=T,las=1,cex.lab=1.2,cex.axis=0.8,font.lab=2,horizontal=F,legend.width=1.5,legend.shrink = 0.9,main="P")  
  
    par(op)
    
# plot DEM and R2 for each vegetation type
  print("plot DEM and R2")
    
  Data_ZGN_COR<-Data_ZGN_mean[1:3]
  Data_ZGN_COR$R_PRE_SGS<-Cor_result[1,]

  #plot(raster(matrix(Data_ZGN_mean$DEM,nrow = 294)))
  Data_ZGN_COR_NONA<-na.omit(Data_ZGN_COR)
  
  # set ratio
  ratio.display <- 4/3
  ratio.values <- (max(Data_ZGN_COR_NONA$DEM)-min(Data_ZGN_COR_NONA$DEM))/(max(Data_ZGN_COR_NONA$R_PRE_SGS)-min(Data_ZGN_COR_NONA$R_PRE_SGS))
 
  a<-ggplot(data=Data_ZGN_COR_NONA,aes(x=DEM,y=R_PRE_SGS,color=VEG))+ geom_point(na.rm = T,alpha=0.05)+facet_wrap(~VEG,nrow = 3) + #+geom_smooth(method=lm,se=FALSE) #+ggtitle("SGS")
#  coord_fixed(ratio.values / ratio.display)
  coord_fixed(ratio=400)
  print(a)
  return(Cor_result)
  
}

f_Parallel_set("d")

## [1] "using Zeus for processing"
## 573 MB
## [1] 16

COR_PRE_SGS<-f_cor_var(PRE_ann[,,1:15],SGS,Data_ZGN_mean)

## [1] "plot DEM and R2"

COR_TEM_SGS<-f_cor_var(TEM_ann[,,1:15],SGS,Data_ZGN_mean)

## [1] "plot DEM and R2"

COR_PRE_GSL<-f_cor_var(PRE_ann[,,1:15],GSL,Data_ZGN_mean)

## [1] "plot DEM and R2"

COR_TEM_GSL<-f_cor_var(TEM_ann[,,1:15],GSL,Data_ZGN_mean)

## [1] "plot DEM and R2"

COR_PRE_END<-f_cor_var(PRE_ann[,,1:15],END,Data_ZGN_mean)

## [1] "plot DEM and R2"

COR_TEM_END<-f_cor_var(TEM_ann[,,1:15],END,Data_ZGN_mean)

## [1] "plot DEM and R2"

a<-t(rbind(COR_PRE_SGS,COR_TEM_SGS,COR_PRE_GSL,COR_TEM_GSL,COR_PRE_END,COR_TEM_END))
dim(a)<-c(294,256,12)
write.ENVI(a,"COR_Zagunao")

 ind<-seq(-1,1,by=0.2)
    arg <- list(at=ind, labels=as.character(ind))
    breakpoints<-ind
  #we will select the first 4 colors in the Set1 palette
    n<-length(breakpoints)-1
    color<-brewer.pal(n=n,name="RdYlGn")
    a[VEG==0]<-NA
#plot(raster(a[,,1],xmn=102.5809, xmx=103.2294, ymn=31.1837, ymx=31.9289),xlab=expression(paste("Latitude (",degree,")")),ylab=expression(paste("Longitude (",degree,")")),breaks=breakpoints,col=color, axis.arg=arg,legend=T,las=1,cex.lab=1.2,cex.axis=0.8,font.lab=2,horizontal=F,legend.width=1.5,legend.shrink = 0.9,main="R2")

7.2 Seasonal relationship

f_Parallel_set("d")

## [1] "using Zeus for processing"

## Warning: closing unused connection 19 (<-localhost:11458)

## Warning: closing unused connection 18 (<-localhost:11458)

## Warning: closing unused connection 17 (<-localhost:11458)

## Warning: closing unused connection 16 (<-localhost:11458)

## Warning: closing unused connection 15 (<-localhost:11458)

## Warning: closing unused connection 14 (<-localhost:11458)

## Warning: closing unused connection 13 (<-localhost:11458)

## Warning: closing unused connection 12 (<-localhost:11458)

## Warning: closing unused connection 11 (<-localhost:11458)

## Warning: closing unused connection 10 (<-localhost:11458)

## Warning: closing unused connection 9 (<-localhost:11458)

## Warning: closing unused connection 8 (<-localhost:11458)

## Warning: closing unused connection 7 (<-localhost:11458)

## Warning: closing unused connection 6 (<-localhost:11458)

## Warning: closing unused connection 5 (<-localhost:11458)

## 511 MB
## [1] 16

COR_PRE_pre_SGS<-f_cor_var(PRE_season$preseason[,,1:15],SGS,Data_ZGN_mean)

## [1] "plot DEM and R2"

COR_TEM_pre_SGS<-f_cor_var(TEM_season$preseason[,,1:15],SGS,Data_ZGN_mean)

## [1] "plot DEM and R2"

COR_PRE_grow_GSL<-f_cor_var(PRE_season$growseason[,,1:15],GSL,Data_ZGN_mean)

## [1] "plot DEM and R2"

COR_TEM_grow_GSL<-f_cor_var(TEM_season$growseason[,,1:15],GSL,Data_ZGN_mean)

## [1] "plot DEM and R2"

COR_PRE_post_END<-f_cor_var(PRE_season$postseason[,,1:15],END,Data_ZGN_mean)

## [1] "plot DEM and R2"

COR_TEM_post_END<-f_cor_var(TEM_season$postseason[,,1:15],END,Data_ZGN_mean)

## [1] "plot DEM and R2"

a<-t(rbind(COR_PRE_pre_SGS,COR_TEM_pre_SGS,COR_PRE_grow_GSL,COR_TEM_grow_GSL,COR_PRE_post_END,COR_TEM_post_END))
dim(a)<-c(294,256,12)
write.ENVI(a,"COR_Zagunao_seasonal")

Zagunao Phenology

Ning Liu

24/01/2017