Zagunao Phenology
Ning Liu
24/01/2017
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")
TEM_mon<-read.ENVI("ZGNPhen/TEM")
PRE_mon<-read.ENVI("ZGNPhen/PRE")
nrows<-nrow(DEM)
ncols<-ncol(DEM)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 :3396754 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] 11summary(PRE_ann)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0 941.4 1035.0 986.4 1125.0 1507.0TEM_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] 11summary(TEM_ann)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -7.3550 -0.5007 1.1770 1.8130 3.9160 13.10005 Staistic all variable for each vegetation type and plot
5.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")
5.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()
6 Grid cor analysis of climate and Phe.
# 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)
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")
# 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)
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")
print(a)
return(Cor_result)
}
f_Parallel_set("d")## [1] "using Zeus for processing"
## 435 MB
## [1] 16COR_PRE_SGS<-f_cor_var(PRE_ann,SGS[,,1:11],Data_ZGN_mean)
## [1] "plot DEM and R2"
COR_TEM_SGS<-f_cor_var(TEM_ann,SGS[,,1:11],Data_ZGN_mean)
## [1] "plot DEM and R2"
COR_PRE_GSL<-f_cor_var(PRE_ann,GSL[,,1:11],Data_ZGN_mean)
## [1] "plot DEM and R2"
COR_TEM_GSL<-f_cor_var(TEM_ann,GSL[,,1:11],Data_ZGN_mean)
## [1] "plot DEM and R2"
COR_PRE_END<-f_cor_var(PRE_ann,END[,,1:11],Data_ZGN_mean)
## [1] "plot DEM and R2"
COR_TEM_END<-f_cor_var(TEM_ann,END[,,1:11],Data_ZGN_mean)
## [1] "plot DEM and R2"