Skip to content
GitLab
Select Git revision
  • master default
  • BRIEFCOM
  • ctfs
  • test.dplyr
  • package
Find file BlameHistoryPermalink
merge.data.SPAIN.R 8.03 KiB
1 
### MERGE spain DATA Edited by FH
2 
rm(list = ls())
3 
source("./R/format.function.R")
4 
library(reshape)
5
6 
######################### READ DATA read individuals tree data data.spain <-
7 
######################### read.table('./data/raw/DataSpain/Tree_data_SFI.txt',header=TRUE,stringsAsFactors=FALSE,sep
8 
######################### = '\t')
9 
data.spain <- read.table("./data/raw/DataSpain/Tree_data_SFI_aug13_alldata.txt",
10 
    header = TRUE, stringsAsFactors = FALSE, sep = "\t")
11
12 
###################################### MASSAGE TRAIT DATA Compute maximum height per species plus sd from observed
13 
###################################### height to add variables to the traits data base Because we have two heights,
14 
###################################### then take the max of the two heights and then bootstrap
15 
res.quant.boot <- t(sapply(levels(factor(data.spain[["SP_code"]])), FUN = f.quantile.boot,
16 
    R = 1000, x = log10(apply(data.spain[, c("ht1", "ht2")], 1, max, na.rm = T)),
17 
    fac = factor(data.spain[["SP_code"]])))
18
19 
## create data base
20 
data.max.height <- data.frame(code = rownames(res.quant.boot), Max.height.mean = res.quant.boot[,
21 
    1], Max.height.sd = res.quant.boot[, 2], Max.height.nobs = res.quant.boot[, 3])
