diff --git a/R/CLIMATE.FRANCE.R b/R/CLIMATE.FRANCE.R
new file mode 100644
index 0000000000000000000000000000000000000000..7e046dee7e11515feb25ddd45b2b7d2ba133e874
--- /dev/null
+++ b/R/CLIMATE.FRANCE.R
@@ -0,0 +1,243 @@
+################################################
+################################################
+################################################
+################################################
+### 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="")
+print(names(data.CLIM))
+
+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()
+
+
diff --git a/R/READ.DATA.NFI.FRANCE.R b/R/READ.DATA.NFI.FRANCE.R
new file mode 100644
index 0000000000000000000000000000000000000000..bb84aa0aa4e17119d65923f6a38d755a57ada4d7
--- /dev/null
+++ b/R/READ.DATA.NFI.FRANCE.R
@@ -0,0 +1,443 @@
+############################################################
+############################################################
+########## READ, MERGE AND CLEAN ALL NFI DATA NEW METHODS
+
+### read TREE table downloaded from the web
+arbre2005 <- read.csv("./data/raw/DataFrance/2005/arbres_foret_2005.csv",sep=";", stringsAsFactors=FALSE)
+summary(arbre2005)
+arbre2006 <- read.csv("./data/raw/DataFrance/2006/arbres_foret_2006.csv",sep=";", stringsAsFactors=FALSE)
+summary(arbre2006)
+arbre2007 <- read.csv("./data/raw/DataFrance/2007/arbres_foret_2007.csv",sep=";", stringsAsFactors=FALSE)
+summary(arbre2007)
+arbre2008 <- read.csv("./data/raw/DataFrance/2008/arbres_foret_2008.csv",sep=";", stringsAsFactors=FALSE)
+summary(arbre2005)
+arbre2009 <- read.csv("./data/raw/DataFrance/2009/arbres_foret_2009.csv",sep=";", stringsAsFactors=FALSE)
+summary(arbre2009)
+arbre2010 <- read.csv("./data/raw/DataFrance/2010/arbres_foret_2010.csv",sep=";", stringsAsFactors=FALSE)
+summary(arbre2010)
+arbre2011 <- read.csv("./data/raw/DataFrance/2011/arbres_foret_2011.csv",sep=";", stringsAsFactors=FALSE)
+
+
+
+
+
+#### ALL TREE after 2007 don't have the veget value as in 2005 aned 2006 because new variable ACCI for tree with accident (trunk broken ...)
+### NEED TO UPDATE VEGET FROM ACCI
+arbre2007$veget[arbre2007$acci>0] <- 1
+arbre2008$veget[arbre2008$acci>0] <- 1
+arbre2009$veget[arbre2009$acci>0] <- 1
+arbre2010$veget[arbre2010$acci>0] <- 1
+arbre2011$veget[arbre2011$acci>0] <- 1
+##
+arbre2005$veget <- unclass(arbre2005$veget)-1
+arbre2006$veget <- unclass(arbre2006$veget)-1
+
+######
+## ORI NEED TO BE DEGRADED SINCE 2007 to two level from resprout or from seed
+arbre2007$ori[arbre2007$ori==2] <- 0
+arbre2008$ori[arbre2008$ori==2] <- 0
+arbre2009$ori[arbre2009$ori==2] <- 0
+arbre2010$ori[arbre2010$ori==2] <- 0
+arbre2011$ori[arbre2011$ori==2] <- 0
+
+
+##############################
+### merge all table adding NA when no variable for that year
+arbre.tot <- data.frame(idp=c(arbre2005$idp,arbre2006$idp,arbre2007$idp,arbre2008$idp,arbre2009$idp,arbre2010$idp,arbre2011$idp),
+ a=c(arbre2005$a,arbre2006$a,arbre2007$a,arbre2008$a,arbre2009$a,arbre2010$a,arbre2011$a),
+ veget=c(arbre2005$veget,arbre2006$veget,arbre2007$veget,arbre2008$veget,arbre2009$veget,arbre2010$veget,arbre2011$veget),
+ simplif=c(rep(NA,length(arbre2005$idp)),rep(NA,length(arbre2006$idp)),rep(NA,length(arbre2007$idp)),
+ rep(NA,length(arbre2008$idp)),arbre2009$simplif,arbre2010$simplif,arbre2011$simplif),
+ acci=c(rep(NA,length(arbre2005$idp)),rep(NA,length(arbre2006$idp)),arbre2007$acci,arbre2008$acci,arbre2009$acci,
+ arbre2010$acci,arbre2011$acci),
+ espar=c(as.character(arbre2005$espar),as.character(arbre2006$espar),as.character(arbre2007$espar),
+ as.character(arbre2008$espar),as.character(arbre2009$espar),as.character(arbre2010$espar),as.character(arbre2011$espar)),
+ ori=c(arbre2005$ori,arbre2006$ori,arbre2007$ori,arbre2008$ori,arbre2009$ori,arbre2010$ori,arbre2011$ori),
+ lib=c(arbre2005$lib,arbre2006$lib,arbre2007$lib,arbre2008$lib,arbre2009$lib,arbre2010$lib,arbre2011$lib),
+ forme=c(arbre2005$forme,arbre2006$forme,arbre2007$forme,arbre2008$forme,arbre2009$forme,arbre2010$forme,arbre2011$forme),
+ tige=c(arbre2005$tige,arbre2006$tige,arbre2007$tige,arbre2008$tige,arbre2009$tige,arbre2010$tige,arbre2011$tige),
+ mortb=c(rep(NA,length(arbre2005$idp)),arbre2006$mortb,arbre2007$mortb,arbre2008$mortb,arbre2009$mortb,arbre2010$mortb,
+ arbre2011$mortb),
+ sfgui=c(rep(NA,length(arbre2005$idp)),rep(NA,length(arbre2006$idp)),rep(NA,length(arbre2007$idp)),arbre2008$sfgui,
+ arbre2009$sfgui,arbre2010$sfgui,arbre2011$sfgui),
+ sfgeliv=c(rep(NA,length(arbre2005$idp)),rep(NA,length(arbre2006$idp)),rep(NA,length(arbre2007$idp)),arbre2008$sfgeliv,
+ arbre2009$sfgeliv,arbre2010$sfgeliv,arbre2011$sfgeliv),
+ sfpied=c(rep(NA,length(arbre2005$idp)),rep(NA,length(arbre2006$idp)),rep(NA,length(arbre2007$idp)),arbre2008$sfpied,
+ arbre2009$sfpied,arbre2010$sfpied,arbre2011$sfpied),
+ sfdorge=c(rep(NA,length(arbre2005$idp)),rep(NA,length(arbre2006$idp)),rep(NA,length(arbre2007$idp)),arbre2008$sfdorge,
+ arbre2009$sfdorge,arbre2010$sfdorge,arbre2011$sfdorge),
+ sfcoeur=c(rep(NA,length(arbre2005$idp)),rep(NA,length(arbre2006$idp)),rep(NA,length(arbre2007$idp)),
+ rep(NA,length(arbre2008$idp)),arbre2009$sfcoeur,arbre2010$sfcoeur,arbre2011$sfcoeur),
+ c13=c(arbre2005$c13,arbre2006$c13,arbre2007$c13,arbre2008$c13,arbre2009$c13,arbre2010$c13,arbre2011$c13),
+ ir5=c(arbre2005$ir5,arbre2006$ir5,arbre2007$ir5,arbre2008$ir5,arbre2009$ir5,arbre2010$ir5,arbre2011$ir5),
+ htot=c(arbre2005$htot,arbre2006$htot,arbre2007$htot,arbre2008$htot,arbre2009$htot,arbre2010$htot,arbre2011$htot),
+ hdec=c(rep(NA,length(arbre2005$idp)),rep(NA,length(arbre2006$idp)),rep(NA,length(arbre2007$idp)),arbre2008$hdec,
+ arbre2009$hdec,arbre2010$hdec,arbre2011$hdec),
+ decoupe=c(rep(NA,length(arbre2005$idp)),rep(NA,length(arbre2006$idp)),rep(NA,length(arbre2007$idp)),arbre2008$decoupe,
+ arbre2009$decoupe,arbre2010$decoupe,arbre2011$decoupe),
+ q1=c(arbre2005$q1,arbre2006$q1,arbre2007$q1,arbre2008$q1,arbre2009$q1,arbre2010$q1,arbre2011$q1),
+ q2=c(arbre2005$q2,arbre2006$q2,arbre2007$q2,arbre2008$q2,arbre2009$q2,arbre2010$q2,arbre2011$q2),
+ q3=c(arbre2005$q3,arbre2006$q3,arbre2007$q3,arbre2008$q3,arbre2009$q3,arbre2010$q3,arbre2011$q3),
+ r=c(arbre2005$r,arbre2006$r,arbre2007$r,arbre2008$r,arbre2009$r,arbre2010$r,arbre2011$r),
+ lfsd=c(arbre2005$lfsd,arbre2006$lfsd,arbre2007$lfsd,arbre2008$lfsd,arbre2009$lfsd,arbre2010$lfsd,arbre2011$lfsd),
+ age=c(rep(NA,length(arbre2005$idp)),rep(NA,length(arbre2006$idp)),rep(NA,length(arbre2007$idp)),arbre2008$age,
+ arbre2009$age,arbre2010$age,arbre2011$age),
+ v=c(arbre2005$v,arbre2006$v,arbre2007$v,arbre2008$v,arbre2009$v,arbre2010$v,arbre2011$v),
+ w=c(arbre2005$w,arbre2006$w,arbre2007$w,arbre2008$w,arbre2009$w,arbre2010$w,arbre2011$w),
+ YEAR=c(rep(2005,length(arbre2005$idp)),rep(2006,length(arbre2006$idp)),rep(2007,length(arbre2007$idp)),
+ rep(2008,length(arbre2008$idp)),rep(2009,length(arbre2009$simplif)),rep(2010,length(arbre2010$simplif)),
+ rep(2011,length(arbre2011$simplif))))
+
+rm(arbre2005,arbre2006,arbre2007,arbre2008,arbre2009,arbre2010,arbre2011)
+gc()
+########################################################
+## #### check problem of unit for c13 ir5 and htot by plotting the data
+## plot(arbre.tot$c13,arbre.tot$ir5,col=unclass(factor(arbre.tot$YEAR)),cex=0.1)
+## boxplot(arbre.tot$c13~arbre.tot$YEAR,ylab="c13")
+## boxplot(arbre.tot$ir5~arbre.tot$YEAR,ylab="ir5")
+## boxplot(arbre.tot$htot~arbre.tot$YEAR,ylab="htot")
+## boxplot(arbre.tot$age~arbre.tot$YEAR,ylab="age")
+## ##### SOMETHING VERY STRANGE FOR THE AGE WITH SEVERAL TREE OVER 1000 YEARS OLD
+## boxplot(arbre.tot$mortb~arbre.tot$YEAR,ylab="mortality branche")
+## boxplot(arbre.tot$veget~arbre.tot$YEAR,ylab="accident")
+## boxplot(arbre.tot$simplif~arbre.tot$YEAR,ylab="accident")
+## boxplot(arbre.tot$w~arbre.tot$YEAR,ylab="accident")
+## boxplot(arbre.tot$ori~arbre.tot$YEAR,ylab="accident")
+## boxplot(arbre.tot$lib~arbre.tot$YEAR,ylab="accident")
+## boxplot(arbre.tot$forme~arbre.tot$YEAR,ylab="accident")
+## boxplot(arbre.tot$tige~arbre.tot$YEAR,ylab="accident")
+## ### NEED TO USE ONLY THE TIGE == 1 in the ANALYSIS
+## boxplot(arbre.tot$sfgui~arbre.tot$YEAR,ylab="accident")
+## boxplot(arbre.tot$sfgeliv~arbre.tot$YEAR,ylab="accident")
+
+
+## ##### CHECK OTHER VARIABLE OK DONE
+## ### USE BRANCH MORTALITY AS AN INDICATOR OF MORTALITY ?? ABIOTIC STRESS ?
+
+## x11()
+## plot(arbre.tot$c13,arbre.tot$htot,col=unclass(factor(arbre.tot$YEAR)))
+
+save(arbre.tot,file="./data/process/arbre.tot.Rdata")
+
+
+
+#########################################
+###### DEAD
+#########################################
+
+
+##############################################################
+#############################################################
+## READ AND MERGE DEAD DATA
+### MERGE WITH DEAD TREE and MERGE WITH PLOT DATA!!
+
+### read DEAD TREE table downloaded from the web
+arbre_mort2005 <- read.csv("./data/raw/DataFrance/2005/arbres_morts_foret_2005.csv",sep=";"
+ , stringsAsFactors=FALSE)
+## summary(arbre_mort2005)
+arbre_mort2006 <- read.csv("./data/raw/DataFrance/2006/arbres_morts_foret_2006.csv",sep=";"
+ , stringsAsFactors=FALSE)
+## summary(arbre_mort2006)
+arbre_mort2007 <- read.csv("./data/raw/DataFrance/2007/arbres_morts_foret_2007.csv",sep=";"
+ , stringsAsFactors=FALSE)
+## summary(arbre_mort2007)
+arbre_mort2008 <- read.csv("./data/raw/DataFrance/2008/arbres_morts_foret_2008.csv",sep=";"
+ , stringsAsFactors=FALSE)
+## summary(arbre_mort2005)
+arbre_mort2009 <- read.csv("./data/raw/DataFrance/2009/arbres_morts_foret_2009.csv",sep=";"
+ , stringsAsFactors=FALSE)
+## summary(arbre_mort2009)
+arbre_mort2010 <- read.csv("./data/raw/DataFrance/2010/arbres_morts_foret_2010.csv",sep=";"
+ , stringsAsFactors=FALSE)
+## summary(arbre_mort2010)
+arbre_mort2011 <- read.csv("./data/raw/DataFrance/2011/arbres_morts_foret_2011.csv",sep=";"
+ , stringsAsFactors=FALSE)
+
+## names(arbre_mort2005)
+## names(arbre_mort2006)
+## names(arbre_mort2007)
+## names(arbre_mort2008)
+## names(arbre_mort2009)
+## names(arbre_mort2010)
+## names(arbre_mort2011)
+
+### merge 2005 2006 2007 to compute c13
+
+arbre_mort05_07 <- rbind(arbre_mort2005, arbre_mort2006, arbre_mort2007)
+arbre_mort05_07$c13 <- rep(NA,length(arbre_mort05_07$c0))
+arbre_mort05_07$espar2 <- as.numeric(substr(arbre_mort05_07$espar,1,2))
+arbre_mort05_07$year <- c(rep(2005,length=length(arbre_mort2005[,1])),
+ rep(2006,length=length(arbre_mort2006[,1])),
+ rep(2007,length=length(arbre_mort2007[,1])))
+
+#### NEED TO CONVERT c0 into c13 before 2008
+### before 2008 no date dead but all tree died less than 5 years ago.
+### need to do convertion between c0 and c13.
+## for that use the NFI data from previous inventory that were reporting both c0 and c13 for all species
+## fit an allometric relationship and then use it to predict c13 in this data base
+
+### READ DATA CYCLE 3 ORGINAL /// NEED TO CONVERT C0 and C13 in cm *100
+arbre.cycle3 <- read.table("./data/raw/DataFrance/cycle3/data.arbre.tot.txt",sep=" ", stringsAsFactors=FALSE)
+### change the C from m to cm
+arbre.cycle3$C0 <- arbre.cycle3$C0*100
+arbre.cycle3$C13 <- arbre.cycle3$C13*100
+
+
+## LOAD library RMA regression
+library(lmodel2)
+
+
+## the regression between C13 and C0 vary between species, but not same species in cycle 3 because less details (the classification is only based on number and no letters
+## remove the letters and apply the same model to all species that have the same number code
+species.list <- c(as.numeric(names(table(substr(arbre_mort05_07$espar,1,2))))[-1])
+## length(table(arbre.cycle3$ESS))
+
+for (i in species.list)
+{
+if (sum(arbre.cycle3$ESS==i)>50)
+ {
+lmodel2.res <- lmodel2(C13~C0 ,data=arbre.cycle3[arbre.cycle3$ESS==i,],range.x="relative",range.y="relative")
+arbre_mort05_07$c13[arbre_mort05_07$espar2==i & !is.na(arbre_mort05_07$espar2)] <-
+ lmodel2.res$regression.results[4,2] + arbre_mort05_07$c0[arbre_mort05_07$espar2==i &
+ !is.na(arbre_mort05_07$espar2)]*lmodel2.res$regression.results[4,3]
+print(i)
+print(lmodel2.res$regression.results[4,2:3])
+print(range(arbre_mort05_07$c13[arbre_mort05_07$espar2==i &
+ !is.na(arbre_mort05_07$espar2)],na.rm=T))
+ }
+else
+ {
+lmodel2.res <- lmodel2(C13~C0 ,data=arbre.cycle3,range.x="relative",range.y="relative")
+arbre_mort05_07$c13[arbre_mort05_07$espar2==i & !is.na(arbre_mort05_07$espar2)] <-
+ lmodel2.res$regression.results[4,2] + arbre_mort05_07$c0[arbre_mort05_07$espar2==i &
+ !is.na(arbre_mort05_07$espar2)]*lmodel2.res$regression.results[4,3]
+print(i)
+print(lmodel2.res$regression.results[4,2:3])
+print(range(arbre_mort05_07$c13[arbre_mort05_07$espar2==i & !is.na(arbre_mort05_07$espar2)],na.rm=T))
+
+ }
+}
+
+### for species with NO DATA in cycle 3 apply mean model over all species
+lmodel2.res <- lmodel2(C13~C0 ,data=arbre.cycle3,range.x="relative",range.y="relative")
+arbre_mort05_07$c13[is.na(arbre_mort05_07$espar2)] <-
+ lmodel2.res$regression.results[4,2] + arbre_mort05_07$c0[is.na(arbre_mort05_07$espar2)]*lmodel2.res$regression.results[4,3]
+
+## check predicted C13 from C0
+head(cbind(arbre_mort05_07$c13,arbre_mort05_07$c0,arbre_mort05_07$espar) )
+
+### ok donne c13 added
+
+########################
+## MERGE WITH other dead data
+
+arbre_mort_tot <- data.frame(idp=c(arbre_mort05_07$idp,arbre_mort2008$idp,arbre_mort2009$idp,arbre_mort2010$idp,arbre_mort2011$idp),
+ a=c(arbre_mort05_07$a,arbre_mort2008$a,arbre_mort2009$a,arbre_mort2010$a,arbre_mort2011$a),
+ espar=c(as.character(arbre_mort05_07$espar),as.character(arbre_mort2008$espar),as.character(arbre_mort2009$espar),
+ as.character(arbre_mort2010$espar),as.character(arbre_mort2011$espar)),
+ ori=c(arbre_mort05_07$ori,arbre_mort2008$ori,arbre_mort2009$ori,arbre_mort2010$ori,arbre_mort2011$ori),
+ veget=c(arbre_mort05_07$veget,arbre_mort2008$veget,arbre_mort2009$veget,arbre_mort2010$veget,arbre_mort2011$veget),
+ datemort=c(rep(NA,length(arbre_mort05_07$ori)),arbre_mort2008$datemort,arbre_mort2009$datemort,arbre_mort2010$datemort
+ ,arbre_mort2011$datemort),
+ c13=c(arbre_mort05_07$c13,arbre_mort2008$c13,arbre_mort2009$c13,arbre_mort2010$c13,arbre_mort2011$c13),
+ v=c(arbre_mort05_07$v,arbre_mort2008$v,arbre_mort2009$v,arbre_mort2010$v,arbre_mort2011$v),
+ w=c(arbre_mort05_07$w,arbre_mort2008$w,arbre_mort2009$w,arbre_mort2010$w,arbre_mort2011$w),
+ YEAR=c(rep(2005,length(arbre_mort2005$idp)),rep(2006,length(arbre_mort2006$idp)),rep(2007,length(arbre_mort2007$idp))
+ ,rep(2008,length(arbre_mort2008$idp)),rep(2009,length(arbre_mort2009$idp)),rep(2010,length(arbre_mort2010$idp))
+ ,rep(2011,length(arbre_mort2011$idp)))
+ )
+rm(arbre.cycle3,arbre_mort2005,arbre_mort2006,arbre_mort2007,arbre_mort2008,arbre_mort2009,arbre_mort2010,arbre_mort2011)
+gc()
+
+save(arbre_mort_tot,file="./data/process/arbre_mort_tot.Rdata")
+
+
+
+################################################################################
+#### MERGE DEAD AND ALIVE TREE
+
+
+head(arbre_mort_tot)
+head(arbre.tot)
+
+arbre.ALIVE.DEAD <- data.frame(
+ idp=c(arbre.tot$idp,arbre_mort_tot$idp),
+ a=c(arbre.tot$a,arbre_mort_tot$a),
+ veget=c(arbre.tot$veget,arbre_mort_tot$veget),
+ simplif=c(arbre.tot$simplif,rep(NA,length=length(arbre_mort_tot$idp))),
+ acci=c(arbre.tot$acci,rep(NA,length=length(arbre_mort_tot$idp))),
+ espar=c(as.character(arbre.tot$espar),as.character(arbre_mort_tot$espar)),
+ ori=c(arbre.tot$ori,arbre_mort_tot$ori),
+ lib=c(arbre.tot$lib,rep(NA,length=length(arbre_mort_tot$idp))),
+ forme=c(arbre.tot$forme,rep(NA,length=length(arbre_mort_tot$idp))),
+ tige=c(arbre.tot$tige,rep(NA,length=length(arbre_mort_tot$idp))),
+ mortb=c(arbre.tot$mortb,rep(NA,length=length(arbre_mort_tot$idp))),
+ sfgui=c(arbre.tot$sfgui,rep(NA,length=length(arbre_mort_tot$idp))),
+ sfgeliv=c(arbre.tot$sfgeliv,rep(NA,length=length(arbre_mort_tot$idp))),
+ sfpied=c(arbre.tot$sfpied,rep(NA,length=length(arbre_mort_tot$idp))),
+ sfdorge=c(arbre.tot$sfdorge,rep(NA,length=length(arbre_mort_tot$idp))),
+ sfcoeur=c(arbre.tot$sfcoeur,rep(NA,length=length(arbre_mort_tot$idp))),
+ c13=c(arbre.tot$c13,arbre_mort_tot$c13),
+ ir5=c(arbre.tot$ir5,rep(NA,length=length(arbre_mort_tot$idp))),
+ htot=c(arbre.tot$htot,rep(NA,length=length(arbre_mort_tot$idp))),
+ hdec=c(arbre.tot$hdec,rep(NA,length=length(arbre_mort_tot$idp))),
+ decoupe=c(arbre.tot$decoupe,rep(NA,length=length(arbre_mort_tot$idp))),
+ q1=c(arbre.tot$q1,rep(NA,length=length(arbre_mort_tot$idp))),
+ q2=c(arbre.tot$q2,rep(NA,length=length(arbre_mort_tot$idp))),
+ q3=c(arbre.tot$q3,rep(NA,length=length(arbre_mort_tot$idp))),
+ r=c(arbre.tot$r,rep(NA,length=length(arbre_mort_tot$idp))),
+ lfsd=c(arbre.tot$lfsd,rep(NA,length=length(arbre_mort_tot$idp))),
+ age=c(arbre.tot$age,rep(NA,length=length(arbre_mort_tot$idp))),
+ v=c(arbre.tot$v,arbre_mort_tot$v),
+ w=c(arbre.tot$w,rep(10000/(pi*(c(15))^2),length=length(arbre_mort_tot$w))),## assume that all dead tree are sampled on the whole plot as explained in the method
+ YEAR=c(arbre.tot$YEAR,arbre_mort_tot$YEAR),
+ datemort=c(rep(NA,length=length(arbre.tot$YEAR)),arbre_mort_tot$datemort),
+ dead=c(rep(0,length=length(arbre.tot$YEAR)),rep(1,length=length(arbre_mort_tot$idp))))## 1 = dead
+
+## delete plot with DEAD tree missing because of no C13 or no espar
+arbre.ALIVE.DEAD2 <- arbre.ALIVE.DEAD[!(arbre.ALIVE.DEAD$idp %in% unique(c(names(tapply(is.na(arbre.ALIVE.DEAD$c13),
+ INDEX=arbre.ALIVE.DEAD$idp,FUN=sum))[tapply(is.na(arbre.ALIVE.DEAD$c13),INDEX=arbre.ALIVE.DEAD$idp,FUN=sum)>0],
+ names(tapply(is.na(arbre.ALIVE.DEAD$espar),INDEX=arbre.ALIVE.DEAD$idp,FUN=sum))[tapply(is.na(arbre.ALIVE.DEAD$espar),
+ INDEX=arbre.ALIVE.DEAD$idp,FUN=sum)>0]))),]
+
+save(arbre.ALIVE.DEAD2,file='./data/process/arbre.ALIVE.DEAD2.Rdata')
+
+
+########################################################################################
+#########################################
+#########################################
+#########################################
+##### LOAD DATA FOR PLOT INFO
+### read DEAD TREE table downloaded from the web
+placette2005 <- read.csv("./data/raw/DataFrance/2005/placettes_foret_2005.csv",sep=";"
+ , stringsAsFactors=FALSE)
+## summary(placette2005)
+placette2006 <- read.csv("./data/raw/DataFrance/2006/placettes_foret_2006.csv",sep=";"
+ , stringsAsFactors=FALSE)
+## summary(placette2006)
+placette2007 <- read.csv("./data/raw/DataFrance/2007/placettes_foret_2007.csv",sep=";"
+ , stringsAsFactors=FALSE)
+## summary(placette2007)
+placette2008 <- read.csv("./data/raw/DataFrance/2008/placettes_foret_2008.csv",sep=";"
+ , stringsAsFactors=FALSE)
+## summary(placette2005)
+placette2009 <- read.csv("./data/raw/DataFrance/2009/placettes_foret_2009.csv",sep=";"
+ , stringsAsFactors=FALSE)
+## summary(placette2009)
+placette2010 <- read.csv("./data/raw/DataFrance/2010/placettes_foret_2010.csv",sep=";"
+ , stringsAsFactors=FALSE)
+## summary(placette2010)
+placette2011 <- read.csv("./data/raw/DataFrance/2011/placettes_foret_2011.csv",sep=";"
+ , stringsAsFactors=FALSE)
+
+## names(placette2005)
+## ## uta -> uta1
+## ## sfo NA
+## ## plisi NA
+## names(placette2006)
+## ## plisi NA
+## names(placette2007)
+## ## dcespar1 dcespar2 tpespar1 tpespar2 iti pentexp NA
+## ## deleted acces
+## names(placette2008)
+## ## gest incid peupnr portance asperite
+## names(placette2009)
+## names(placette2010)
+## names(placette2011)
+
+placette2005$tplant <- as.character(placette2005$tplant)
+placette2005$tplant[placette2005$tplant==""] <-0
+placette2006$tplant[placette2006$tplant==""] <-0
+
+
+
+### for selection of plot use plisi=0 dc=0 tplant=0
+## incid indicateur d incident récent utilisé ou pas ?? avec ou sans perturbatrion naturelle
+
+placette_tot <-data.frame(idp=c(placette2005$idp,placette2006$idp,placette2007$idp,placette2008$idp,placette2009$idp,placette2010$idp,
+ placette2011$idp),
+ xl93=c(placette2005$xl93,placette2006$xl93,placette2007$xl93,placette2008$xl93,placette2009$xl93,placette2010$xl93,
+ placette2011$xl93),
+ yl93=c(placette2005$yl93,placette2006$yl93,placette2007$yl93,placette2008$yl93,placette2009$yl93,placette2010$yl93,
+ placette2011$yl93),
+ dep=c(placette2005$dep,placette2006$dep,placette2007$dep,placette2008$dep,placette2009$dep,placette2010$dep,
+ placette2011$dep),
+ csa=c(placette2005$csa,placette2006$csa,placette2007$csa,placette2008$csa,placette2009$csa,placette2010$csa,
+ placette2011$csa),
+ plisi=c(rep(NA,length(placette2005$tm2)),rep(NA,length(placette2006$idp)),placette2007$plisi,placette2008$plisi,
+ placette2009$plisi,placette2010$plisi,placette2011$plisi),
+ uta1=c(placette2005$uta,placette2006$uta1,placette2007$uta1,placette2008$uta1,placette2009$uta1,placette2010$uta1,
+ placette2011$uta1),
+ tm2=c(placette2005$tm2,placette2006$tm2,placette2007$tm2,placette2008$tm2,placette2009$tm2,placette2010$tm2,
+ placette2011$tm2),
+ sfo=c(rep(NA,length(placette2005$tm2)),placette2006$sfo,placette2007$sfo,placette2008$sfo,placette2009$sfo,
+ placette2010$sfo,placette2011$sfo),
+ incid=c(rep(NA,length(placette2005$idp)),rep(NA,length(placette2006$idp)),rep(NA,length(placette2007$idp)),
+ rep(NA,length(placette2008$idp)),placette2009$incid,placette2010$incid,placette2011$incid),
+ dc=c(placette2005$dc,placette2006$dc,placette2007$dc,placette2008$dc,placette2009$dc,placette2010$dc,placette2011$dc),
+ tplant=c(placette2005$tplant,placette2006$tplant,placette2007$tplant,placette2008$tplant,placette2009$tplant,
+ placette2010$tplant,placette2011$tplant),
+ esspre=c(placette2005$esspre,placette2006$esspre,placette2007$esspre,placette2008$esspre,placette2009$esspre,
+ placette2010$esspre,placette2011$esspre),
+ cac=c(placette2005$cac,placette2006$cac,placette2007$cac,placette2008$cac,placette2009$cac,placette2010$cac,
+ placette2011$cac),
+ ess_age_1=c(placette2005$ess_age_1,placette2006$ess_age_1,placette2007$ess_age_1,placette2008$ess_age_1,
+ placette2009$ess_age_1,placette2010$ess_age_1,placette2011$ess_age_1),
+ YEAR=c(rep(2005,length(placette2005$idp)),rep(2006,length(placette2006$idp)),rep(2007,length(placette2007$idp)),
+ rep(2008,length(placette2008$idp)),rep(2009,length(placette2009$idp)),rep(2010,length(placette2010$idp)),
+ rep(2011,length(placette2011$idp))))
+
+rm(placette2005,placette2006,placette2007,placette2008,placette2009,placette2010,placette2011)
+
+
+save(placette_tot,file="./data/process/placette_tot.Rdata")
+
+#####################################
+## LOAD elevation data
+##### LOAD ALTITUDE DATA
+## load("./data/process/placette_tot.Rdata")
+alti <- read.csv("./data/raw/DataFrance/altitude/SER_alti.csv",header=T,sep=";", stringsAsFactors=FALSE)
+alti2011 <- read.csv("./data/raw/DataFrance/altitude/SER_alti_2011.csv",header=T,sep=";",stringsAsFactors=FALSE)
+names(alti2011) <- names(alti)
+alti.tot <- rbind(alti,alti2011)
+
+## sum( placette_tot$idp %in% alti.tot$IDP)/length(placette_tot$idp)
+## table(placette_tot$YEAR[! placette_tot$idp %in% alti.tot$IDP])
+
+placette_tot.alti <- merge(placette_tot,alti.tot,by.x="idp",by.y="IDP")
+
+### write csv file for Piedallu
+write.csv(placette_tot.alti,"./data/process/placette_tot.alti.csv")
+
+### write new fiel for 2011
+write.csv(placette_tot.alti[placette_tot.alti$YEAR==2011,],"./data/process/placette_tot.alti.2011.csv")
+
+
+#########################################
+#########################################
+#########################################
+##### LOAD DATA ECOLOGIE
+### read DEAD TREE table downloaded from the web
+ecologie2005 <- read.csv("./data/raw/DataFrance/2005/ecologie_2005.csv",sep=";", stringsAsFactors=FALSE)
+## summary(ecologie2005)
+ecologie2006 <- read.csv("./data/raw/DataFrance/2006/ecologie_2006.csv",sep=";", stringsAsFactors=FALSE)
+## summary(ecologie2006)
+ecologie2007 <- read.csv("./data/raw/DataFrance/2007/ecologie_2007.csv",sep=";", stringsAsFactors=FALSE)
+## summary(ecologie2007)
+ecologie2008 <- read.csv("./data/raw/DataFrance/2008/ecologie_2008.csv",sep=";", stringsAsFactors=FALSE)
+## summary(ecologie2005)
+ecologie2009 <- read.csv("./data/raw/DataFrance/2009/ecologie_2009.csv",sep=";", stringsAsFactors=FALSE)
+## summary(ecologie2009)
+ecologie2010 <- read.csv("./data/raw/DataFrance/2010/ecologie_2010.csv",sep=";", stringsAsFactors=FALSE)
+## summary(ecologie2010)
+ecologie2011 <- read.csv("./data/raw/DataFrance/2011/ecologie_2011.csv",sep=";", stringsAsFactors=FALSE)
+
+
+ecologie_tot <- rbind(ecologie2005,ecologie2006,ecologie2007,ecologie2008,ecologie2009,ecologie2010,ecologie2011)
+head(ecologie_tot)
+rm(ecologie2005,ecologie2006,ecologie2007,ecologie2008,ecologie2009,ecologie2010,ecologie2011)
+
+save(ecologie_tot,file="./data/process/ecologie_tot.Rdata")
+
+
+