Skip to content
GitLab
Select Git revision
  • master default
  • BRIEFCOM
  • ctfs
  • test.dplyr
  • package
Find file BlameHistoryPermalink
  • fhui28's avatar
    test · 319aad84
    fhui28 authored
    319aad84
merge.data.FRANCE.R 8.88 KiB
1 
#############################################
2 
#############################################
3 
#############################################
4 
### MERGE FRENCH DATA
5 
source("./R/format.function.R")
6 
## source("./R/READ.DATA.NFI.FRANCE.R") need to be run if no already doen (long ...)
7 
## source("./R/CLIMATE.FRANCE.R")
8
9 
###########################################################################
10 
######################### READ DATA (test)
11
12 
### read individuals tree data
13 
load('./data/process/arbre.ALIVE.DEAD2.Rdata')
14
15 
### read climate
16 
ecologie.clim.data <- readRDS("./data/process/ecologie.clim.data.rds")
17 
#### merge
18 
dataIFN.FRANCE.t <- merge(arbre.ALIVE.DEAD2,
19 
                          subset(ecologie.clim.data,select=
20 
                                 c("idp","SER","sgdd","WB.y","WB.s","WS.y","WS.s","MAT","SAP")),by="idp")
21 
rm(arbre.ALIVE.DEAD2,ecologie.clim.data)
22
23 
#### load plot data and merge
24 
load("./data/process/placette_tot.Rdata")
25 
dataIFN.FRANCE <- merge(dataIFN.FRANCE.t,placette_tot[,names(placette_tot) != "YEAR"],by="idp")
26 
rm(placette_tot,dataIFN.FRANCE.t)
27
28 
### read IFN species names and clean
29 
species.clean <- fun.clean.species.tab(read.csv("./data/species.list/species.csv",
30 
                                                sep="\t",stringsAsFactors=FALSE))
31
32 
### read TRY data
33 
data.TRY.sd.update <-  readRDS("./data/process/data.TRY.sd.update.rds")
34 
data.frame.TRY <- data.frame(Latin_name=rownames(data.TRY.sd.update ),data.TRY.sd.update)
35 
rm(data.TRY.sd.update)
36 
### merge with code and species name
37 
merge.TRY <- merge(species.clean,data.frame.TRY,by="Latin_name")
38 
rm(species.clean,data.frame.TRY)
39
40
41
42
43
44 
########################################################################
45 
### Compute maximum height per species plus sd from observed height in NFI data to add variables to the traits data base
46 
res.quant.boot <- t(sapply(levels(factor(dataIFN.FRANCE[["espar"]])),
47 
                         FUN=f.quantile.boot,
48 
                         R=1000,
49 
                         x=log10(dataIFN.FRANCE[["htot"]]),
50 
                         fac=factor(dataIFN.FRANCE[["espar"]])))
51 
## create data base
52 
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])
53 
rm(res.quant.boot)
54 
## write.csv(data.max.height,file="./data/process/data.max.height.csv")
55
56 
##############
57 
## merge TRY with max height
58 
merge.TRY <- merge(merge.TRY,data.max.height,by="code")
59 
rm(data.max.height)
60 
# use mean sd of max tree height over all species
61 
merge.TRY$Max.height.sd.1 <- rep(mean(merge.TRY[["Max.height.sd"]],na.rm=TRUE),length=nrow(merge.TRY))
62
63 
### keep only variables needed in traits data
64 
names.traits.data <- c("code","Latin_name","Leaf.N.mean","Seed.mass.mean","SLA.mean","Wood.Density.mean",
65 
  "Leaf.Lifespan.mean","Max.height.mean","Leaf.N.sd.1","Seed.mass.sd.1","SLA.sd.1", "Wood.Density.sd.1",
66 
  "Leaf.Lifespan.sd.1","Max.height.sd.1")
67
68 
data.traits <- merge.TRY[,names.traits.data]
69 
names(data.traits) <- c("sp","Latin_name","Leaf.N.mean","Seed.mass.mean","SLA.mean","Wood.Density.mean",
70 
  "Leaf.Lifespan.mean","Max.height.mean","Leaf.N.sd","Seed.mass.sd","SLA.sd", "Wood.Density.sd",
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
"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), sort=FALSE)
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203
########################################################################################### ### 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") #### 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")