France.R

#!/usr/bin/env Rscript
############################################# MERGE FRENCH DATA
rm(list = ls()); 
source("R/format.data/format.fun.R")
library(reshape)
dir.create("output/formatted/France", recursive=TRUE,showWarnings=FALSE)


################################ READ DATA
data.france <- read.csv("data/raw/France/dataIFN.FRANCE.csv", stringsAsFactors = FALSE)
### read IFN species names and clean
species.clean <- fun.clean.species.tab(read.csv("data/raw/France/species.csv", stringsAsFactors = FALSE))

######## MASSAGE TRAIT DATA Compute maximum height per species plus sd from observed
######## height to add variables to the traits data base Because we have two heights,
######## then take the max of the two heights and then bootstrap
res.quant.boot <- t(sapply(levels(factor(data.france[["espar"]])), FUN = f.quantile.boot, 
    R = 1000, x = log10(data.france[["htot"]]), fac = factor(data.france[["espar"]])))
## create data base
data.max.height <- data.frame(sp = 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='output/formatted/France/max.height.csv')
rm(data.max.height)

########################################## FORMAT INDIVIDUAL TREE DATA change unit and names of variables to be the same
########################################## in all data for the tree
data.france$G <- data.france[["ir5"]]/5 * 2  ## diameter growth in  mm per year
data.france$year <- rep(5, nrow(data.france))  ## number of year between measurement
data.france$D <- data.france[["c13"]]/pi  ## diameter in cm
data.france$sp <- as.character(data.france[["espar"]]); data.france$espar <- NULL ## species code
data.france$sp.name <- species.clean[match(data.france$sp,species.clean$sp),"Latin_name"]
data.france$plot <- (data.france[["idp"]]); data.france$idp <- NULL ## plot code
data.france$cluster <- rep(NA,nrow(data.france))
data.france$tree.id <- paste(data.france[["plot.id"]], data.france[["a"]], sep = "_")  ## tree unique id
data.france$weights <- data.france[["w"]]/10000
data.france$obs.id <- 1:nrow(data.france) ## There is only obs per tree.id, so this is superfluous
data.france$census <- rep(1,nrow(data.france)) ## only one census
######################## change coordinates system of x y to be in lat long WGS84
library(sp)
library(dismo)
library(rgdal)
data.sp <- data.france[, c("tree.id", "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
data.france$Lon <- coordinates(data.sp2)[, "xl93"]
data.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')
###################### ECOREGION - 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
GRECO.temp <- substr(data.france[["SER"]], 1, 1)  ## get GRECO from SER (smaller division by keeping only the first letter)
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(data.france[['xl93']],data.france[['yl93']],col=unclass(factor(GRECO.temp)))
data.france$ecocode <- GRECO.temp  ## a single code for each ecoregion
###################### PERCENT DEAD variable percent dead/cannot do with since dead variable is
###################### missing compute numer of dead per plot to remove plot with disturbance
perc.dead <- tapply(data.france[["dead"]], INDEX = data.france[["plot"]], FUN = function.perc.dead2)
data.france <- merge(data.france, data.frame(plot = as.numeric(names(perc.dead)), 
    perc.dead = perc.dead),by="plot", sort = FALSE)
########################################################################################### PLOT SELECTION FOR THE ANALYSIS