Skip to content

GitLab

Commit db81af37 authored by Georges Kunstler's avatar Georges Kunstler
Browse files

new processimng of the data

parent 3b4e3196
Expand all Hide whitespace changes
Inline Side-by-side
Makefile
View file @ db81af37
......@@ -223,6 +223,15 @@ $(D6)/Done.txt: R/process.data/test.tree.CWM.R $(D5)/Done.txt $(D3Done)
#-------------------------------------------------------
#-------------------------------------------------------
MERGE.ALL: output/processed/data.all.csv
output/processed/data.all.csv: R/process.data/merge.all.processed.data.R
Rscript $<
#-------------------------------------------------------
# This susbtitution rule should work as rule, but not, why not?
# docs/output/formatted/%/tree.csv: $(D1)/%/* %.R
# Rscript %.R
......
R/analysis/lmer.run.nolog.all.R 0 → 100644
View file @ db81af37
###########################
###########################
### FUNCTION TO RUN LMER ESTIMATION WITH NO logBA for all in one big model
library(lme4)
run.models.for.set.all.traits <- function(model.file,fun.model, traits =
c("SLA", "Wood.density","Max.height","Leaf.N","Seed.mass"),type.filling, ...){
for(trait in traits)
run.multiple.model.for.set.one.trait(model.file,fun.model, trait, type.filling=type.filling,...)
}
run.multiple.model.for.set.one.trait <- function(model.files,fun.model, trait, type.filling, ecoregions = get.ecoregions.for.set(set), ...){
for (m in model.files)
try(run.model.for.set.one.trait (m, fun.model,trait, type.filling=type.filling, ...))
}
run.model.for.set.one.trait <- function(model.file,fun.model, trait, set,type.filling, ecoregions = get.ecoregions.for.set(set), ...){
fun.model <- match.fun(fun.model)
try(fun.model(model.file, trait, type.filling=type.filling,...))
}
#=====================================================================
# Run lmer() (in package lme4) for one ecoregion, one trait, one model
#=====================================================================
model.files.lmer.1 <- c("R/analysis/model.lmer.all/model.lmer.LOGLIN.E.R",
"R/analysis/model.lmer.all/model.lmer.LOGLIN.R.R",
"R/analysis/model.lmer.all/model.lmer.LOGLIN.ER.R")
model.files.lmer.2 <- c("R/analysis/model.lmer.all/model.lmer.LOGLIN.nocomp.R",
"R/analysis/model.lmer.all/model.lmer.LOGLIN.AD.R",
"R/analysis/model.lmer.all/model.lmer.LOGLIN.simplecomp.R")
model.files.lmer.Tf.1 <- c("R/analysis/model.lmer.all/model.lmer.LOGLIN.E.Tf.R",
"R/analysis/model.lmer.all/model.lmer.LOGLIN.R.Tf.R",
"R/analysis/model.lmer.all/model.lmer.LOGLIN.ER.Tf.R")
model.files.lmer.Tf.2 <- c("R/analysis/model.lmer.all/model.lmer.LOGLIN.nocomp.Tf.R",
"R/analysis/model.lmer.all/model.lmer.LOGLIN.AD.Tf.R",
"R/analysis/model.lmer.all/model.lmer.LOGLIN.HD.Tf.R",
"R/analysis/model.lmer.all/model.lmer.LOGLIN.simplecomp.Tf.R")
fun.call.lmer <- function(formula,df.lmer){
lmer.output <- lmer(formula=formula,data=df.lmer,REML = FALSE)
return(lmer.output)
}
fun.call.lmer.and.save <- function(formula,df.lmer,path.out){
lmer.output <- lmer(formula=formula,data=df.lmer,REML = FALSE)
print(summary(lmer.output))
saveRDS(lmer.output,file=file.path(path.out, "results.nolog.all.rds"))
}
run.lmer <- function (model.file, trait,
min.obs=10, sample.size=NA,
type.filling) {
require(lme4)
source(model.file, local = TRUE)
model <- load.model()
#= Path for output
path.out <- output.dir.lmer(model$name, trait, 'all',
type.filling=type.filling)
print(path.out)
dir.create(path.out, recursive=TRUE, showWarnings=FALSE)
cat("run lmer for model",model.file," for trait",
trait,"\n")
df.lmer <- load.and.prepare.data.for.lmer(trait,
min.obs, sample.size,
type.filling=type.filling) # return a DF
cat("Ok data with Nobs",nrow(df.lmer),
"\n")
#= Run model
fun.call.lmer.and.save(formula=model$lmer.formula.tree.id,df.lmer,path.out,std=std)
}
#========================================================================
output.dir.lmer <- function (model, trait, set,type.filling) {
file.path("output/lmer", set,trait,type.filling,model)
}
#============================================================
# Function to prepare data for lmer
#============================================================
load.and.prepare.data.for.lmer <- function(trait,
min.obs, sample.size, type.filling,
base.dir = "output/processed/"){
### load data
data.tree.tot <- read.csv(file.path(base.dir,"data.all.csv"),
stringsAsFactors = FALSE)
fun.data.for.lmer(data.tree.tot,trait,type.filling=type.filling)
}
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){
## remove tree with NA
data.tree <- subset(data.tree,subset=(!is.na(data.tree[["BA.G"]])) &
(!is.na(data.tree[["D"]])) )
## 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)
## select species with no missing traits
data.tree <- data.tree[(!is.na(data.tree[[abs.CWM.tntf]]) &
!is.na(data.tree[[BATOT]])),]
# select species with minimum obs
data.tree <- subset(data.tree,subset=data.tree[["sp"]] %in%
names(table(factor(data.tree[["sp"]])))[table(factor(data.tree[["sp"]]))>min.obs])
# select obs abov min perc.neigh
data.tree <- subset(data.tree,subset=data.tree[[perc.neigh]] > min.perc.neigh & !is.na(data.tree[[perc.neigh]]))
return(data.tree)
}
fun.center.and.standardized.var <- function(x){
return((x-mean(x))/sd(x))
}
### TODO NEED TO CONTIMNUE TO CHANGE FROM HERE
### 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){
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"]]))
tree.id <- unclass(factor(data.tree[["tree.id"]]))
plot.species.id <- unclass(factor(paste(data.tree[["plot"]],data.tree[["sp"]],sep="")))
#= multiply CWMs by BATOT
sumTnTfBn.abs <- data.tree[[abs.CWM.tntf]]*data.tree[[BATOT]]
sumTnBn <- data.tree[[CWM.tn]]*data.tree[[BATOT]]
sumTfBn <- data.tree[[tf]]*data.tree[[BATOT]]
sumTnTfBn.diff <- sumTnBn-sumTfBn
sumBn <- data.tree[[BATOT]]
return(data.frame(logG=logG,
logD=logD,
species.id=species.id,
tree.id=tree.id,
plot.species.id=plot.species.id,
sumTnTfBn.diff=sumTnTfBn.diff,
sumTnTfBn.abs=sumTnTfBn.abs,
Tf=data.tree[[tf]],
sumTnBn=sumTnBn,
sumTfBn=sumTfBn,
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){
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"]]))
tree.id <- unclass(factor(data.tree[["tree.id"]]))
plot.species.id <- unclass(factor(paste(data.tree[["plot"]],data.tree[["sp"]],sep="")))
#= multiply CWMs by BATOT
sumTnTfBn.abs <- data.tree[[abs.CWM.tntf]]*data.tree[[BATOT]]
sumTnBn <- data.tree[[CWM.tn]]*data.tree[[BATOT]]
sumTfBn <- data.tree[[tf]]*data.tree[[BATOT]]
sumTnTfBn.diff <- sumTnBn-sumTfBn
sumBn <- data.tree[[BATOT]]
return(data.frame(logG=logG,
logD=logD,
species.id=species.id,
plot.species.id=plot.species.id,
sumTnTfBn.diff=sumTnTfBn.diff,
sumTnTfBn.abs=sumTnTfBn.abs,
Tf=data.tree[[tf]],
sumTnBn=sumTnBn,
sumTfBn=sumTfBn,
sumBn=sumBn))
}
#============================================================
# Function to prepare data for lmer with new CWM data
# that will be used in a simple linear model
#============================================================
fun.data.for.lmer <- function(data.tree,trait,min.obs=10,type.filling='species') {
if(! trait %in% c("SLA", "Leaf.N","Seed.mass","SLA","Wood.density","Max.height")) stop("need trait to be in SLA Leaf.N Seed.mass SLA Wood.density or Max.height ")
# get vars names
CWM.tn <- paste(trait,"CWM",'fill',sep=".")
abs.CWM.tntf <- paste(trait,"abs.CWM",'fill',sep=".")
tf <- paste(trait,"focal",sep=".")
BATOT <- "BATOT"
perc.neigh <- paste(trait,"perc",type.filling,sep=".")
data.tree <- fun.select.data.for.analysis(data.tree,abs.CWM.tntf,perc.neigh,BATOT,min.obs)
#= DATA LIST FOR JAGS
if (length(table(table(data.tree[["tree.id"]])))>1){
lmer.data <- fun.get.the.variables.for.lmer.tree.id(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf)
}
if (length(table(table(data.tree[["tree.id"]])))<2){
lmer.data <- fun.get.the.variables.for.lmer.no.tree.id(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf)
}
return(lmer.data)
}
R/find.trait/France.R
View file @ db81af37
......@@ -2,7 +2,7 @@
##### FORMAT TRAIT FOR France
source("R/find.trait/trait-fun.R")
source("R/format.data/format.fun.R")
source("R/format.data/format-fun.R")
### read species names
......
R/find.trait/Luquillo.R
View file @ db81af37
......@@ -62,6 +62,8 @@ data.trait$Max.height.sd <- NA
data.trait$h <- NULL
data.TRAITS.std <- data.trait
write.csv(data.TRAITS.std,file="output/formatted/Luquillo/traits.std.csv",row.names = FALSE)
rm(data.trait)
## extract
data.traits <- fun.extract.format.sp.traits.NOT.TRY(sp=species.clean$Latin_name, Latin_name=species.clean$Latin_name, data=data.TRAITS.std,name.match.traits="Latin_name")
......
R/find.trait/Mbaiki.R
View file @ db81af37
......@@ -9,9 +9,12 @@
#
#===================================================#
# The following files must be in the Mbaiki raw data folder (../../data/raw/Mbaiki/)
# 1. "MBaiki_1995_2000_2005.csv" --> Inventory file (year 1995, 2000 and 2005) for Mbaiki with species list
# 2. "liaison_spmbaiki_bois-niv-sp.csv" --> Extract of the Cirad wood density data-base for Mbaiki species
# The following files must be in the Mbaiki raw data folder
# (../../data/raw/Mbaiki/)
# 1. "MBaiki_1995_2000_2005.csv" --> Inventory file (year 1995, 2000 and 2005)
# for Mbaiki with species list
# 2. "liaison_spmbaiki_bois-niv-sp.csv" --> Extract of the Cirad wood density
# data-base for Mbaiki species
# 3. "DensiteBoisSimpleMbaiki.csv" --> An additional file for wood density data
# 4. "Autour-de-Mbaiki-Releves-par-trait-et-taxon.txt" --> PlanTrait data
......@@ -33,10 +36,11 @@ nspecies <- length(List.species)
# Building the trait data-set
data.species <- data.frame(Species=List.species,
Leaf.N.mean=NA,Leaf.N.sd=NA,Leaf.N.n=NA,
Seed.mass.mean=NA,Seed.mass.sd=NA,Seed.mass.n=NA,
SLA.mean=NA,SLA.sd=NA,SLA.n=NA,
Wood.density.mean=NA,Wood.density.sd=NA,Wood.density.n=NA)
Leaf.N.mean=NA, Leaf.N.sd=NA, Leaf.N.n=NA,
Seed.mass.mean=NA, Seed.mass.sd=NA, Seed.mass.n=NA,
SLA.mean=NA, SLA.sd=NA, SLA.n=NA,
Wood.density.mean=NA, Wood.density.sd=NA,
Wood.density.n=NA)
#= Wood density
......@@ -44,24 +48,33 @@ data.species <- data.frame(Species=List.species,
data.wd <- read.csv2("data/raw/Mbaiki/liaison_spmbaiki_bois-niv-sp.csv")
index <- which(duplicated(data.wd$Scientifique.Name.1995))
data.wd.2 <- data.wd[-index,] # remove duplicates
# Load another Cirad wood density data-base (from Cindy Gidoin and Sylvie Gourlet-Fleury)
data.wd.Cindy <- read.table("data/raw/Mbaiki/DensiteBoisSimpleMbaiki.csv",header=TRUE,sep="\t")
# Load another Cirad wood density data-base
# (from Cindy Gidoin and Sylvie Gourlet-Fleury)
data.wd.Cindy <- read.table("data/raw/Mbaiki/DensiteBoisSimpleMbaiki.csv",
header=TRUE,sep="\t")
data.wd.Cindy.2 <- data.wd.Cindy[data.wd.Cindy$NivMoyID=="E",]
data.wd.Cindy.2$LatinName <- paste(data.wd.Cindy.2$Genre,data.wd.Cindy.2$Espece,sep=" ")
data.wd.Cindy.2$LatinName <- paste(data.wd.Cindy.2$Genre,
data.wd.Cindy.2$Espece,sep=" ")
ListSpeciesCindy <- levels(as.factor(data.wd.Cindy.2$LatinName))
# Merge
data.species.2 <- merge(data.species,data.wd.2,by.x="Species",by.y="Scientifique.Name.1995",all.x=TRUE)
data.species.3 <- merge(data.species.2,data.wd.Cindy.2,by.x="Species",by.y="LatinName",all.x=TRUE)
data.species.2 <- merge(data.species, data.wd.2,by.x="Species",
by.y="Scientifique.Name.1995",all.x=TRUE)
data.species.3 <- merge(data.species.2,data.wd.Cindy.2,by.x="Species",
by.y="LatinName",all.x=TRUE)
# Fill-up wood density variables
data.species.3$Wood.density.mean <- as.numeric(as.character(data.species.3$moy_infra_densite))
data.species.3$Wood.density.sd <- sqrt(as.numeric(as.character(data.species.3$var_infra_densite)))
data.species.3$Wood.density.n <- as.numeric(as.character(data.species.3$nb_arbre))
data.species.3$Wood.density.mean <- as.numeric(as.character(
data.species.3$moy_infra_densite))
data.species.3$Wood.density.sd <- sqrt(as.numeric(as.character(
data.species.3$var_infra_densite)))
data.species.3$Wood.density.n <- as.numeric(as.character(
data.species.3$nb_arbre))
sum(!is.na(data.species.3$Wood.density.mean))
# Complete with Cindy data
for (i in 1:nrow(data.species.3)) {
sp.i <- data.species.3$Species[i]
b <- ifelse(sp.i %in% ListSpeciesCindy, data.wd.Cindy.2$MoyID[data.wd.Cindy.2$LatinName==sp.i], NA)
b <- ifelse(sp.i %in% ListSpeciesCindy,
data.wd.Cindy.2$MoyID[data.wd.Cindy.2$LatinName==sp.i], NA)
c <- data.wd.Cindy.2$VarID[data.wd.Cindy.2$LatinName==sp.i]
d <- data.wd.Cindy.2$NivVarID[data.wd.Cindy.2$LatinName==sp.i]
if (is.na(data.species.3$Wood.density.mean[i]) & !is.na(b)) {
......@@ -152,11 +165,29 @@ species.clean <- data.frame(sp=data.tree[!duplicated(data.tree[["sp"]]),"sp"],
## extract
data.traits <- fun.extract.format.sp.traits.NOT.TRY(sp=species.clean$Latin_name, Latin_name=species.clean$Latin_name, data=data.TRAITS.std,name.match.traits="Species")
data.traits <- fun.extract.format.sp.traits.NOT.TRY(sp=species.clean$Latin_name,
Latin_name=species.clean$Latin_name,
data=data.TRAITS.std,
name.match.traits="Species")
## change sp
data.traits$sp <- species.clean$sp
## EXTRACTION FROM TRY too fill missing value
## read in data
data.TRY.std <- readRDS("output/formatted/TRY/data.TRY.std.rds")
data.traits2 <- fun.extract.format.sp.traits.TRY(sp=species.clean[["Latin_name"]],
sp.syno.table=species.clean,
data=data.TRY.std)
## TODO NEED TO COMBINE BOTH
data.traits <- fun.combine.nontry.and.try(trait = c('Leaf.N', 'SLA', 'Seed.mass',
'Wood.density'),
data1 = data.traits,
data2 = data.traits2)
#### GET THE ANGIO/CONIF AND EVERGREEN/DECIDUOUS
# read try categrocial data
......@@ -172,29 +203,15 @@ 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.traits <- merge(data.traits,data.cat.extract[,c("sp","Phylo.group","Pheno.T")],by="sp")
data.traits <- merge(data.traits,
data.cat.extract[,c("sp","Phylo.group","Pheno.T")],
by="sp")
## compute perc of traits cover per species
print(sapply(c('Leaf.N','SLA','Seed.mass','Wood.density','Max.height'),
fun.compute.perc.cover.one.trait,data.traits))
### Export
write.csv(data.traits,file="output/formatted/Mbaiki/traits.csv",row.names = FALSE)
# Export
#data.species.4$sp <- data.species.4$Species
#write.csv(data.species.4,file="output/formatted/Mbaiki/traits.csv",row.names=FALSE)
## #===================================
## # How many species have trait values
## nspecies
## # Obs at the species level
## nsp.LeafN <- sum(!is.na(data.species.4$Leaf.N.mean))
## nsp.SeedMass <- sum(!is.na(data.species.4$Seed.mass.mean))
## nsp.SLA <- sum(!is.na(data.species.4$SLA.mean))
## nsp.WSG <- sum(!is.na(data.species.4$Wood.density.mean))
## nsp.AllTraits <- sum(!is.na(data.species.4$Leaf.N.mean) & !is.na(data.species.4$Seed.mass.mean) &
## !is.na(data.species.4$SLA.mean) & !is.na(data.species.4$Wood.density.mean))
## # Summary in a matrix
## matsum <- as.data.frame(matrix(nrow=1,ncol=6))
## names(matsum) <- c("Total","LeafN","SeedMass","SLA","WSG","AllTraits")
## matsum[1,] <- c(nspecies,nsp.LeafN,nsp.SeedMass,nsp.SLA,nsp.WSG,nsp.AllTraits)
## sink("Summary_Traits_Mbaiki.txt")
## matsum
## sink()
write.csv(data.traits,file="output/formatted/Mbaiki/traits.csv",
row.names = FALSE)
R/find.trait/Paracou.R
View file @ db81af37
This diff is collapsed.
R/find.trait/US.R
View file @ db81af37
......@@ -13,16 +13,22 @@ species.clean <- data.frame(sp=data.tree[!duplicated(data.tree[["sp"]]),"sp"],
## read in data
data.TRY.std <- readRDS("output/formatted/TRY/data.TRY.std.rds")
max.height <- read.csv(file="output/formatted/US/max.height.csv", stringsAsFactors = FALSE)
max.height <- read.csv(file="output/formatted/US/max.height.csv",
stringsAsFactors = FALSE)
max.height$sp <- paste("sp",max.height$sp,sep=".")
data.traits <- fun.extract.format.sp.traits.TRY(sp=species.clean[["sp"]],sp.syno.table=species.clean,data=data.TRY.std)
data.traits <- fun.extract.format.sp.traits.TRY(sp=species.clean[["sp"]],
sp.syno.table=species.clean,
data=data.TRY.std)
data.traits <- merge(data.traits,subset(max.height,select=c("sp","Max.height.mean","Max.height.sd")),by="sp",all.x=TRUE,all.y=FALSE)
data.traits <- merge(data.traits,subset(max.height,
select=c("sp","Max.height.mean","Max.height.sd")),
by="sp",all.x=TRUE,all.y=FALSE)
data.traits$Max.height.genus <- FALSE
# genus mean for height
height.genus.DF <- do.call("rbind",lapply(data.traits$Latin_name,fun.compute.mean.genus,data.traits,"Max.height.mean"))
height.genus.DF <- do.call("rbind",lapply(data.traits$Latin_name,fun.compute.mean.genus,
data.traits,"Max.height.mean"))
data.traits[is.na(data.traits[["Max.height.mean"]]),
c("Max.height.mean","Max.height.sd","Max.height.genus")] <- height.genus.DF[is.na(data.traits[["Max.height.mean"]]),]
......
R/find.trait/trait-fun.R
View file @ db81af37
......@@ -354,3 +354,25 @@ 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']
return(data.cat.extract)
}
## compute perc of traits cover per species
fun.compute.perc.cover.one.trait <- function(trait, data) {
t.mean <- paste(trait,"mean",sep=".")
t.genus <- paste(trait,"genus",sep=".")
sp.perc <- sum(!is.na(data[[t.mean]]) & !data[[t.genus]])/nrow(data)
genus.perc <- sum(!is.na(data[[t.mean]]))/nrow(data)
return( c(sp.perc=sp.perc,genus.perc=genus.perc) )
}
#### function to combine nontry with try data
fun.combine.nontry.and.try <- function(trait, data1, data2) {
t.mean <- paste(trait,"mean",sep=".")
t.genus <- paste(trait,"genus",sep=".")
for (i in 1:length(t.mean)) {
data1[is.na(data1[[t.mean[i]]]),
c(t.mean[i],t.genus[i])] <- data2[is.na(data1[[t.mean[i]]]),
c(t.mean[i],t.genus[i])]
}
return( data1 )
}
R/format.data/France.R
View file @ db81af37
......@@ -77,7 +77,15 @@ data.france <- subset(data.france, subset = data.france[["plisi"]] == 0) ## no
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
data.france$MAP <- data.france$SAP
### get wc climate
data.france$MAT <- NULL
source("R/utils/climate.R")
clim <- GetClimate(data.france$Lat,data.france$Lon)
data.france$MAT <- clim$MAT
data.france$MAP <- clim$MAP
## SELECT GOOD COLUMNS
## names of variables for abiotic conditions
vec.abio.var.names <- c("MAT", "MAP")
......
R/format.data/Luquillo.R
View file @ db81af37
......@@ -116,10 +116,10 @@ data.luq[["ecocode"]] <- "tropical"
###################### 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.luq, select = c(vec.basic.var))
y
data.tree <- subset(data.tree,subset=!is.na(data.tree$x) & !is.na(data.tree$y))
## convert var factor in character or numeric
data.tree <- fun.convert.type.B(data.tree)
......
R/format.data/Mbaiki.R
View file @ db81af37
......@@ -24,17 +24,31 @@ data.mbaiki$tree.id <- apply(data.mbaiki[,c("cluster","plot","tree")],1,paste,co
# At Mbaiki, a cluster is made of four plots of 100m x 100m arrange as 200m x 200m forest clusters
# The coordinates of the trees are given for each cluster (origin at the bottom of the left side).
Levels.cluster <- levels(as.factor(data.mbaiki$cluster))
n.cluster <- length(Levels.cluster) # 10 clusters => 10 "plots"
for (i in Levels.cluster){
dat <- subset(data.mbaiki, subset=data.mbaiki[["cluster"]]==i)
x11()
symbols(dat$x,dat$y,circles=dat$circum2000,inches=0.2,main=i)
}
## Ghislain's comment: OK, tree location in clusters looks good
## create small plot
square.s.t <- 20
make.quad <- do.call("rbind",lapply(unique(data.mbaiki$cluster),
FUN=fun.quadrat.cluster,
cluster=data.mbaiki$cluster,
tree.id=data.mbaiki$tree.id,
x=data.mbaiki$x,
y=data.mbaiki$y,
square.s=square.s.t))
data.mbaiki <- merge(data.mbaiki, make.quad,by='tree.id')
## Levels.cluster <- levels(as.factor(data.mbaiki$cluster))
## n.cluster <- length(Levels.cluster) # 10 clusters => 10 "plots"
## for (i in Levels.cluster){
## dat <- subset(data.mbaiki, subset=data.mbaiki[["cluster"]]==i)
## x11()
## symbols(dat$x,dat$y,circles=dat$circum2000,inches=0.2,main=i,
## fg=unclass(factor(dat$make.quad)))
## }
### read species names
species.clean <- read.csv("data/raw/Mbaiki/Mbaiki.traits.csv",stringsAsFactors=FALSE, header = T, sep = ",")
species.clean$sp.name <- species.clean$Species
......@@ -107,7 +121,7 @@ data.mbaiki$D <- data.mbaiki[["dbh1"]]; data.mbaiki$D[data.mbaiki$D == 0] <- NA
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)
data.mbaiki$plot <- paste(data.mbaiki$cluster,data.mbaiki$plot)
data.mbaiki$plot <- data.mbaiki$make.quad
### delete recruit in 2000 or 2005 for first census
data.mbaiki <- subset(data.mbaiki,subset=!is.na(data.mbaiki$D))
......
R/format.data/NSW.R
View file @ db81af37
......@@ -60,65 +60,13 @@ data.nswtnd$plot <- rep(NA,nrow(data.nswtnd))
data.nsw <- rbind(data.nswbrc, data.nswbrt, data.nswbs1, data.nswbs2, data.nswtnd)
data.nsw$Plot <- NULL
##################### CREATE PLOT BASED ON square.sxsquare.s m cell from X Y
fun.quadrat <- function(x,y,square.s) {
if(sum(!is.na(x))>10){
vec.x <- seq(0,max(x,na.rm=T)+9,by=square.s)
vec.y <- seq(0,max(y,na.rm=T)+9,by=square.s)
x.cut <- cut(x,breaks=vec.x)
y.cut <- cut(y,breaks=vec.y)
out <- apply(cbind(x.cut,y.cut),1,paste,collapse=".")
out[is.na(x.cut)] <- NA
out[is.na(y.cut)] <- NA
return(unclass(out))
}else{
return(rep(NA,length(x)))
}
}
## function to apply per cluster
fun.quadrat.cluster <- function(cluster.id,cluster,x,y,square.s){
temp <- fun.quadrat(x[cluster==cluster.id], y[cluster==cluster.id],square.s=square.s)
temp[!is.na(temp)] <- paste((cluster[cluster==cluster.id])[!is.na(temp)],temp[!is.na(temp)],sep=".")
return(as.vector(temp))
}
square.s.t <- 20
make.quad <- unlist(sapply(unique(data.nsw$cluster),FUN=fun.quadrat.cluster,cluster=data.nsw$cluster,x=data.nsw$x, y=data.nsw$y,square.s=square.s.t))
data.nsw <- cbind(data.nsw, make.quad)
## plot to check problem
nn <- 0
for (i in unique(data.nsw$cluster)){
if (sum(!is.na(data.nsw$x[data.nsw$cluster == i]))>3){
x11()
nn <- nn+1
cat("\n\n\n",nn,"\n")
print(i)
print(range(data.nsw$x[data.nsw$cluster == i],na.rm=TRUE))
print(table(factor(data.nsw$make.quad[data.nsw$cluster == i])))
plot(data.nsw$x[data.nsw$cluster == i],data.nsw$y[data.nsw$cluster == i],col=unclass(factor(data.nsw$make.quad[data.nsw$cluster == i])),main=i)
points(data.nsw$x[data.nsw$cluster == i],data.nsw$y[data.nsw$cluster == i],cex=0.2)
vec.x <- seq(0,max(data.nsw$x[data.nsw$cluster == i],na.rm=T)+9,by=square.s.t)
vec.y <- seq(0,max(data.nsw$y[data.nsw$cluster == i],na.rm=T)+9,by=square.s.t)
abline(v=vec.x)
abline(h=vec.y)
}
}
## subplot becomes plot; plot becomes combination of x and y
data.nsw$plot[is.na(data.nsw$plot)] <- as.vector(unclass(data.nsw$make.quad))[is.na(data.nsw$plot)] ## plot code
data.nsw$make.quad <- NULL
########################################## FORMAT INDIVIDUAL TREE DATA - Each tree has at most 3 observations
data.nsw$tree.id <- apply(data.nsw[, c("cluster","Tree.No")], 1, paste, collapse = "_")
table(table(data.nsw$tree.id))
data.nsw2 <- data.frame(data.nsw[1, ], ye