Skip to content
GitLab
Select Git revision
  • master default
  • BRIEFCOM
  • ctfs
  • test.dplyr
  • package
Find file BlameHistoryPermalink
merge.data.PARACOU.R 10.84 KiB
1 
### MERGE paracou DATA
2 
rm(list = ls())
3 
source("./R/format.function.R")
4 
library(reshape)
5
6 
############################ read individuals tree data
7 
data.paracou <- read.table("./data/raw/DataParacou/20130717_paracou_1984_2012.csv",
8 
                           header=TRUE,stringsAsFactors=FALSE,sep = ";", na.strings = "NULL")
9 
#barplot(apply(!is.na(data.paracou[,paste("circ_",1984:2012,sep="")]),MARGIN=2,FUN=sum),las=3)
10
11 
# select good columns
12 
data.paracou <- data.paracou[,c("foret","parcelle","carre","arbre","vernaculaire","idtaxon",
13 
                                "x","y","circ_2001","code_2001","circ_2005","code_2005",
14 
                                "circ_2009","code_2009","campagne_mort","type_mort")]
15 
colnames(data.paracou) <- c("forest","plot","subplot","tree","vernacular","taxonid","x","y",
16 
                            "circum2001","code2001","circum2005","code2005","circum2009",
17 
                            "code2009","yeardied","typedeath")
18
19 
### change numeric separator
20 
numeric.col.name <-   c("x","y","circum2001","code2001","circum2005","code2005","circum2009","code2009")
21 
for(k in numeric.col.name){
22 
	data.paracou[,k] <- gsub(",",".",data.paracou[,k]); data.paracou[,k] <- as.numeric(data.paracou[,k])
23 
    } ## Replace all , in decimals with .
