#############################################
#############################################
#############################################
### MERGE FRENCH DATA
source("./R/format.function.R")

###########################################################################
######################### READ DATA

### read individuals tree data
load('./data/process/arbre.ALIVE.DEAD2.Rdata')

### read climate
ecologie.clim.data <- readRDS("./data/process/ecologie.clim.data.rds")
#### merge
dataIFN.FRANCE.t <- merge(arbre.ALIVE.DEAD2,
                          subset(ecologie.clim.data,select=
                                 c("idp","SER","sgdd","WB.y","WB.s","WS.y","WS.s","MAT","SAP")),by="idp")
rm(arbre.ALIVE.DEAD2,ecologie.clim.data)

#### load plot data and merge
load("./data/process/placette_tot.Rdata")
dataIFN.FRANCE <- merge(dataIFN.FRANCE.t,placette_tot[,names(placette_tot) != "YEAR"],by="idp")
rm(placette_tot,dataIFN.FRANCE.t)

### read IFN species names and clean
species.clean <- fun.clean.species.tab(read.csv("./data/species.list/species.csv",
                                                sep="\t",stringsAsFactors=FALSE)) 

### read TRY data
data.TRY.sd.update <-  readRDS("./data/process/data.TRY.sd.update.rds")
data.frame.TRY <- data.frame(Latin_name=rownames(data.TRY.sd.update ),data.TRY.sd.update)
rm(data.TRY.sd.update)
### merge with code and species name
merge.TRY <- merge(species.clean,data.frame.TRY,by="Latin_name")
rm(species.clean,data.frame.TRY)





########################################################################
### Compute maximum height per species plus sd from observed height in NFI data to add variables to the traits data base
res.quant.boot <- t(sapply(levels(factor(dataIFN.FRANCE[["espar"]])),
                         FUN=f.quantile.boot,
                         R=1000,
                         x=log10(dataIFN.FRANCE[["htot"]]),
                         fac=factor(dataIFN.FRANCE[["espar"]])))
## create data base
data.max.height <- data.frame(code=rownames(res.quant.boot),Max.height.mean=res.quant.boot[,1],Max.height.sd=res.quant.boot[,2],Max.height.nobs=res.quant.boot[,3])
rm(res.quant.boot)
## write.csv(data.max.height,file="./data/process/data.max.height.csv")

##############
## merge TRY with max height
merge.TRY <- merge(merge.TRY,data.max.height,by="code")
rm(data.max.height)
# use mean sd of max tree height over all species
merge.TRY$Max.height.sd.1 <- rep(mean(merge.TRY[["Max.height.sd"]],na.rm=TRUE),length=nrow(merge.TRY))

### keep only variables needed in traits data
names.traits.data <- c("code","Latin_name","Leaf.N.mean","Seed.mass.mean","SLA.mean","Wood.Density.mean",
  "Leaf.Lifespan.mean","Max.height.mean","Leaf.N.sd.1","Seed.mass.sd.1","SLA.sd.1", "Wood.Density.sd.1",
  "Leaf.Lifespan.sd.1","Max.height.sd.1")

data.traits <- merge.TRY[,names.traits.data]
names(data.traits) <- c("sp","Latin_name","Leaf.N.mean","Seed.mass.mean","SLA.mean","Wood.Density.mean",
  "Leaf.Lifespan.mean","Max.height.mean","Leaf.N.sd","Seed.mass.sd","SLA.sd", "Wood.Density.sd",
  "Leaf.Lifespan.sd","Max.height.sd") ## rename to have standard variables name
rm(merge.TRY,names.traits.data)







################################################################
#############  FORMAT INDIVIDUAL TREE DATA
#############

## change unit and names of variables to be the same in all data for the tree 
dataIFN.FRANCE$G <- dataIFN.FRANCE[["ir5"]]/5*2 ##diameter growth in  mm per year
dataIFN.FRANCE$year <- rep(5,length(dataIFN.FRANCE[["ir5"]])) ## number of year between measurement
dataIFN.FRANCE$D <- dataIFN.FRANCE[["c13"]]/pi ## diameter in cm
dataIFN.FRANCE$dead <- dataIFN.FRANCE[["dead"]] ## dummy variable for dead tree 0 alive 1 dead
dataIFN.FRANCE$sp <- as.character(dataIFN.FRANCE[["espar"]]) ## species code
dataIFN.FRANCE$plot <- (dataIFN.FRANCE[["idp"]]) ## plot code
dataIFN.FRANCE$htot <- (dataIFN.FRANCE[["htot"]]) ## height of tree in m
dataIFN.FRANCE$tree.id <- paste(dataIFN.FRANCE[["idp"]],dataIFN.FRANCE[["a"]],sep=".") ## tree unique id


#### change coordinates system of x y to be in lat long WGS84
library(sp)
library(dismo)
library(rgdal)
data.sp <-  dataIFN.FRANCE[,c("idp","xl93","yl93")]
coordinates(data.sp) <- c("xl93", "yl93") ## EPSG CODE 2154
proj4string(data.sp) <- CRS("+init=epsg:2154")  # define projection system of our data ## EPSG CODE 2154
summary(data.sp)

data.sp2 <- spTransform(data.sp,CRS("+init=epsg:4326")) ## change projection in WGS84 lat lon
dataIFN.FRANCE$Lon <- coordinates(data.sp2)[,"xl93"]
dataIFN.FRANCE$Lat <- coordinates(data.sp2)[,"yl93"]
rm(data.sp,data.sp2)
## ## plot on world map
## library(rworldmap)
## newmap <- getMap(resolution = "coarse")  # different resolutions available
## plot(newmap)
## points(data.sp2,cex=0.2,col="red")

