TRY.R 10.97 KiB
########################################################
########################################################
###### READ TRY AND FORMAT DATA CHECK ERROR
################
#### use AccSpeciesName because not author name
source("./R/FUN.TRY.R")
library(MASS)
library(doParallel)
library(mvoutlier)
## read TRY data
TRY.DATA <- read.table("./data/raw/DataTRY/TRY_Proposal_177_DataRelease_2013_04_01.txt",
                       sep = "\t",header=TRUE,na.strings="", stringsAsFactors=FALSE)
TRY.DATA2 <- read.table("./data/raw/DataTRY/TRY_Proposal_177_DataRelease_2013_07_23.txt",
                       sep = "\t",header=TRUE,na.strings="", stringsAsFactors=FALSE)
### combine both data set
TRY.DATA <- rbind(TRY.DATA,TRY.DATA2)
rm(TRY.DATA2)
##################################
### ERROR FOUND IN THE DATA BASE
########################
### problem with the seed mass of this obs seed mass = 0 DELETE
TRY.DATA <- TRY.DATA[!(TRY.DATA$ObservationID==1034196 & TRY.DATA$DataName=="Seed dry mass"),]
#### IS "Quercuscrispla sp" an error standing for Quercus crispula synonym of Quercus mongolica subsp. crispula (Blume) Menitsky ? ask Jens
## TRY.DATA[TRY.DATA$AccSpeciesName=="Quercuscrispla sp" ,]
########################
########################
### first create a table with one row per Observation.id and column for each traits and variable
Non.Trait.Data <- c("Latitude", "Longitude", "Reference", "Date of harvest / measurement",
"Altitude", "Mean annual temperature (MAT)","Mean sum of annual precipitation (PPT)",
  "Plant developmental status / plant age","Maximum height reference",
  "Source in Glopnet",  "Number of replicates", "Sun vers. shade leaf qualifier" )
Trait.Data <- sort(names(((table(TRY.DATA$TraitName)))))
##########################
#### REFORMAT DATA from TRY
registerDoParallel(cores=5) ## affect automaticaly half of the core detected to the foreach here I decide to affect 4 cores
getDoParWorkers() ## here 8 core so 4 core if want to use more registerDoParallel(cores=6)
 TRY.DATA.FORMATED <- foreach(ObservationID.t=unique(TRY.DATA$ObservationID), .combine=rbind) %dopar%
            fun.extract.try(ObservationID.t,data=TRY.DATA,Non.Trait.Data,Trait.Data)