22 
rm(res.quant.boot)
23 
write.csv(data.max.height, file = "./data/process/data.max.height.spain.csv")  # I was planning to save processed data in that folder
24 
# ## merge TRY with max height merge.TRY <-
25 
# merge(merge.TRY,data.max.height,by='code') rm(data.max.height) ## use mean sd
26 
# of max tree height over all species merge.TRY$Max.height.sd.1 <-
27 
# rep(mean(merge.TRY[['Max.height.sd']],na.rm=TRUE),length=nrow(merge.TRY)) ###
28 
# keep only variables needed in traits data names.traits.data <-
29 
# c('code','Latin_name','Leaf.N.mean','Seed.mass.mean','SLA.mean','Wood.Density.mean',
30 
# 'Leaf.Lifespan.mean','Max.height.mean','Leaf.N.sd.1','Seed.mass.sd.1','SLA.sd.1',
31 
# 'Wood.Density.sd.1', 'Leaf.Lifespan.sd.1','Max.height.sd.1') data.traits <-
32 
# merge.TRY[,names.traits.data] names(data.traits) <-
33 
# c('sp','Latin_name','Leaf.N.mean','Seed.mass.mean','SLA.mean','Wood.Density.mean',
34 
# 'Leaf.Lifespan.mean','Max.height.mean','Leaf.N.sd','Seed.mass.sd','SLA.sd',
35 
# 'Wood.Density.sd', 'Leaf.Lifespan.sd','Max.height.sd') ## rename to have
36 
# standard variables name rm(merge.TRY,names.traits.data)
37
38 
################################################################ FORMAT INDIVIDUAL TREE DATA
39
40 
## change unit and names of variables to be the same in all data for the tree
41 
data.spain$G <- data.spain[["adbh"]]  ## diameter growth in mm per year
42 
data.spain$year <- data.spain[["years"]]  ## number of year between measurement - MISSING
43 
data.spain$D <- data.spain[["dbh1"]]/10  ## diameter in mm convert to cm
44 
data.spain$dead <- as.numeric(data.spain[["Life_status"]] == "dead")  ## dummy variable for dead tree 0 alive 1 dead - MIGHT WANT TO CHANGE THIS TO BE BASED ON MORTALITY_CUT
45 
data.spain$sp <- as.character(data.spain[["SP_code"]])  ## species code
46 
data.spain$plot <- (data.spain[["Plot_ID_SFI"]])  ## plot code
47 
data.spain$htot <- data.spain[["ht1"]]  ## height of tree in m
48 
data.spain$tree.id <- data.spain$Tree_ID_SFI  ## tree unique id
49
50 
#### change coordinates system of x y to be in lat long WGS84/don't know how to do
51 
#### this
52 
library(sp)
53 
library(dismo)
54 
library(rgdal)
55
56 
data.sp <- data.spain[, c("Tree_ID_SFI", "CX", "CY")]
57 
coordinates(data.sp) <- c("CX", "CY")  # define x y
58 
proj4string(data.sp) <- CRS("+init=epsg:23030")  # define projection system of our data ## EPSG CODE 23030 ED50 / UTM zone 30N
59 
summary(data.sp)
60
61 
detach(package:rgdal)
62 
data.sp2 <- spTransform(data.sp, CRS("+init=epsg:4326"))  ## change projection in WGS84 lat lon
63 
data.spain$Lon <- coordinates(data.sp2)[, "CX"]
64 
data.spain$Lat <- coordinates(data.sp2)[, "CY"]
65 
## ## plot on world map library(rworldmap) newmap <- getMap(resolution = 'coarse')
66 
## # different resolutions available plot(newmap)
67 
## points(data.sp2,cex=0.2,col='red')
68 
rm(data.sp, data.sp2)
69
70 
###################### ECOREGION merge greco to have no ecoregion with low number of observation
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
greco <- read.csv(file = "./data/raw/DataSpain/R_Ecoregion.csv", header = T) greco <- greco[, c("Plot_ID_SFI", "BIOME", "eco_code")] greco2 <- greco[!duplicated(greco$Plot), ] rm(greco) data.spain <- merge(data.spain, greco2, by = "Plot_ID_SFI") rm(greco2) table(data.spain$eco_code) ## There's an eco-region with no code, and one with < 1000 sites The former we ## could drop as they were on the border of Spain library(RColorBrewer) mycols <- brewer.pal(10, "Set3") ecoreg <- unclass(data.spain$eco_code) plot(data.spain[["CX"]][order(ecoreg)], data.spain[["CY"]][order(ecoreg)], pty = ".", cex = 0.2, col = rep(mycols, as.vector(table(ecoreg)))) legend("topleft", col = mycols, legend = levels(data.spain$eco_code), pch = rep(19, length(levels(ecoreg))), cex = 2) points(data.spain[["CX"]][ecoreg == 9], data.spain[["CY"]][ecoreg == 9], pty = ".", cex = 0.5, col = "black") ## Highlight the 'rare' ecoregions points(data.spain[["CX"]][ecoreg == 1], data.spain[["CY"]][ecoreg == 1], pty = ".", cex = 0.5, col = "black") ## Highlight the 'rare' ecoregions PA1219 looks to be similar to PA1209, merge ## them together data.spain$eco_codemerged <- combine_factor(data.spain$eco_code, c(1:8, 6, 9)) data.spain <- data.spain[-which(data.spain$eco_codemerged == ""), ] ###################### 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.spain[["dead"]], INDEX = data.spain[["plot"]], FUN = function.perc.dead) table(data.spain$dead) ## VARIABLE TO SELECT PLOT WITH NOT BIG DISTURBANCE KEEP OFTHER VARIABLES IF ## AVAILABLE (disturbance record) data.spain <- merge(data.spain, data.frame(plot = as.numeric(names(perc.dead)), perc.dead = perc.dead), sort = FALSE, by = "plot") ########################################################### PLOT SELECTION FOR THE ANALYSIS Remove data with mortality == 1 or 2 table(data.spain$Mortality_Cut) data.spain <- subset(data.spain, subset = (data.spain[["Mortality_Cut"]] == 0 | data.spain[["Mortality_Cut"]] == "")) colnames(data.spain)[colnames(data.spain) %in% c("mat", "pp", "PET")] <- c("MAT", "PP", "PET") colnames(data.spain)[names(data.spain) == "eco_codemerged"] <- c("ecocode") vec.abio.var.names <- c("MAT", "PP", "PET") vec.basic.var <- c("tree.id", "sp", "sp.name", "plot", "ecocode", "D", "G", "dead", "year", "htot", "Lon", "Lat", "perc.dead") data.tree <- subset(data.spain, select = c(vec.basic.var, vec.abio.var.names)) save(data.spain, file = "./data/process/datspain.RData") ############################################## 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.spain[["tree.id"]]), diam = as.vector(data.spain[["D"]]), sp = as.vector(data.spain[["sp"]]), id.plot = as.vector(data.spain[["plot"]]), weights = as.vector(1/(pi * (data.spain[["R1"]])^2)), weight.full.plot = 1/(pi * (25)^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.spain[["sp"]]))) > 0) stop("competition index sp name not the same as in data.tree")
141142143144145146147148149150151152153154
#### 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.spain[["tree.id"]], ecocode = dataIFN.spain[["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.spain <- list(data.tree = data.tree, data.BA.SP = data.BA.sp, data.traits = data.traits) save(list.spain, file = "./data/process/list.spain.Rdata")