############################
## merge greco to have no ecoregion with low number of observation

## merge A and B Grand Ouest cristallin and oceanique and Center semi-oceanique 
## merge G D E Vosges Jura massif cemtral (low mountain)
## merge H and I Alpes and Pyrenees
## Merge J and K Corse and Mediteraneen
dataIFN.FRANCE$GRECO <- substr( dataIFN.FRANCE[["SER"]],1,1) ## get GRECO from SER (smaller division by keeping only the first letter

GRECO.temp <- dataIFN.FRANCE[["GRECO"]]
GRECO.temp <- sub("[AB]","AB",GRECO.temp)
GRECO.temp <- sub("[GDE]","GDE",GRECO.temp)
GRECO.temp <- sub("[HI]","HI",GRECO.temp)
GRECO.temp <- sub("[JK]","JK",GRECO.temp)
## plot(dataIFN.FRANCE[["xl93"]],dataIFN.FRANCE[["yl93"]],col=unclass(factor(GRECO.temp)))
## add NEW GRECO variable to data base
dataIFN.FRANCE$ecocode <- GRECO.temp ## a single code for each ecoregion

## variable percent dead
#######################
###compute numer of dead per plot to remove plot with disturbance
perc.dead <- tapply(dataIFN.FRANCE[["dead"]],INDEX=dataIFN.FRANCE[["idp"]],FUN=function.perc.dead)
## VARIABLE TO SELECT PLOT WITH NOT BIG DISTURBANCE KEEP OFTHER VARIABLES IF AVAILABLE (disturbance record)
dataIFN.FRANCE <- merge(dataIFN.FRANCE,data.frame(idp=as.numeric(names(perc.dead)),perc.dead=perc.dead))


###########################################################################################
### PLOT SELECTION FOR THE ANALYSIS

## dataIFN.FRANCE <- subset(dataIFN.FRANCE,subset= dataIFN.FRANCE[["YEAR"]] != 2005)## year 2005 bad data according to IFN
dataIFN.FRANCE <- subset(dataIFN.FRANCE,subset= dataIFN.FRANCE[["plisi"]] == 0) # no plot on forest edge
dataIFN.FRANCE <- subset(dataIFN.FRANCE,subset= dataIFN.FRANCE[["dc"]] == 0) # no harvesting
dataIFN.FRANCE <- subset(dataIFN.FRANCE,subset= dataIFN.FRANCE[["tplant"]] == 0) # no plantation
dataIFN.FRANCE <- subset(dataIFN.FRANCE,subset= !is.na(dataIFN.FRANCE[["SER"]])) # missing SER

####################################
###################################
## SELECT GOOD COLUMNS

##names of variables for abiotic conditions
vec.abio.var.names <-  c("MAT","SAP","sgdd","WB.s","WB.y","WS.s","WS.y")
## other var
vec.basic.var <-  c("tree.id","sp","plot","ecocode","D","G","dead","year","htot","Lon","Lat","perc.dead")

data.tree <- subset(dataIFN.FRANCE,select=c(vec.basic.var,vec.abio.var.names))







##################################################################################################################
##################################################################################################################
##################################################################################################################
##############################################
## COMPUTE MATRIX OF COMPETITION INDEX WITH SUM OF BA PER SPECIES IN EACH PLOT in m^2/ha without the target species

data.BA.SP <-  BA.SP.FUN(id.tree=as.vector(dataIFN.FRANCE[["tree.id"]]),
diam=as.vector(dataIFN.FRANCE[["D"]]),
sp=as.vector(dataIFN.FRANCE[["sp"]]),
id.plot=as.vector(dataIFN.FRANCE[["idp"]]),
weights=as.vector(dataIFN.FRANCE[["w"]])/10000,
weight.full.plot=1/(pi*(c(15))^2))

## change NA and <0 data for 0
data.BA.SP[which(is.na(data.BA.SP),arr.ind=TRUE)] <- 0
data.BA.SP[,-1][which(data.BA.SP[,-1]<0,arr.ind=TRUE)] <- 0

### CHECK IF sp and sp name for column are the same
if(sum(! (names(data.BA.SP)[-1] %in% unique(dataIFN.FRANCE[["sp"]]))) >0) stop("competition index sp name not the same as in data.tree")

## if greater than 0 PROBLEM

#### compute BA tot for all competitors
BATOT.COMPET <- apply(data.BA.SP[,-1],MARGIN=1,FUN=sum,na.rm=TRUE)
data.BA.SP$BATOT.COMPET <- BATOT.COMPET
### create data frame
names(data.BA.SP) <- c("tree.id",names(data.BA.SP)[-1])
data.BA.sp <- merge(data.frame(tree.id=dataIFN.FRANCE[["tree.id"]],ecocode=dataIFN.FRANCE[["ecocode"]]),data.BA.SP,by="tree.id",sort=FALSE)
## test
if( sum(! data.BA.sp[["tree.id"]] == data.tree[["tree.id"]]) >0) stop("competition index not in the same order than data.tree")



## save everything as a list
list.FRANCE <- list(data.tree=data.tree,data.BA.SP=data.BA.sp,data.traits=data.traits)
save(list.FRANCE,file="./data/process/list.FRANCE.Rdata")