Commit 3a5a957d authored by Georges Kunstler's avatar Georges Kunstler
Browse files

start adding leaftype

parent 1c7099dc
......@@ -138,11 +138,33 @@ if(std) { data.tree.tot <- read.csv(file.path(base.dir, set,ecoregion,"data.tr
fun.data.for.lmer(data.tree.tot,trait,type.filling=type.filling)
}
## remove outlier following Condit approach
fun.remove.outlier.Condit <- function(data.tree){
data.tree[['G']][!(trim.positive.growth(data.tree[['G']]) &
trim.negative.growth(dbh1=data.tree[['D']]*10,
dbh2=data.tree[['D']]*10 +
data.tree[['year']]*data.tree[['G']]))] <- NA
data.tree[['BA.G']][!(trim.positive.growth(data.tree[['G']]) &
trim.negative.growth(data.tree[['D']]*10,dbh2=data.tree[['D']]*10 +
data.tree[['year']]*data.tree[['G']]))] <- NA
return(data.tree)
}
fun.select.data.for.analysis <- function(data.tree,abs.CWM.tntf,perc.neigh,BATOT,min.obs,
min.perc.neigh=0.90,min.BA.G=-50,max.BA.G=150){
min.perc.neigh=0.90,min.BA.G=-60,max.BA.G=500){
## remove tree with NA
data.tree <- subset(data.tree,subset=(!is.na(data.tree[["BA.G"]])) &
(!is.na(data.tree[["D"]])) )
## remove outlier
data.tree <- fun.remove.outlier.Condit(data.tree)
## remove tree with zero
data.tree <- subset(data.tree,subset=data.tree[["BA.G"]]>min.BA.G & data.tree[["BA.G"]]<max.BA.G &
data.tree[["D"]]>9.9)
......@@ -162,7 +184,7 @@ return((x-mean(x))/sd(x))
}
### get variables for lmer
fun.get.the.variables.for.lmer.tree.id <- function(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf,min.BA.G=50){
fun.get.the.variables.for.lmer.tree.id <- function(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf,min.BA.G=60){
logG <- fun.center.and.standardized.var(log(data.tree[["BA.G"]]+min.BA.G))
logD <- log(data.tree[["D"]])
species.id <- unclass(factor(data.tree[["sp"]]))
......@@ -187,7 +209,7 @@ return(data.frame(logG=logG,
sumBn=sumBn))
}
fun.get.the.variables.for.lmer.no.tree.id <- function(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf,min.BA.G=50){
fun.get.the.variables.for.lmer.no.tree.id <- function(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf,min.BA.G=60){
logG <- fun.center.and.standardized.var(log(data.tree[["BA.G"]]+min.BA.G))
logD <- log(data.tree[["D"]])
species.id <- unclass(factor(data.tree[["sp"]]))
......@@ -236,3 +258,42 @@ lmer.data <- fun.get.the.variables.for.lmer.no.tree.id(data.tree,BATOT,CWM.tn,ab
## remove growth outliers
##' .. remove too negative growth based on Condit R package with param fitted to BCI ..
##'
##' .. taken from trim.growth function in CTFS.R ..
##' @title trim.negative.growth
##' @param dbh1 in mm
##' @param dbh2 in mm
##' @param slope not to be changed
##' @param intercept
##' @param err.limit
##' @return a vector TRUE FALSE with FALSE outlier to be removed
##' @author Kunstler
trim.negative.growth <- function(dbh1,dbh2,slope=0.006214,
intercept=.9036,err.limit=5){
stdev.dbh1 <- slope*dbh1+intercept
bad.grow <- which(dbh2<=(dbh1-err.limit*stdev.dbh1))
accept <- rep(TRUE,length(dbh1))
accept[bad.grow] <- FALSE
accept[is.na(dbh1) | is.na(dbh2) | dbh2<=0 | dbh1<=0] <- FALSE
return(accept)
}
##' .. remove too high growth ..
##'
##' .. taken from trim.growth in Condit CTFS R package ..
##' @title trim.positive.growth
##' @param growth in mm
##' @param maxgrow in mm
##' @return TRUE FALSE vector with FALSE outlier
##' @author Kunstler
trim.positive.growth <- function(growth,maxgrow=75){
bad.grow <- which(growth>maxgrow)
accept <- rep(TRUE,length(growth))
accept[bad.grow] <- FALSE
accept[is.na(growth)] <- FALSE
return(accept)
}
......@@ -100,7 +100,7 @@ data.tree.tot <- read.csv(file=file.path(base.dir, "data.all.csv"),
}
fun.select.data.for.analysis <- function(data.tree,abs.CWM.tntf,perc.neigh,BATOT,min.obs,
min.perc.neigh=0.90,min.BA.G=-50,max.BA.G=150){
min.perc.neigh=0.90,min.BA.G=-60,max.BA.G=500){
## remove tree with NA
data.tree <- subset(data.tree,subset=(!is.na(data.tree[["BA.G"]])) &
(!is.na(data.tree[["D"]])) )
......@@ -129,7 +129,7 @@ return((x-mean(x))/sd(x))
### get variables for lmer
fun.get.the.variables.for.lmer.tree.id <- function(data.tree,BATOT, CWM.tn,
abs.CWM.tntf, tf,
min.BA.G=50){
min.BA.G=60){
logG <- fun.center.and.standardized.var(log(data.tree[["BA.G"]]+min.BA.G))
logD <- fun.center.and.standardized.var(log(data.tree[["D"]]))
species.id <- factor(data.tree[["sp"]])
......
......@@ -9,7 +9,7 @@ species.clean <- data.frame(sp=data.tree[!duplicated(data.tree[["sp"]]),"sp"],
Latin_name=data.tree[!duplicated(data.tree[["sp"]]),"sp.name"],
Latin_name_syn=data.tree[!duplicated(data.tree[["sp"]]),"sp.name"],
stringsAsFactors =FALSE)
rm(data.tree)
## read in data
data.TRY.std <- readRDS("output/formatted/TRY/data.TRY.std.rds")
......@@ -65,6 +65,7 @@ data.cat.extract <- do.call("rbind",lapply(data.traits$sp ,fun.get.cat.var.from.
data.cat.extract <- fun.change.factor.pheno.try(data.cat.extract)
data.cat.extract <- fun.change.factor.angio.try(data.cat.extract)
data.cat.extract <- fun.fill.pheno.try.with.zanne(data.cat.extract)
data.cat.extract <- fun.change.factor.leaftype.try(data.cat.extract)
## fix pheno for species with issue
data.cat.extract[data.cat.extract$Latin_name %in% c('Betula spp.','Crataegus spp.','Fraxinus spp.',
......
......@@ -291,6 +291,8 @@ fun.get.cat.var.from.try <- function(sp,data.traits,try.cat,Pheno.Zanne){
try.cat$AccSpeciesName,"PhylogeneticGroup"][1],
Pheno.T=try.cat[data.traits[data.traits$sp==sp,"Latin_name"] ==
try.cat$AccSpeciesName,"LeafPhenology"][1],
LeafType.T=try.cat[data.traits[data.traits$sp==sp,"Latin_name"] ==
try.cat$AccSpeciesName,"LeafType"][1],
Pheno.Z=Pheno.Zanne[data.traits[data.traits$sp==sp,"Latin_name"] ==
Pheno.Zanne$Binomial,'Phenology'][1] ,
stringsAsFactors=FALSE)
......@@ -301,6 +303,8 @@ fun.get.cat.var.from.try <- function(sp,data.traits,try.cat,Pheno.Zanne){
try.cat$AccSpeciesName,"PhylogeneticGroup"][1],
Pheno.T=try.cat[data.traits[data.traits$sp==sp,"Latin_name"] ==
try.cat$AccSpeciesName,"LeafPhenology"][1],
LeafType.T=try.cat[data.traits[data.traits$sp==sp,"Latin_name"] ==
try.cat$AccSpeciesName,"LeafType"][1],
Pheno.Z=NA ,
stringsAsFactors=FALSE)
}
......@@ -311,6 +315,7 @@ fun.get.cat.var.from.try <- function(sp,data.traits,try.cat,Pheno.Zanne){
Latin_name=data.traits[data.traits$sp==sp,"Latin_name"],
Phylo.group=NA,
Pheno.T=NA,
LeafType.T=NA,
Pheno.Z=Pheno.Zanne[data.traits[data.traits$sp==sp,"Latin_name"] ==
Pheno.Zanne$Binomial,'Phenology'][1] ,
stringsAsFactors=FALSE)
......@@ -319,11 +324,19 @@ fun.get.cat.var.from.try <- function(sp,data.traits,try.cat,Pheno.Zanne){
Latin_name=data.traits[data.traits$sp==sp,"Latin_name"],
Phylo.group=NA,
Pheno.T=NA,
LeafType.T=NA,
Pheno.Z=NA ,
stringsAsFactors=FALSE)
}
}
return(data.res)
## if missing value for Phylo.group check genus
if(is.na(data.res$Phylo.group)) {
genus <- sub(" .*", "", data.traits[data.traits$sp==sp,"Latin_name"])
genus.vec <- sub(" .*", "", try.cat$AccSpeciesName)
data.res$Phylo.group <- try.cat[genus == genus.vec,"PhylogeneticGroup"][1]
}
return(data.res)
}
......@@ -349,6 +362,15 @@ data.cat.extract[data.cat.extract$Phylo.group=='Angiosperm_Monocotyl'
return(data.cat.extract)
}
fun.change.factor.leaftype.try <- function(data.cat.extract){
data.cat.extract[data.cat.extract$LeafType.T=='broadleaved'
& !is.na(data.cat.extract$LeafType.T),'LeafType.T'] <- 'broadleaved'
data.cat.extract[data.cat.extract$LeafType.T!='broadleaved'
& !is.na(data.cat.extract$LeafType.T),'LeafType.T'] <- 'No.broadleaved'
return(data.cat.extract)
}
# fill missing TRY deciduous/evergreen by Zanne
fun.fill.pheno.try.with.zanne <- function(data.cat.extract){
data.cat.extract[is.na(data.cat.extract$Pheno.T),'Pheno.T'] <- data.cat.extract[is.na(data.cat.extract$Pheno.T),'Pheno.Z']
......
......@@ -87,6 +87,8 @@ data.japan <- merge(data.japan, data.frame(plot = site.data$PlotID, Lat = site.d
## x11()
## plot(site.data$Temp30yr,site.data$Rain30yr)
################################## FORMAT INDIVIDUAL TREE DATA
data.japan$year <- as.numeric(difftime(data.japan$date2, data.japan$date1, units = "weeks")/52) ## number of year between measurement
data.japan$G <- 10*(data.japan[["dbh2"]] - data.japan[["dbh1"]])/data.japan$year ## diameter growth in mm per year
......@@ -130,14 +132,14 @@ make.quad <- do.call("rbind",lapply(unique(data.japan$cluster),
data.japan <- merge(data.japan, make.quad,by.x='obs.id',by.y="tree.id")
## #plot each plot
## ## #plot each plot
## pdf("figs/plots.japan.pdf")
## lapply(unique(data.japan[["cluster"]]),FUN=fun.circles.plot,
## data.japan[['x']],data.japan[['y']],
## data.japan[["cluster"]],data.japan[["D"]],
## inches=0.2,fg.l=data.japan$make.quad)
## inches=0.1,fg.l=data.japan$make.quad)
## dev.off()
#
## #
data.japan$plot <- data.japan$make.quad
......
......@@ -2,6 +2,8 @@
### MERGE Mbaiki DATA
data.mbaiki$D <- data.mbaiki[["dbh1"]];
data.mbaiki$D[data.mbaiki$D == 0] <- NA;
rm(list = ls())
# setwd("/home/ghislain/Documents/Ghislain-CIRAD/Traits_Competition_Georges/trait.competition.workshop")
......@@ -51,12 +53,14 @@ data.mbaiki <- merge(data.mbaiki, make.quad,by='tree.id')
### read species names
species.clean <- read.csv("data/raw/Mbaiki/Mbaiki.traits.csv",stringsAsFactors=FALSE, header = T, sep = ",")
species.clean <- read.csv("data/raw/Mbaiki/Mbaiki.traits.csv",
stringsAsFactors=FALSE, header = T, sep = ",")
species.clean$sp.name <- species.clean$Species
#data.mbaiki$spname %in% species.clean$sp.name
## All species in data.mbaiki are in species.clean. Good to go!
### Based on the familly provided by the tnrs web site from teh species list seems to be no tree fern or palm
### Based on the familly provided by the tnrs web site from teh species list
### seems to be no tree fern or palm
### see files "data/raw/Mbaiki/tnrs_results_Mbaiki.txt" NEED to ckeck with ghislain
......@@ -67,19 +71,27 @@ data.mbaiki2 <- as.data.frame(data.mbaiki2)
data.mbaiki2$yr1 <- rep(c(1995,1995+5),nrow(data.mbaiki));
data.mbaiki2$yr2 <- rep(c(2000,2000+5),nrow(data.mbaiki))
data.mbaiki2$year <- rep(c(5,5),nrow(data.mbaiki))
data.mbaiki2$dbh1 <- c(rbind(data.mbaiki$circum1995/pi,data.mbaiki$circum2000/pi))
data.mbaiki2$dbh2 <- c(rbind(data.mbaiki$circum2000/pi,data.mbaiki$circum2005/pi))
data.mbaiki2$code1 <- c(as.numeric(rbind(data.mbaiki$code1995,data.mbaiki$code2000)))
data.mbaiki2$code2 <- c(as.numeric(rbind(data.mbaiki$code2000,data.mbaiki$code2005)))
data.mbaiki2$dbh1 <- c(rbind(data.mbaiki$circum1995/pi,
data.mbaiki$circum2000/pi))
data.mbaiki2$dbh2 <- c(rbind(data.mbaiki$circum2000/pi,
data.mbaiki$circum2005/pi))
data.mbaiki2$code1 <- c(as.numeric(rbind(data.mbaiki$code1995,
data.mbaiki$code2000)))
data.mbaiki2$code2 <- c(as.numeric(rbind(data.mbaiki$code2000,
data.mbaiki$code2005)))
data.mbaiki2$dead <- rep(0,nrow(data.mbaiki)*2)
data.mbaiki2$dead[c(as.numeric(data.mbaiki[["yeardied"]]) %in% 1995:2000 & (!is.na(data.mbaiki[["yeardied"]])),
as.numeric(data.mbaiki[["yeardied"]]) %in% 2000:2005 & (!is.na(data.mbaiki[["yeardied"]])))] <- 1
data.mbaiki2$dead[c(as.numeric(data.mbaiki[["yeardied"]]) %in% 1995:2000 &
(!is.na(data.mbaiki[["yeardied"]])),
as.numeric(data.mbaiki[["yeardied"]]) %in% 2000:2005 &
(!is.na(data.mbaiki[["yeardied"]])))] <- 1
## remove tree dead at first census for both date - as in dead before 1995?
data.mbaiki <- subset(data.mbaiki2, yeardied > 1995)
data.mbaiki <- data.mbaiki2[data.mbaiki2$yeardied > 1995 |
is.na(data.mbaiki2$yeardied),]
## change unit and names of variables to be the same in all data for the tree
data.mbaiki$G <- 10*(data.mbaiki$dbh2-data.mbaiki$dbh1)/data.mbaiki$year ## diameter growth in mm per year
data.mbaiki$G <- 10*(data.mbaiki$dbh2-data.mbaiki$dbh1)/data.mbaiki$year
## diameter growth in mm per year
## Ghislain's comment
## ------------------
......@@ -110,15 +122,22 @@ data.mbaiki$G <- 10*(data.mbaiki$dbh2-data.mbaiki$dbh1)/data.mbaiki$year ## diam
## TRUE;9;"mesure au relascope"
##
## So (1) for the competition (BA computation, all trees can be kept).
## (2) for a good measure of growth: remove trees with codes 1 and 9 but keep the others.
data.mbaiki$G[data.mbaiki$code1 %in% c(1,9)] <- NA ## indivs with code indicating problem in dbh measurment at dbh1
data.mbaiki$G[data.mbaiki$code2 %in% c(1,9)] <- NA ## indivs with code indicating problem in dbh measurment at dbh2
data.mbaiki$BA.G <- (pi*(data.mbaiki$dbh2/2)^2-pi*(data.mbaiki$dbh1/2)^2)/data.mbaiki$year ## BA growth in cm2 per year
data.mbaiki$BA.G[data.mbaiki$code1 %in% c(1,9)] <- NA ## indivs with code indicating problem in dbh measurment at dbh1
data.mbaiki$BA.G[data.mbaiki$code2 %in% c(1,9)] <- NA ## indivs with code indicating problem in dbh measurment at dbh2
## (2) for a good measure of growth: remove trees with codes 1 and 9
## but keep the others.
data.mbaiki$G[data.mbaiki$code1 %in% c(1,9)] <- NA
## indivs with code indicating problem in dbh measurment at dbh1
data.mbaiki$G[data.mbaiki$code2 %in% c(1,9)] <- NA
## indivs with code indicating problem in dbh measurment at dbh2
data.mbaiki$BA.G <- (pi*(data.mbaiki$dbh2/2)^2-pi*(data.mbaiki$dbh1/2)^2)/
data.mbaiki$year ## BA growth in cm2 pBer year
data.mbaiki$BA.G[data.mbaiki$code1 %in% c(1,9)] <- NA
## indivs with code indicating problem in dbh measurment at dbh1
data.mbaiki$BA.G[data.mbaiki$code2 %in% c(1,9)] <- NA
B## indivs with code indicating problem in dbh measurment at dbh2
data.mbaiki$census <- data.mbaiki$yr1
data.mbaiki$D <- data.mbaiki[["dbh1"]]; data.mbaiki$D[data.mbaiki$D == 0] <- NA ;## diameter in cm
data.mbaiki$D <- data.mbaiki[["dbh1"]];
data.mbaiki$D[data.mbaiki$D == 0] <- NA ;## diameter in cm
data.mbaiki$cluster <- paste("p",data.mbaiki$cluster,sep=".")
data.mbaiki$htot <- rep(NA,length(data.mbaiki[["G"]])) ## height of tree in m
data.mbaiki$obs.id <- 1:nrow(data.mbaiki)
......@@ -150,7 +169,8 @@ data.mbaiki$MAT <- clim$MAT
data.mbaiki$MAP <- clim$MAP
###################### PLOT SELECTION FOR THE ANALYSIS - NEEDS REDOING
vec.basic.var <- c("obs.id","tree.id", "sp", "sp.name", "cluster","plot", "ecocode", "D", "G","BA.G","year", "dead",
vec.basic.var <- c("obs.id","tree.id", "sp", "sp.name", "cluster","plot",
"ecocode", "D", "G","BA.G","year", "dead",
'Lon','Lat',"x", "y", "census",'MAT','MAP')
data.tree <- subset(data.mbaiki, select = c(vec.basic.var))
......@@ -166,7 +186,8 @@ write.csv(data.tree,file="output/formatted/Mbaiki/tree.csv",row.names = FALSE)
### write data plot with variables only at the plot level.
vec.basic.var.p <- c("plot", "cluster", "Lon","Lat","ecocode",'MAT','MAP')
data.plot <- subset(data.tree, subset=!duplicated(data.tree$cluster),select = c(vec.basic.var.p))
data.plot <- subset(data.tree, subset=!duplicated(data.tree$cluster),
select = c(vec.basic.var.p))
write.csv(data.plot,file="output/formatted/Mbaiki/plot.csv",row.names = FALSE)
......@@ -35,7 +35,7 @@ data.nz$year <- (data.nz[["t1"]] - data.nz[["t0"]]) ## number of year between m
data.nz$D <- data.nz[["D0"]] ## diameter in cm
data.nz$dead <- as.numeric(is.na(data.nz[["D1"]])) ## dummy variable for dead tree 0 alive 1 dead
data.nz$plot <- data.nz$plid; data.nz$plid <- NULL
data.nz$cluster <- rep(NA, nrow(data.nz)) ## no plot cluster
data.nz$cluster <- data.nz$surv ## no plot cluster
data.nz$htot <- rep(NA, nrow(data.nz)) ## Max height is already available so have as missing
data.nz$tree.id <- data.nz$tag
head(subset(data.nz,subset=data.nz$tag %in% names(table(data.nz$tag))[table(data.nz$tag)==2])) ### 3 individuals with two observation with strange value check with Daniel and David add an issue
......
......@@ -58,20 +58,20 @@ data.paracou <- subset(data.paracou,subset=(!is.na(data.paracou[["y"]])) &
data.paracou[["y"]]<251)
#### REMOVE PLOTs 16 17 18 ACCORDING TO GHSILAIN
data.paracou <- subset(data.paracou,
subset=! data.paracou[["cluster"]] %in% 17:18)
## keep only tree alive in 2001 and 1998 for cluster 16
VEC.select <- data.paracou$yeardied
VEC.select <- NA
VEC.select[data.paracou$cluster != 16] <-
!(data.paracou$yeardied[data.paracou$cluster != 16] <= 2001 &
!is.na(data.paracou$yeardied[data.paracou$cluster != 16]))
VEC.select[data.paracou$cluster == 16] <-
!(data.paracou$yeardied[data.paracou$cluster == 16] <= 1998 &
!is.na(data.paracou$yeardied[data.paracou$cluster == 16]))
subset=! data.paracou[["cluster"]] %in% 16:18)
## ## keep only tree alive in 2001 and 1998 for cluster 16
## VEC.select <- data.paracou$yeardied
## VEC.select <- NA
## VEC.select[data.paracou$cluster != 16] <-
## !(data.paracou$yeardied[data.paracou$cluster != 16] <= 2001 &
## !is.na(data.paracou$yeardied[data.paracou$cluster != 16]))
data.paracou <- subset(data.paracou,
subset=VEC.select)
## VEC.select[data.paracou$cluster == 16] <-
## !(data.paracou$yeardied[data.paracou$cluster == 16] <= 1998 &
## !is.na(data.paracou$yeardied[data.paracou$cluster == 16]))
## data.paracou <- subset(data.paracou,
## subset=VEC.select)
## plot each plot
## pdf("figs/plots.paracou.pdf")
......@@ -93,15 +93,15 @@ make.quad <- do.call("rbind",lapply(unique(data.paracou$cluster),
square.s=square.s.t))
data.paracou <- merge(data.paracou, make.quad,by='tree.id')
## Levels.cluster <- levels(as.factor(data.paracou$cluster))
## n.cluster <- length(Levels.cluster) # 10 clusters => 10 "plots"
Levels.cluster <- levels(as.factor(data.paracou$cluster))
n.cluster <- length(Levels.cluster) # 10 clusters => 10 "plots"
## for (i in Levels.cluster){
## dat <- subset(data.paracou, subset=data.paracou[["cluster"]]==i)
## x11()
## symbols(dat$x,dat$y,circles=dat$circum2001,inches=0.2,main=i,
## fg=unclass(factor(dat$make.quad)))
## }
for (i in Levels.cluster){
dat <- subset(data.paracou, subset=data.paracou[["cluster"]]==i)
x11()
symbols(dat$x,dat$y,circles=dat$circum2001,inches=0.2,main=i,
fg=unclass(factor(dat$make.quad)))
}
......@@ -141,10 +141,10 @@ 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$yr1[data.paracou2$cluster==16] <- rep(c(1998,1998+7),
nrow(data.paracou[data.paracou2$cluster==16,]))
data.paracou2$yr2[data.paracou2$cluster==16] <- rep(c(2005,2005+5),
nrow(data.paracou[data.paracou2$cluster==16,])) # For 16
## data.paracou2$yr1[data.paracou2$cluster==16] <- rep(c(1998,1998+7),
## nrow(data.paracou[data.paracou2$cluster==16,]))
## data.paracou2$yr2[data.paracou2$cluster==16] <- rep(c(2005,2005+5),
## nrow(data.paracou[data.paracou2$cluster==16,])) # For 16
data.paracou2$year <- rep(c(4,4),nrow(data.paracou))
data.paracou2$year[data.paracou2$cluster==16] <- rep(c(7,5),
......
......@@ -28,7 +28,6 @@ if(dim(data.all)[1] != sum(mat.perc[['num.obs']]))
#- pattern of CWM
#- pattern ED angio/conif
## look at BATOT per ecocode
MAP.ECOCODE <- tapply(data.all$MAP,INDEX=data.all$ecocode,FUN=mean)
SET.ECOCODE <- tapply(data.all$set,INDEX=data.all$ecocode,FUN=unique)
......@@ -37,13 +36,78 @@ boxplot(data.all$BATOT ~ data.all$ecocode, las = 3,
at = log(MAP.ECOCODE), boxwex = 0.05, outline = FALSE,
col=col.vec[SET.ECOCODE])
pdf("figs/test.processed/boxplot.BATOT.D.G.BA.G.pdf")
par(mfrow=c(2,3))
boxplot(data.all$BATOT ~ data.all$set, las = 3)
boxplot(data.all$D ~ data.all$set, las = 3)
boxplot(data.all$G ~ data.all$set, las = 3)
boxplot(data.all$BA.G ~ data.all$set, las = 3)
dev.off()
pdf("figs/test.processed/boxplot.traits.pdf")
par(mfrow=c(2,3))
boxplot(data.all$Leaf.N.focal ~ data.all$set, las = 3,main="Leaf.N")
boxplot(data.all$SLA.focal ~ data.all$set, las = 3,main="SLA")
boxplot(data.all$Wood.density.focal ~ data.all$set, las = 3,main="Wood.density")
boxplot(data.all$Seed.mass.focal ~ data.all$set, las = 3,main="Seed.mass")
boxplot(data.all$Max.height.focal ~ data.all$set, las = 3,main="Max.height")
dev.off()
par(mfrow=c(2,2))
boxplot(data.all$Leaf.N.perc.species ~ data.all$set, las = 3,main="Leaf.N")
boxplot(data.all$SLA.perc.species ~ data.all$set, las = 3,main="SLA")
boxplot(data.all$Wood.density.perc.species ~ data.all$set, las = 3,main="Wood.density")
boxplot(data.all$Seed.mass.perc.species ~ data.all$set, las = 3,main="Seed.mass")
par(mfrow=c(2,2))
boxplot(data.all$Leaf.N.perc.genus ~ data.all$set, las = 3,main="Leaf.N")
boxplot(data.all$SLA.perc.genus ~ data.all$set, las = 3,main="SLA")
boxplot(data.all$Wood.density.perc.genus ~ data.all$set, las = 3,main="Wood.density")
boxplot(data.all$Seed.mass.perc.genus ~ data.all$set, las = 3,main="Seed.mass")
### compute percentage of evergreen and conif per plot
### compute mean BATOT, number of species, traits and VAR OF TRAITS per cluster
system.time(data.summarise <- fun.compute.all.var.cluster(data.all))
### NEED TO CHECK WHY JAPAN REACH 300 of BATOT
## angio
boxplot(data.summarise$per.angio~data.summarise$set,las=3)
per.angio.ecocode.mean <- tapply(data.summarise$per.angio,data.summarise$ecocode,mean,na.rm=TRUE)
per.angio.ecocode.quant <- do.call('rbind',tapply(data.summarise$per.angio,data.summarise$ecocode,quantile,
probs=c(0.1,0.8),na.rm=TRUE))
names.ecocode <- tapply(data.summarise$ecocode,data.summarise$ecocode,unique)
names(per.angio.ecocode.mean) <- names.ecocode
par(mar=c(8.1,4.1,4.1,2.1))
mp <- barplot(per.angio.ecocode.mean,las=3)
segments(mp,per.angio.ecocode.quant[,1],
mp,per.angio.ecocode.quant[,2])
## deciduous
boxplot(data.summarise$per.deciduous~data.summarise$set,las=3)
per.deciduous.ecocode.mean <- tapply(data.summarise$per.deciduous,data.summarise$ecocode,mean,na.rm=TRUE)
per.deciduous.ecocode.quant <- do.call('rbind',tapply(data.summarise$per.deciduous,data.summarise$ecocode,quantile,
probs=c(0.1,0.8),na.rm=TRUE))
names.ecocode <- tapply(data.summarise$ecocode,data.summarise$ecocode,unique)
names(per.deciduous.ecocode.mean) <- names.ecocode
par(mar=c(8.1,4.1,4.1,2.1))
mp <- barplot(per.deciduous.ecocode.mean,las=3)
segments(mp,per.deciduous.ecocode.quant[,1],
mp,per.deciduous.ecocode.quant[,2])
par(mfrow=c(1,2))
plot(data.summarise$MAP,data.summarise$BATOT,
,col=col.vec[data.summarise$set],cex=0.1)
......@@ -51,9 +115,11 @@ fun.boxplot.breaks((data.summarise$MAP),data.summarise$BATOT,Nclass=15,add=TRUE)
legend("topright",legend=names(col.vec),col=col.vec,pch=1)
plot(log(data.summarise$MAP),data.summarise$BATOT,
,col=col.vec[data.summarise$set],cex=0.1)
fun.boxplot.breaks(log(data.summarise$MAP),data.summarise$BATOT,Nclass=15,add=TRUE)
fun.boxplot.breaks(log(data.summarise$MAP),data.summarise$BATOT,
Nclass=15,add=TRUE)
by(data.summarise,INDICES=data.summarise$set,function(data,col.vec) {x11();plot(data$MAP,data$BATOT,
by(data.summarise,INDICES=data.summarise$set,
function(data,col.vec) {x11();plot(data$MAP,data$BATOT,
main=unique(data$set),
col=col.vec[data$set])},col.vec)
......
......@@ -232,14 +232,14 @@ fun.test.value.one.ecoregion.B <- function(data.CWM,set,ecocode.select,path.form
#### compute description of each table
fun.perc.sup.0.9.sp <- function(i,var.name,var.names.perc,data){
return( sum(data[[var.names.perc[i]]]>0.9 &
fun.perc.sup.treshold.sp <- function(i,var.name,var.names.perc,data,treshold){
return( sum(data[[var.names.perc[i]]]>treshold &
!is.na(data[[var.names.perc[i]]]) &
!is.na(data[[var.name[i]]]))/length(data[[var.names.perc[i]]]))
}
fun.compute.percentage.species.genus <- function(data){
fun.compute.percentage.species.genus <- function(data,treshold=0.8){
perc.gymno <- sum( data[['Phylo.group']]=='Gymnosperm' &
!is.na(data[['Phylo.group']]))/sum(!is.na(data[['Phylo.group']]))
perc.ev <- sum(data[['Pheno.T']]=='EV'
......@@ -248,17 +248,17 @@ traits.name <- c("Leaf.N","Seed.mass","SLA","Wood.density","Max.height")
var.name.fill <- paste(traits.name,"abs.CWM.fill",sep=".")
var.name.perc.genus <- paste(traits.name,"perc.genus",sep=".")
var.name.perc.species <- paste(traits.name,"perc.species",sep=".")
perc0.9.non.missing.sp <- lapply(1:length(var.name.fill),fun.perc.sup.0.9.sp,var.name.fill,
var.name.perc.species,data)
perc0.9.non.missing.genus <- lapply(1:length(var.name.fill),fun.perc.sup.0.9.sp,var.name.fill,
var.name.perc.genus,data)
perc.T.non.missing.sp <- lapply(1:length(var.name.fill),fun.perc.sup.treshold.sp,var.name.fill,
var.name.perc.species,data,treshold=treshold)
perc.T.non.missing.genus <- lapply(1:length(var.name.fill),fun.perc.sup.treshold.sp,var.name.fill,
var.name.perc.genus,data,treshold=treshold)
perc.genus <- lapply(var.name.perc.genus,function(i,data) mean(data[[i]],na.rm=TRUE),data)
perc.species <- lapply(var.name.perc.species,function(i,data) mean(data[[i]],na.rm=TRUE),data)
vec.res <- c(nrow(data),perc.gymno,perc.ev,perc0.9.non.missing.sp,
perc0.9.non.missing.genus,
vec.res <- c(nrow(data),perc.gymno,perc.ev,percT.non.missing.sp,
percT.non.missing.genus,
perc.genus,perc.species)
names(vec.res) <- c('num.obs','perc.gymno','perc.ev',paste(traits.name,"perc.above0.9.sp",sep="."),
paste(traits.name,"perc.above0.9.genus",sep="."),
names(vec.res) <- c('num.obs','perc.gymno','perc.ev',paste(traits.name,"perc.above",threshold,"sp",sep="."),
paste(traits.name,"perc.above",treshold,"genus",sep="."),
var.name.perc.genus,var.name.perc.species)
return(vec.res)
}
......
......@@ -6,29 +6,18 @@ mkdir -p trait.workshop
for site in "'BCI'" "'Canada'" "'France'" "'Fushan'" "'NVS'" "'Paracou'" "'Spain'" "'US'" "'Swiss'" "'Sweden'" "'NSW'" "'Mbaiki'" "'Luquillo'" "'Japan'"; do
echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.nolog.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.1, run.lmer,type.filling='species',std=TRUE);print('done')\"" > trait.workshop/species1$site.sh
echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.nolog.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.1, run.lmer,type.filling='species',std=FALSE);print('done')\"" > trait.workshop/species1$site.sh
qsub trait.workshop/species1$site.sh -l nodes=1:ppn=1 -N "lmerspecies1$site" -q opt32G -j oe
echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.nolog.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.2, run.lmer,type.filling='species',std=TRUE);print('done')\"" > trait.workshop/species2$site.sh
echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.nolog.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.2, run.lmer,type.filling='species',std=FALSE);print('done')\"" > trait.workshop/species2$site.sh
qsub trait.workshop/species2$site.sh -l nodes=1:ppn=1 -N "lmerspecies2$site" -q opt32G -j oe
echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.nolog.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.1, run.lmer,type.filling='genus',std=TRUE);print('done')\"" > trait.workshop/genus1$site.sh
echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.nolog.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.1, run.lmer,type.filling='genus',std=FALSE);print('done')\"" > trait.workshop/genus1$site.sh
qsub trait.workshop/genus1$site.sh -l nodes=1:ppn=1 -N "lmergenus1$site" -q opt32G -j oe
echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.nolog.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.2, run.lmer,type.filling='genus',std=TRUE);print('done')\"" > trait.workshop/genus2$site.sh
echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.nolog.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.2, run.lmer,type.filling='genus',std=FALSE);print('done')\"" > trait.workshop/genus2$site.sh