Skip to content
GitLab
Select Git revision
  • master default
  • BRIEFCOM
  • ctfs
  • test.dplyr
  • package
Find file BlameHistoryPermalink
merge.data.FRANCE.R 7.67 KiB
1 
############################################# MERGE FRENCH DATA
2 
rm(list = ls());
3 
source("./R/format.function.R")
4 
library(reshape)
5
6 
################################ READ DATA
7 
data.france <- read.csv("./data/raw/DataFrance/dataIFN.FRANCE.csv", stringsAsFactors = FALSE)
8 
### read IFN species names and clean
9 
species.clean <- fun.clean.species.tab(read.csv("./data/raw/DataFrance/species.csv", stringsAsFactors = FALSE))
10
11 
### read TRY data
12 
data.TRY.sd.update <- readRDS("./data/process/data.TRY.sd.update.rds")
13 
data.frame.TRY <- data.frame(Latin_name = rownames(data.TRY.sd.update), data.TRY.sd.update)
14 
rm(data.TRY.sd.update)
15 
### merge with code and species name
16 
merge.TRY <- merge(species.clean, data.frame.TRY, by = "Latin_name")
17 
rm(species.clean, data.frame.TRY)
18
19
20 
###################################### MASSAGE TRAIT DATA Compute maximum height per species plus sd from observed
21 
###################################### height to add variables to the traits data base Because we have two heights,
22 
###################################### then take the max of the two heights and then bootstrap
23 
res.quant.boot <- t(sapply(levels(factor(data.france[["espar"]])), FUN = f.quantile.boot,
24 
    R = 1000, x = log10(data.france[["htot"]]), fac = factor(data.france[["espar"]])))
25
26 
## create data base
27 
data.max.height <- data.frame(code = rownames(res.quant.boot), Max.height.mean = res.quant.boot[,
28 
    1], Max.height.sd = res.quant.boot[, 2], Max.height.nobs = res.quant.boot[, 3])
29 
rm(res.quant.boot)
30 
## write.csv(data.max.height,file='./data/process/data.max.height.csv')
31
32 
## merge TRY with max height
33 
merge.TRY <- merge(merge.TRY, data.max.height, by = "code")
34 
rm(data.max.height)
35 
## use mean sd of max tree height over all species
36 
merge.TRY$Max.height.sd.1 <- rep(mean(merge.TRY[["Max.height.sd"]], na.rm = TRUE),
37 
    length = nrow(merge.TRY))
38
39 
### keep only variables needed in traits data
40 
names.traits.data <- c("code", "Latin_name", "Leaf.N.mean", "Seed.mass.mean", "SLA.mean",
41 
    "Wood.Density.mean", "Leaf.Lifespan.mean", "Max.height.mean", "Leaf.N.sd.1",
42 
    "Seed.mass.sd.1", "SLA.sd.1", "Wood.Density.sd.1", "Leaf.Lifespan.sd.1", "Max.height.sd.1")
43
44 
data.traits <- merge.TRY[, names.traits.data]
45 
names(data.traits) <- c("espar", "latin_name", "leafN.mean", "seedmass.mean", "SLA.mean",
46 
    "wooddensity.mean", "leaflifespan.mean", "maxheight.mean", "leafN.sd", "seedmass.sd",
47 
    "SLA.sd", "wooddensity.sd", "leaflifespan.sd", "maxheight.sd")
48 
rm(merge.TRY, names.traits.data)
49
50
51 
########################################## FORMAT INDIVIDUAL TREE DATA change unit and names of variables to be the same
52 
########################################## in all data for the tree
53
54 
data.france$G <- data.france[["ir5"]]/5 * 2  ## diameter growth in  mm per year
55 
data.france$year <- rep(5, nrow(data.france))  ## number of year between measurement
56 
data.france$D <- data.france[["c13"]]/pi  ## diameter in cm
57 
data.france$sp <- as.character(data.france[["espar"]]); data.france$espar <- NULL ## species code
58 
data.france$plot.id <- (data.france[["idp"]]); data.france$idp <- NULL ## plot code
59 
data.france$tree.id <- paste(data.france[["plot.id"]], data.france[["a"]], sep = "_")  ## tree unique id
60 
data.france$weights <- data.france[["w"]]/10000
61 
data.france$obs.id <- 1:nrow(data.france) ## There is only obs per tree.id, so this is superfluous
62
63
64 
######################## change coordinates system of x y to be in lat long WGS84
65 
library(sp)
66 
library(dismo)
67 
library(rgdal)
68 
data.sp <- data.france[, c("idp", "xl93", "yl93")]
69 
coordinates(data.sp) <- c("xl93", "yl93")  ## EPSG CODE 2154
70 
proj4string(data.sp) <- CRS("+init=epsg:2154")  # define projection system of our data ## EPSG CODE 2154
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
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.id"]], FUN = function.perc.dead2) data.france <- merge(data.france, data.frame(plot.id = as.numeric(names(perc.dead)), perc.dead = perc.dead), sort = FALSE) ########################################################################################### PLOT SELECTION FOR THE ANALYSIS ## data.france <- subset(data.france,subset= data.france[['YEAR']] != 2005) ## year 2005 bad data according to IFN data.france <- subset(data.france, subset = data.france[["plisi"]] == 0) ## no plot on forest edge data.france <- subset(data.france, subset = data.france[["dc"]] == 0) ## no harvesting data.france <- subset(data.france, subset = data.france[["tplant"]] == 0) ## no plantation data.france <- subset(data.france, subset = !is.na(data.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(data.france, select = c(vec.basic.var, vec.abio.var.names)) ############################################## #################### GENERATE ONE OBJECT PER ECOREGION - 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(data.france[["tree.id"]]), diam = as.vector(data.france[["D"]]), sp = as.vector(data.france[["sp"]]), id.plot = as.vector(data.france[["idp"]]), weights = as.vector(data.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(data.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
141142143144145146147148149150151
names(data.BA.SP) <- c("tree.id", names(data.BA.SP)[-1]) data.BA.sp <- merge(data.frame(tree.id = data.france[["tree.id"]], ecocode = data.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")