## head(TRY.DATA.FORMATED) 
## dim(TRY.DATA.FORMATED) 
saveRDS(TRY.DATA.FORMATED,file="./data/process/TRY.DATA.FORMATED.rds")
########################
########## READ RDS
TRY.DATA.FORMATED <- readRDS("./data/process/TRY.DATA.FORMATED.rds")
## TRY.DATA.FORMATED[TRY.DATA.FORMATED$ObservationID==1034196,"StdValue.Seed.mass"] <- NA
## head(TRY.DATA.FORMATED)
### export species list to check on tnrs web site
species.TRY <- (unique(TRY.DATA.FORMATED[["AccSpeciesName"]]))
write.csv(as.matrix(species.TRY),file="./data/process/species.TRY.csv",row.names=FALSE)
#######
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
## read data from TNRS tnrs.TRY <- read.delim("./output/tnrs_results.TRY.txt",sep="\t", na.strings="",stringsAsFactors=FALSE,header=TRUE) head(tnrs.TRY) fix(tnrs.TRY) #################### #################### ## COMPUTE MEAN AND SD FOR SPECIES from FRENCH NFI for 6 key traits key.main.traits2 <- c("StdValue.Leaf.nitrogen..N..content.per.dry.mass", "StdValue.Seed.mass", "StdValue.Leaf.specific.area..SLA.", "StdValue.Stem.specific.density..SSD.", "StdValue.Stem.conduit.area..vessel.and.tracheid.", "StdValue.Leaf.lifespan") ############################### ############################## ## READ CSV TABLE WITH LATIN NAME and CODE FOR FRENCH NFI DATA species.tab <- read.csv("./data/species.list/species.csv",sep="\t") species.tab2 <- species.tab[!is.na(species.tab$Latin_name),] rm(species.tab) gc() ### species IFN reformat names ## clean species names and synonyme names species.tab2$Latin_name <- (gsub("_", " ", species.tab2$Latin_name)) species.tab2$Latin_name_syn<- (gsub("_", " ", species.tab2$Latin_name_syn)) ## THIS TABLE HAS ALREADY THE SYNONYME FOR THE FRENCH SPECIES ## remove trailing white space species.tab2$Latin_name_syn<- trim.trailing(species.tab2$Latin_name_syn) ## create vector of species name species.IFN <- unique(pecies.tab2$Latin_name ) ## ###################################################################################### ## ####################################################################################### ## #### CHECKING SPECIES NAME TO GET ALL SYNONYMES ## ## NOT DONE YET !!!!!!!!!!!!!!!! ## ### export name to check in http://tnrs.iplantcollaborative.org/quick_start.html ## old.names <- unique(species.tab2$Latin_name) ## write.csv(as.matrix(old.names),file="./data/process/old.names.csv",row.names=FALSE) ## ## need to remove first raw with column name to submit to teh website ## ## read data from TNRS iPLANT ## tnrs.FRANCE <- read.delim("./output/tnrs_results.IFN.FRANCE.txt",sep="\t", na.strings="",stringsAsFactors=FALSE,header=TRUE) ## (cbind(test.tnrs$Name_submitted,test.tnrs$Accepted_name_species))[test.tnrs$Selected=="true",] ## ## need to do the same for TRY to have same match ## ## ACCORDING TO WILL THE BEST SOURCE IS http://www.theplantlist.org/ BUT NOT EASY TO ACCESS AND NOT WITH SYNO ########################################################################### ########################################################################### ##### EXTRACT SPECIES MEAN AND SD ### change format try species names TRY.DATA.FORMATED$AccSpeciesName <- as.character(TRY.DATA.FORMATED$AccSpeciesName) #### extract mean and sd per species without experimental data and detection of outlier when enough data or include experimental data. If no data compute mean of genus. ### The detection of outlier is based on the method in Kattge et al. 2011 only for univariate outlier. res.list <- lapply(species.IFN,FUN=fun.species.traits,species.table=species.tab2,traits=key.main.traits2,data=TRY.DATA.FORMATED) names(res.list) <- species.IFN ##### TRANSFORM LIST IN A TABLE
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
data.mean <- t(sapply(species.IFN,FUN=function(i,res.list) res.list[[i]]$mean ,res.list=res.list)) data.sd <- t(sapply(species.IFN,FUN=function(i,res.list) res.list[[i]]$sd,res.list=res.list)) data.exp <- t(sapply(species.IFN,FUN=function(i,res.list) res.list[[i]]$exp ,res.list=res.list)) data.genus <- t(sapply(species.IFN,FUN=function(i,res.list) res.list[[i]]$genus ,res.list=res.list)) data.nobs <- t(sapply(species.IFN,FUN=function(i,res.list) res.list[[i]]$nobs ,res.list=res.list)) #### data frame with all variables data.ifn.species.try.noout <- data.frame(data.mean,data.sd,data.exp,data.genus,data.nobs) names(data.ifn.species.try.noout) <- c(paste(c("Leaf.N","Seed.mass","SLA","Wood.Density" ,"Vessel.area","Leaf.Lifespan"),"mean",sep=".") ,paste(c("Leaf.N","Seed.mass","SLA","Wood.Density" ,"Vessel.area","Leaf.Lifespan"),"sd",sep=".") ,paste(c("Leaf.N","Seed.mass","SLA","Wood.Density" ,"Vessel.area","Leaf.Lifespan"),"exp",sep=".") ,paste(c("Leaf.N","Seed.mass","SLA","Wood.Density" ,"Vessel.area","Leaf.Lifespan"),"genus",sep=".") ,paste(c("Leaf.N","Seed.mass","SLA","Wood.Density" ,"Vessel.area","Leaf.Lifespan"),"nobs",sep=".")) saveRDS(data.ifn.species.try.noout ,file="./data/process/data.ifn.species.try.noout.rds") #### data.ifn.species.try.noout <- readRDS("./data/process/data.ifn.species.try.noout.rds") ##################################################################### #### assume that the SD is equal mean SD species ######## ## The table 5 in Kattge et al. 2011 GCB provides estimation of mean species sd ### SLA species sd log 0.09 ### Nmass species sd log 0.08 ### Seed Mass sd log 0.13 ## # SEE ## sd.log.SLA <- 0.09 ### based on Kattge et al. 2011 ## sd.log.Nmass <- 0.08 ### based on Kattge et al. 2011 ## sd.log.Seed.Mass <- 0.13 ### based on Kattge et al. 2011 ## sd.log.LL <- 0.03 ### based on Kattge et al. 2011 ###################### ### Computed sd over the data in log10 we have under the assumption sd constant over species with lm per species sd.vec.sp <- rep(NA,6) for(i in 1:length(key.main.traits2)){ eval(parse(text=paste("lm.obj <-lm(log10(TRY.DATA.FORMATED$",key.main.traits2[i],")~TRY.DATA.FORMATED$AccSpeciesName)"))) sd.vec.sp[i] <- sd(residuals(lm.obj)) } ## higher than tthe one reported in Kattge et al. 2011 #######################333 ### Computed sd over teh data we have under the assumption sd constant over genus sd.vec.genus <- rep(NA,6) for(i in 1:length(key.main.traits2)){ eval(parse(text=paste("lm.obj <-lm(log10(TRY.DATA.FORMATED$",key.main.traits2[i],")~sapply(TRY.DATA.FORMATED$AccSpeciesName,fun.get.genus))"))) sd.vec.genus[i] <- sd(residuals(lm.obj)) } ###################################################################################################### ### add columns with mean sd per species or per genus depending on whether species or genus data ## create vector nobs.names <- paste(c("Leaf.N","Seed.mass","SLA","Wood.Density",
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
"Vessel.area","Leaf.Lifespan") ,"nobs",sep=".") sd.names <- paste(c("Leaf.N","Seed.mass","SLA","Wood.Density","Vessel.area","Leaf.Lifespan") ,"sd",sep=".") genus.names <- paste(c("Leaf.N","Seed.mass","SLA","Wood.Density", "Vessel.area","Leaf.Lifespan") ,"genus",sep=".") names(sd.vec.sp) <- c("sdlog10.sp.Nmass","sdlog10.sp.Seed.Mass","sdlog10.sp.SLA", "sdlog10.sp.WD","sdlog10.sp.Vessel","sdlog10.sp.LL") names(sd.vec.genus) <- c("sdlog10.gs.Nmass","sdlog10.gs.Seed.Mass","sdlog10.gs.SLA", "sdlog10.gs.WD","sdlog10.gs.Vessel","sdlog10.gs.LL") ## save mean species and genus sd saveRDS(sd.vec.sp,file="./data/process/sd.vec.sp.rds") saveRDS(sd.vec.genus,file="./data/process/sd.vec.genus.rds") sd.vec.sp <- readRDS(file="./data/process/sd.vec.sp.rds") sd.vec.genus <- readRDS(file="./data/process/sd.vec.genus.rds") #### add column with the mean sd species or genus data.TRY.sd.update <- data.frame(data.ifn.species.try.noout, data.ifn.species.try.noout[,sd.names]) sd.names.1 <- paste(sd.names,1,sep=".") for (i in 1:length(sd.names.1)){ data.TRY.sd.update[[sd.names.1[i]]][!data.TRY.sd.update[[genus.names[i]]]] <- sd.vec.sp[i] data.TRY.sd.update[[sd.names.1[i]]][data.TRY.sd.update[[genus.names[i]]]] <- sd.vec.genus[i] } head(data.TRY.sd.update,10) saveRDS(data.TRY.sd.update,file="./data/process/data.TRY.sd.update.rds") ### # plot sd to show mark pdf("./figs/sd.traits.pdf") r <- barplot(sd.vec.sp ,names.arg=c("Leaf.N","SM","SLA","WD","Vessel","LL"),las=2,ylim=c(0,0.9),ylab="sd log10") points(r[,1],sd.vec.genus,col="red",pch=16,cex=2) ## for (i in 1:length(nobs.names)){ ## ## sd.obs <- data.TRY.sd.update[[sd.names[i]]][!data.TRY.sd.update[[genus.names[i]]]] ## ## points(rep(r[i,1],length(sd.obs)),sd.obs) ## ## sd.obs <- data.TRY.sd.update[[sd.names[i]]][data.TRY.sd.update[[genus.names[i]]]] ## ## points(rep(r[i,1],length(sd.obs)),sd.obs,col="red",pch=4) ## print(sd.obs) ## } dev.off()
281282283284