Newer
Older
#############################################
#############################################
#############################################
### MERGE FRENCH DATA
source("./R/format.function.R")
###########################################################################
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 <- 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"]
## ## 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
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
## 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))
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")
#### 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)
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")