################################################
################################################
################################################
################################################
### LOAD and PROCEED CLIMATIC DATA
################################################
################################################
################################################
################################################
source("./R/FUN.climate.R")
## Data extracted from C. Piedallu data base with
## Radition with cloud cover for each plots
## Piedallu, C, and J Gegout. 2008. “Efficient Assessment of Topographic Solar Radiation to Improve Plant Distribution Models.” Agricultural and Forest Meteorology 148 (11) (October): 1696–1706. doi:10.1016/j.agrformet.2008.06.001.
## Temperature and precipitation with temperature corrected by elevation based on
## Piedallu, C., J. C. Gégout, V. Perez, F. Lebourgeois, and R. Field. 2012. “Soil Water Balance Performs Better Than Climatic Water Variables in Tree Species Distribution Modelling.” Global Ecology and Biogeography. http://onlinelibrary.wiley.com/doi/10.1111/geb.12012/full.

######
## LOAD PIDEALLU DATA
data.CLIM <- read.csv("./data/raw/DataFrance/climate_piedallu/placettesGK_avec_2011.csv",
                      header=T,sep='\t',stringsAsFactors=FALSE,dec=",",na.strings="")

data.CLIM$rumkg_500 <- as.numeric(gsub(",",".",data.CLIM$rumkg_500))


#########################################
#########################################
#### COMPUTE CLIMATIC VARIABLES !!!

### NEED TO COMPUTE DDG5
### NEED TO COMPUTE SOIL MAX WATER CONTENT
### NEED TO COMPUTE PET
### NEED TO COMPUTE WATER BUDGET

## remove NA lines
data.CLIM.na <- is.na(data.CLIM[["tmoy6190_1_cor"]])
data.CLIM2 <-  data.CLIM[!data.CLIM.na,]

### apply function sgdd
sgdd.vec <- apply(data.CLIM2[,82:93],MARGIN=1,FUN=fun.sgdd)
MeanT.vec <- apply(data.CLIM2[,82:93],MARGIN=1,FUN=mean)
plot(MeanT.vec,sgdd.vec)
dim(data.CLIM)

#################################################
### Max Soil Water Content
### compute based on
## Piedallu, C., J. C. Gégout, A. Bruand, and I. Seynave. 2011. “Mapping Soil Water Holding Capacity over Large Areas to Predict Potential Production of Forest Stands.” Geoderma 160 (3): 355–366.

### read data texture as in Piedallu et al. 2011
data.texture <- read.table("./data/raw/DataFrance/climate_piedallu/texture.txt",header=T)
### load ecological data
load(file="./data/process/ecologie_tot.Rdata")
head(ecologie_tot)

codeprof <- c(2.5,10,20,30,40,50,60,70,80,92.5) # code prof en cm
codecaillou <- c(2.5,10,20,30,40,50,60,70,80,90,97.5) ## code percentage en %

#compute swhc
swhc <- fun.swhc(affroc=ecologie_tot$affroc,cailloux=ecologie_tot$cailloux,text2=ecologie_tot$text2,text1=ecologie_tot$text1,prof2=ecologie_tot$prof2,prof1=ecologie_tot$prof1,codeprof,codecaillou,data.texture)
swhc[is.na(ecologie_tot$prof2)] <-  NA

### add to ecological data based 
ecologie_tot$SWHC <- swhc


##################################
##################################
#### COMPUTE PET WITH TURC FORMULA
### unit to convert radneb61_ from J/cm2/month in MJ/m2/day /100/30 then to convert in KJ/m2/day as in Nick formula *1000

## RAD and Temperature of the 12 months 
radneb61_1to12 <- 66:77
tmoy6190_1_cor.1to12 <- 82:93

### apply in parallel
library(doParallel)
registerDoParallel(cores=6) ## affect automaticaly half of the core detected to the foreach
getDoParWorkers() ## here 8 core so 4 core if want to use more registerDoParallel(cores=6)

PET.matrix <-
foreach(i=1:length(data.CLIM$idp), .combine=rbind) %dopar%
	{
fun.PET(i,rad=data.CLIM[,radneb61_1to12],temp=data.CLIM[,tmoy6190_1_cor.1to12])
	}

PET.matrix <- PET.matrix[!data.CLIM.na,]
PET.matrix[PET.matrix<0] <- 0 ## affect zero if negative PET

## plot to check PET computed by me vs PET of Christian
## par(mfrow=c(3,4))
## for (i in 1:12) {plot(PET.matrix[,i],data.CLIM2[,i+48]) ; lines(0:100,0:100)}

colnames(PET.matrix) <-  paste("PET.cor.",1:12,sep="")
##############
### MERGE CLIMATE and PET
data.CLIM2 <- cbind(data.CLIM2,PET.matrix,sgdd=sgdd.vec)

