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