diff --git a/merge.data.FRANCE.R b/merge.data.FRANCE.R
index 277893a41e237382a76e042af5627d3bdd6a43b0..3a935444cf156cc2afb2d219b47116139b66b801 100644
--- a/merge.data.FRANCE.R
+++ b/merge.data.FRANCE.R
@@ -1,11 +1,12 @@
############################################# MERGE FRENCH DATA
+rm(list = ls());
source("./R/format.function.R")
+library(reshape)
-################################ READ DATA (test)
-dataIFN.FRANCE <- read.csv("./data/raw/DataFrance/dataIFN.FRANCE.csv", stringsAsFactors = FALSE)
+################################ READ DATA
+data.france <- read.csv("./data/raw/DataFrance/dataIFN.FRANCE.csv", stringsAsFactors = FALSE)
### read IFN species names and clean
-species.clean <- fun.clean.species.tab(read.csv("./data/raw/DataFrance/species.csv",
- sep = "\t", stringsAsFactors = FALSE))
+species.clean <- fun.clean.species.tab(read.csv("./data/raw/DataFrance/species.csv", stringsAsFactors = FALSE))
### read TRY data
data.TRY.sd.update <- readRDS("./data/process/data.TRY.sd.update.rds")
@@ -16,23 +17,22 @@ merge.TRY <- merge(species.clean, data.frame.TRY, by = "Latin_name")
rm(species.clean, data.frame.TRY)
+###################################### 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"]])))
-
-
-######################################################################## 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 with max height
+## merge TRY with max height
merge.TRY <- merge(merge.TRY, data.max.height, by = "code")
rm(data.max.height)
-# use mean sd of max tree height over all species
+## 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))
@@ -42,103 +42,88 @@ names.traits.data <- c("code", "Latin_name", "Leaf.N.mean", "Seed.mass.mean", "S
"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
+names(data.traits) <- c("espar", "latin_name", "leafN.mean", "seedmass.mean", "SLA.mean",
+ "wooddensity.mean", "leaflifespan.mean", "maxheight.mean", "leafN.sd", "seedmass.sd",
+ "SLA.sd", "wooddensity.sd", "leaflifespan.sd", "maxheight.sd")
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
+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$plot.id <- (data.france[["idp"]]); data.france$idp <- NULL ## plot code
+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
-
-
-################################################################ 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
+######################## 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")]
+data.sp <- data.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"]
+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')
-############################ merge greco to have no ecoregion with low number of observation
+###################### 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
-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 <- 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(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)),
+## 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
-## 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
+## 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(dataIFN.FRANCE, select = c(vec.basic.var, vec.abio.var.names))
-
-
-
+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
-
-
-################################################################################################################## 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 *
+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
@@ -146,7 +131,7 @@ 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")
+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
@@ -154,13 +139,12 @@ 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 <- 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")
diff --git a/merge.data.SWISS.R b/merge.data.SWISS.R
index eb7ba34e585922074b14ae36c480002ad88bbb56..1e975b98d7c6bd9b91f4be6c527170565b598b48 100644
--- a/merge.data.SWISS.R
+++ b/merge.data.SWISS.R
@@ -94,7 +94,6 @@ data.swiss <- merge(data.swiss, data.frame(siteid = data.climate$CLNR, swb = dat
MAT = data.climate$tave_68, MAP = data.climate$MAP, stringsAsFactors =FALSE), sort = F, all.x = T)
rm(data.climate)
-##############################################
##############################################
#################### GENERATE ONE OBJECT PER ECOREGION
# vector of ecoregion name