### MERGE WITH ECOLOGICAL DATA
ecologie.clim <- merge(ecologie_tot,data.CLIM2,by="idp",all.x=T)
dim(ecologie_tot)
dim(ecologie.clim)
names(ecologie.clim)


ecologie.clim2 <-  ecologie.clim[!is.na(ecologie.clim$PET.cor.1) &!is.na(ecologie.clim$SWHC),]

###########################################################################################
###### COMPUTE WATER BUDGET

## ### test function
## fun.WaterBudget(i=16,prcp.m=(ecologie.clim2[,paste("prec6190_",1:12,sep="")]),
##                 PET.m=(ecologie.clim2[,paste("PET.cor.",1:12,sep="")]),
##                 Ta.m=(ecologie.clim2[,c("tmoy6190_1_cor",paste("tmoy6190_1_cor.",1:11,sep=""))]),
##                 SWHC.v=(ecologie.clim2[  ,"SWHC"]),n=2)


### apply function in parallel
library(doParallel)
registerDoParallel(cores=6) ## affect automaticaly half of the core detected to the foreach
getDoParWorkers() ## here 8 core so 4 core if want to use more registerDoParallel(cores=6)

WBWS.matrix <-
foreach(i=1:length(ecologie.clim2$idp), .combine=rbind) %dopar%
	{
fun.WaterBudget(i,prcp.m=(ecologie.clim2[,paste("prec6190_",1:12,sep="")]),
                PET.m=(ecologie.clim2[,paste("PET.cor.",1:12,sep="")]),
                Ta.m=(ecologie.clim2[,c("tmoy6190_1_cor",paste("tmoy6190_1_cor.",1:11,sep=""))]),
                SWHC.v=(ecologie.clim2[  ,"SWHC"]),n=2)
         }

### MERGE all and saved
WBWS.matrix2 <- cbind("idp"=ecologie.clim2$idp,WBWS.matrix)
ecologie.clim.data <- merge(ecologie.clim,WBWS.matrix2,by="idp",all.x=T)

#####
# change tplant not good format
load("./data/process/placette_tot.Rdata")
ecologie.clim.data$tplant <-  placette_tot$tplant

### Compute mean T and annual sum of precip

ecologie.clim.data$MAT <- apply(ecologie.clim.data[,114:125],MARGIN=1,FUN=mean)
ecologie.clim.data$SAP <- apply(ecologie.clim.data[,69:80],MARGIN=1,FUN=sum)


saveRDS(ecologie.clim.data,file="./data/process/ecologie.clim.data.rds")






#########################################################
#########################################################
###### FIGURES OF CLIMATIC DATA
#############################
ecologie.clim.data <- readRDS("./data/process/ecologie.clim.data.rds")
 names(ecologie.clim.data)

###check  climatic data

pdf("./figs/sgdd.tmin.map.pdf")
par(mfrow=c(2,2))
plot(ecologie.clim.data$xl93,ecologie.clim.data$yl93,
     col=rev(heat.colors(10))[cut(ecologie.clim.data$sgdd,quantile(ecologie.clim.data$sgdd,probs=(0:10)/10,na.rm=T),labels=F)],
     cex=0.2,pty="s",xlab=NA,ylab=NA,main="SGDD")

plot(ecologie.clim.data$xl93,ecologie.clim.data$yl93,
 col=rev(heat.colors(10))[cut(ecologie.clim.data$tmin6190_min_cor,quantile(ecologie.clim.data$tmin6190_min_cor,probs=(0:10)/10,na.rm=T),labels=F)],
     cex=0.2,pty="s",xlab=NA,ylab=NA,main="Tmin")
 plot( ecologie.clim.data$tmin6190_min_cor,ecologie.clim.data$sgdd,cex=0.2,xlab="Tmin",ylab="SGDD")
dev.off()

pdf("./figs/Water.map.pdf")
  par(mfrow=c(2,2),mar=c(1,1,1,1))
plot(ecologie.clim.data$xl93,ecologie.clim.data$yl93,
 col=rev(colorRampPalette(c("blue", "red"))( 10 ))[cut(ecologie.clim.data$WB.y,quantile(ecologie.clim.data$WB.y,probs=(0:10)/10,na.rm=T),labels=F)],
     cex=0.2,pty="s",xlab=NA,ylab=NA)
plot(ecologie.clim.data$xl93,ecologie.clim.data$yl93,
  col=rev(colorRampPalette(c("blue", "red"))( 10 ))[cut(ecologie.clim.data$WB.s,quantile(ecologie.clim.data$WB.s,probs=(0:10)/10,na.rm=T),labels=F)],
     cex=0.2,pty="s",xlab=NA,ylab=NA)