24
25 
data.paracou$tree.id <- apply(data.paracou[,c("plot","subplot","tree")],1,paste,collapse="_");
26 
data.paracou$sp <- data.paracou[["taxonid"]]
27 
data.paracou <- data.paracou[,c(ncol(data.paracou),1:(ncol(data.paracou)-1))]
28
29 
## ## plot each plot
30 
## pdf("./figs/plots.paracou.pdf")
31 
## lapply(unique(data.paracou[["plot"]]),FUN=fun.circles.plot,data.paracou[['x']],data.paracou[['y']],data.paracou[["plot"]],data.paracou[["circum2009"]],inches=0.2)
32 
## dev.off()
33
34 
############################# SELECT OBSERVATION WITHOUT PROBLEMS
35 
## REMOVE ALL TREES WITH X OR Y >250 m
36 
data.paracou <- subset(data.paracou,subset=(!is.na(data.paracou[["x"]])) & data.paracou[["x"]]<251 &  data.paracou[["y"]]<251)
37 
#### REMOVE PLOTs 16 17 18 ACCORDING TO  GHSILAIN
38 
data.paracou <- subset(data.paracou,subset=! data.paracou[["plot"]] %in% 16:18)
39 
## keep only tree alive in 2001
40 
data.paracou <- subset(data.paracou,subset=!(as.numeric(data.paracou[["yeardied"]])<=2001 & !is.na(data.paracou[["yeardied"]])))
41
42
43 
######################################## MASSAGE TRAIT DATA
44
45 
### read species names
46 
species.clean <- read.csv("./data/raw/DataParacou/20130717_paracou_taxonomie.csv",stringsAsFactors=FALSE, header = T, sep = ";")
47 
species.clean$sp <- species.clean[["idTaxon"]]
48 
species.clean$Latin_name <-  paste(species.clean[["Genre"]],species.clean[["Espece"]],sep=" ")
49 
## keep only one row pers idTaxon
50 
species.clean <- subset(species.clean,subset=!duplicated(species.clean[["sp"]]),select=c("sp","Latin_name","Genre","Espece","Famille","idCIRAD"))
51
52 
## select only species present in data base
53 
species.clean <-  subset(species.clean,subset=species.clean[["sp"]] %in% data.paracou[["sp"]])
54 
## percentage of species with no taxonomic identification
55 
length(grep("Indet",species.clean[["Latin_name"]]))/nrow(species.clean) ## 25%
56
57 
dataWD <- read.csv("./data/raw/DataParacou/WD-Species-Paracou-Ervan_GV.csv",stringsAsFactors=FALSE, header = T,sep=" ")
58 
#dataWD <- merge(dataWD, species.clean, by = "idCIRAD", sort = F)
59 
length(unique(species.clean$idCIRAD)) !=  dim(species.clean)
60 
## dataWD uses idCIRAD as identifier, but this is not a unique identifier in species.clean!
61 
## But wood density seems to also be available from seed.traits
62
63 
### need to read the different traits data based and merge .....
64 
bridge <- read.csv("./data/raw/DataParacou/BridgeDATA.g.csv",stringsAsFactors=FALSE, header = T, sep = ";")
65 
bridge$Latin_name <- paste(bridge[["Genus"]],bridge[["species"]],sep=" ")
66 
### check % of match of the bridg data
67 
sum(species.clean[["Latin_name"]] %in% bridge[["Latin_name"]])/length(species.clean[["Latin_name"]])
68 
## only 307 species /775 are in teh traits data ....
69
70 
seed.traits <- read.csv("./data/raw/DataParacou/Autour-de-Paracou-Releves-par-trait-et-taxon.txt",stringsAsFactors=FALSE, header = T, sep = "\t")
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
## Reformat seed.traits to one row per species, with each trait as a column spp.means <- (cast(seed.traits, LIB_TAXON ~ METHO_LIB, value = "MEASURE", fun = mean)) colnames(spp.means)[-1] <- paste(colnames(spp.means)[-1],".mean",sep="") spp.sds <- (cast(seed.traits, LIB_TAXON ~ METHO_LIB, value = "MEASURE", fun = sd)) colnames(spp.sds) <- paste(colnames(spp.sds),".sd",sep="") seed.traits2 <- cbind(spp.means,spp.sds[,-1])[,c("LIB_TAXON","Leaf nitrogen concentration (standard).mean","Leaf nitrogen concentration (standard).sd", "Specific leaf area (standard).mean","Specific leaf area (standard).sd", "Wood density.mean","Wood density.sd")] colnames(seed.traits2)[1] <- c("Latin_name") seed.traits2 <- seed.traits2[order(seed.traits2$Latin_name),] for(k in 2:ncol(seed.traits2)) seed.traits2[,k][!is.finite(seed.traits2[,k])] <- NA seed.traits2$sp <- species.clean$sp[match(seed.traits2$Latin_name,species.clean$Latin_name)] colnames(seed.traits2) <- c("Latin_name","Leaf.N.mean","Leaf.N.sd","SLA.mean","SLA.sd","Wood.density.mean","Wood.density.sd","sp") seed.traits2$SLA.mean <- seed.traits2$SLA.mean/1 ## Conversion from m2/kg to mm2/mg seed.traits2$Wood.density.mean <- seed.traits2$Wood.density.mean/1 ## conversion from g/cm3 to mg/mm3 # - $sp$ the species code as in previous table # - $Latin\_name$ the latin name of the species # - $Leaf.N.mean$ Leaf Nitrogen per mass in TRY mg/g # - $Seed.mass.mean$ dry mass in TRY mg # - $SLA.mean$ in TRY mm2 mg-1 # - $Wood.density.mean$ in TRY mg/mm3 # - $Max.height.mean$ from NFI data I compute the 99% quantile in m # - and the same columns with $sd$ instead of $mean$ with either the mean sd within species if species mean or the mean sd with genus if genus mean because no species data # - a dummy variable with true or false if genus mean ## still to be completed ############################################ FORMAT INDIVIDUAL TREE DATA data.paracou2 <- data.paracou[rep(1:nrow(data.paracou),each=2),c(1:10,(ncol(data.paracou)-2):ncol(data.paracou))] rownames(data.paracou2) <- 1:nrow(data.paracou2); data.paracou2 <- as.data.frame(data.paracou2) data.paracou2$yr1 <- rep(c(2001,2001+4),nrow(data.paracou)); data.paracou2$yr2 <- rep(c(2005,2005+4),nrow(data.paracou)) data.paracou2$year <- rep(c(4,4),nrow(data.paracou)) data.paracou2$dbh1 <- c(rbind(data.paracou$circum2001/pi,data.paracou$circum2005/pi)) data.paracou2$dbh2 <- c(rbind(data.paracou$circum2005/pi,data.paracou$circum2009/pi)) data.paracou2$code1 <- c(as.numeric(rbind(data.paracou$code2001,data.paracou$code2005))) data.paracou2$code2 <- c(as.numeric(rbind(data.paracou$code2005,data.paracou$code2009))) data.paracou2$dead <- rep(0,nrow(data.paracou)*2) data.paracou2$dead[c(as.numeric(data.paracou[["yeardied"]]) %in% 2002:2005 & (!is.na(data.paracou[["yeardied"]])), as.numeric(data.paracou[["yeardied"]]) %in% 2006:2009 & (!is.na(data.paracou[["yeardied"]])))] <- 1 ## remove tree dead at first census for both date (census 2001-2005 2005-2009) data.paracou <- subset(data.paracou2,subset=!(data.paracou2[['yr1']] ==2005 & (as.numeric(data.paracou[["yeardied"]]) %in% 2002:2005 & (!is.na(data.paracou[["yeardied"]]))))) ## change unit and names of variables to be the same in all data for the tree data.paracou$G <- 10*(data.paracou$dbh2-data.paracou$dbh1)/data.paracou$year ## diameter growth in mm per year data.paracou$G[data.paracou$code1>0] <- NA ## indivs with code indicating problem in dbh measurment at dbh1 data.paracou$G[data.paracou$code2>0] <- NA ## indivs with code indicating problem in dbh measurment at dbh2 data.paracou2$cenus <- data.paracou$yr1 data.paracou$D <- data.paracou[["dbh1"]]; data.paracou$D[data.paracou$D == 0] <- NA ;## diameter in cm data.paracou$plot <- data.paracou$plot#apply(data.paracou[,c("forest","plot","subplot")],1,paste,collapse=".") ## plot code data.paracou$htot <- rep(NA,length(data.paracou[["G"]])) ## height of tree in m data.paracou$obs.id <- 1:nrow(data.paracou) ### delete recruit in 2001 or 2005 for first census data.paracou <- subset(data.paracou,subset=!is.na(data.paracou$D)) ## minimum circumfer 30 delete all tree with a dbh <30/pi, data.paracou <- subset(data.paracou,subset= data.paracou[["D"]]>(30/pi)) ######################## ECOREGION ## paracou has only 1 eco-region ######################## PERCENT DEAD - compute numer of dead per plot to remove plot with disturbance perc.dead <- tapply(data.paracou[["dead"]],INDEX=data.paracou[["plot"]],FUN=function.perc.dead2) data.paracou <- merge(data.paracou,data.frame(plot=names(perc.dead),perc.dead=perc.dead), by = "plot", sort=FALSE)
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199
################################################## VARIABLES SELECTION FOR THE ANALYSIS #vec.abio.var.names <- c("MAT","MAP") ## MISSING NEED OTHER BASED ON TOPOGRAPHY ASK BRUNO vec.basic.var <- c("obs.id","tree.id","sp","plot","D","G","dead","census","year","htot","x","y","perc.dead") data.tree <- subset(data.paracou,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 ## NEED TO COMPUTE BASED ON RADIUS AROUND TARGET TREE ### species as factor because number data.tree[['sp']] <- factor(data.tree[['sp']]) Rlim <- 15 # set size of neighborhood for competition index ## for each census compute competition index ## FOR 2001 data.tree1 <- subset(data.tree,subset=data.tree[["census"]]==2001) library(doParallel) list.BA.SP.data1 <- mclapply(unique(data.tree1[['plot']]),FUN=fun.compute.BA.SP.XY.per.plot,data.tree=data.tree1,Rlim=Rlim,mc.cores=4) data.BA.sp1 <- rbind.fill(list.BA.SP.data1) ## FOR 2005 data.tree2 <- subset(data.tree,subset=data.tree[["census"]]==2005) library(doParallel) list.BA.SP.data2 <- mclapply(unique(data.tree2[['plot']]),FUN=fun.compute.BA.SP.XY.per.plot,data.tree=data.tree2,Rlim=Rlim,mc.cores=4) data.BA.sp2 <- rbind.fill(list.BA.SP.data2) data.tree <- rbind(data.tree1,data.tree2) data.BA.sp <- rbind(data.BA.sp1,data.BA.sp2) rm(data.tree1,data.tree2,data.BA.sp1,data.BA.sp2) ### TEST DATA FORMAT if(sum(! rownames(BA.SP.temp)==data.tree.s[['obs.id']]) >0) stop('rows not in the good order') if(sum(!colnames(BA.SP.temp)==as.character((levels(data.tree.s[['sp']]))))>0) stop('colnames does mot match species name') ## test same order as data.tree if(sum(!data.BA.SP[["obs.id"]] == data.tree[["obs.id"]]) >0) stop("competition index not in the same order than data.tree") ################################################ ## REMOVE TREE IN BUFFER ZONE BUFFER ZONE not.in.buffer.zone <- (data.tree[['x']]<(250-Rlim) & data.tree[['x']]>(0+Rlim) & data.tree[['y']]<(250-Rlim) & data.tree[['y']]>(0+Rlim)) # remove subset data.tree <- subset(data.tree,subset=not.in.buffer.zone) data.BA.sp <- subset(data.BA.sp,subset=not.in.buffer.zone) ## ## save everything as a list ## list.paracou <- list(data.tree=data.tree,data.BA.SP=data.BA.sp,data.traits=data.traits) ## save(list.spain,file="./data/process/list.paracou.Rdata")