plot(ecologie.clim.data$xl93,ecologie.clim.data$yl93,
  col=rev(colorRampPalette(c("blue", "red"))( 10 ))[cut(ecologie.clim.data$WS.y,quantile(ecologie.clim.data$WS.y,probs=(0:10)/10,na.rm=T),labels=F)],
     cex=0.2,pty="s",xlab=NA,ylab=NA)
plot(ecologie.clim.data$xl93,ecologie.clim.data$yl93,
 col=rev(colorRampPalette(c("blue", "red"))( 10 ))[cut(ecologie.clim.data$WS.s,quantile(ecologie.clim.data$WS.s,probs=(0:10)/10,na.rm=T),labels=F)],
     cex=0.2,pty="s",xlab=NA,ylab=NA)
dev.off()

pdf("./figs/Water.var,cor.pdf")
par(mfrow=c(2,2))
plot(ecologie.clim.data[['WB.y']],ecologie.clim.data[['WS.y']],xlab="WB year",ylab="ratio AET/D year",cex=0.2)
plot(ecologie.clim.data[['WB.s']],ecologie.clim.data[['WS.s']],xlab="WB growing season",ylab="ratio AET/D growing season",cex=0.2)
plot(ecologie.clim.data[['WB.s']],ecologie.clim.data[['WB.y']],xlab="WB growing season",ylab="WB growing year",cex=0.2)
plot(ecologie.clim.data[['WS.s']],ecologie.clim.data[['WS.y']],xlab="ratio AET/D growing season",ylab="ratio AET/D year",cex=0.2)
dev.off()
pdf("./figs/sgdd.water.cor.pdf")
plot(ecologie.clim.data[['sgdd']],ecologie.clim.data[['WS.y']],xlab="SGDD",ylab="ratio AET/D year",cex=0.2)
dev.off()

### DO PCA OF CLIMATIC DATA

data.prc <- prcomp(ecologie.clim.data[apply(is.na(ecologie.clim.data[,c('sgdd','tmin6190_min_cor','WB.s','WS.s')]),MARGIN=1,FUN=sum)<1,c('sgdd','tmin6190_min_cor','WB.s','WS.s')],na.action=na.omit,scale=T)
summary(data.prc)
scores <-data.prc$x
pdf("./figs/climate.pca.pdf")
biplot(data.prc,pch=1,xlabs=rep('.',length(data.prc$x[,1])),ylabs=c('SGDD','Tmin','WB','WS'))
 text(150,180,labels=paste("Var PC1 ", ((summary(data.prc)))$importance[2,1],sep=""))
 text(150,165,labels=paste("Var PC2 ", ((summary(data.prc)))$importance[2,2],sep=""))
dev.off()


##### FIGS on elevation correction of temperature done by Christian
### map of diff between mean altitude of the 1x1km cell and actual elevation of the plots
pdf("./figs/map.error.alti.pdf")
plot(data.CLIM$xl2,data.CLIM$yl2,col="grey",cex=0.3,xlab="X",ylab="Y")
data.temp <- data.CLIM[data.CLIM$diff..altIFN...mnt50 < (-100) | (data.CLIM$diff..altIFN...mnt50>100),]
points(data.temp$xl2,data.temp$yl2,col=c("red","coral3","coral","grey","coral","coral3","red")[cut(data.temp$diff..altIFN...mnt50,breaks=c(-500,-300,-200,-100,100,200,300,500),labels=F)],cex=0.5,pch=1)
dev.off()

##### PLOT HIST OF ALL VARIABLES
pdf('./figs/climatvar.hist.pdf')
par(mfrow=c(3,4))
for (i in 37:48) hist(data.CLIM[,i],xlab='Monthly Precip (mm)',main=names(data.CLIM)[i])

par(mfrow=c(3,4))
for (i in 49:60) hist(data.CLIM[,i],xlab='PET (mm)',main=names(data.CLIM)[i])

par(mfrow=c(3,4))
for (i in 66:77) hist(data.CLIM[,i],xlab='Radiation (?)',main=names(data.CLIM)[i])

par(mfrow=c(3,4))
for (i in 82:93) hist(data.CLIM[,i],xlab='temperature (C)',main=names(data.CLIM)[i])
dev.off()

pdf('./figs/xyplot.Tunco.Tcor.pdf')
par(mfrow=c(3,4))
for (i in 1:12) plot(data.CLIM[,(23:34)[i]],data.CLIM[,(82:94)[i]],xlab='temperature (C) uncorrected',ylab='temperature (C) corrected',
                     main=names(data.CLIM)[(23:34)[i]])
dev.off()