Skip to content

GitLab

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

fixe problem with plot size in the US

parent 8baf852a
Hide whitespace changes
Inline Side-by-side
Makefile
View file @ f1d3e6cc
......@@ -27,7 +27,7 @@ $(D2)/TRY/data.TRY.std.rds:
BCI: $(D3)/BCI/Done.txt
$(D3)/BCI/Done.txt: R/process.data/process-fun.R $(D2)/BCI/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('BCI', Rlim=15,std.traits='local'); process_bigplot_dataset('BCI', Rlim=15,std.traits='global'); process_bigplot_dataset('BCI', Rlim=15,std.traits='no');"
Rscript -e "source('$<'); process_bigplot_dataset('BCI', Rlim=15,std.traits='local'); process_bigplot_dataset('BCI', Rlim=15,std.traits='global');"
$(D2)/BCI/traits.csv: R/find.trait/BCI.R R/find.trait/trait-fun.R $(D2)/BCI/tree.csv
Rscript $<
......@@ -40,7 +40,7 @@ $(D2)/BCI/tree.csv: R/format.data/BCI.R $(shell find $(D1)/BCI -type f)
Japan: $(D3)/Japan/Done.txt
$(D3)/Japan/Done.txt: R/process.data/process-fun.R $(D2)/Japan/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Japan', Rlim=15,std.traits='local'); process_bigplot_dataset('Japan', Rlim=15,std.traits='global'); process_bigplot_dataset('Japan', Rlim=15,std.traits='no');"
Rscript -e "source('$<'); process_bigplot_dataset('Japan', Rlim=15,std.traits='local'); process_bigplot_dataset('Japan', Rlim=15,std.traits='global');"
$(D2)/Japan/traits.csv: R/find.trait/Japan.R R/find.trait/trait-fun.R $(D2)/Japan/tree.csv
Rscript $<
......@@ -52,7 +52,7 @@ $(D2)/Japan/tree.csv: R/format.data/Japan.R $(shell find $(D1)/Japan -type f)
Luquillo: $(D3)/Luquillo/Done.txt
$(D3)/Luquillo/Done.txt: R/process.data/process-fun.R $(D2)/Luquillo/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Luquillo', Rlim=15,std.traits='local'); process_bigplot_dataset('Luquillo', Rlim=15,std.traits='no');process_bigplot_dataset('Luquillo', Rlim=15,std.traits='global');"
Rscript -e "source('$<'); process_bigplot_dataset('Luquillo', Rlim=15,std.traits='local'); process_bigplot_dataset('Luquillo', Rlim=15,std.traits='no');"
$(D2)/Luquillo/traits.csv: R/find.trait/Luquillo.R R/find.trait/trait-fun.R $(D2)/Luquillo/tree.csv
Rscript $<
......@@ -65,7 +65,7 @@ $(D2)/Luquillo/tree.csv: R/format.data/Luquillo.R $(shell find $(D1)/Luquillo -t
Mbaiki: $(D3)/Mbaiki/Done.txt
$(D3)/Mbaiki/Done.txt: R/process.data/process-fun.R $(D2)/Mbaiki/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Mbaiki', Rlim=15,std.traits='local'); process_bigplot_dataset('Mbaiki', Rlim=15,std.traits='global'); process_bigplot_dataset('Mbaiki', Rlim=15,std.traits='no');"
Rscript -e "source('$<'); process_bigplot_dataset('Mbaiki', Rlim=15,std.traits='local'); process_bigplot_dataset('Mbaiki', Rlim=15,std.traits='global');"
$(D2)/Mbaiki/traits.csv: R/find.trait/Mbaiki.R R/find.trait/trait-fun.R $(D2)/Mbaiki/tree.csv
Rscript $<
......@@ -78,7 +78,7 @@ $(D2)/Mbaiki/tree.csv: R/format.data/Mbaiki.R $(shell find $(D1)/Mbaiki -type f)
Canada: $(D3)/Canada/Done.txt
$(D3)/Canada/Done.txt: R/process.data/process-fun.R $(D2)/Canada/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('Canada',std.traits='local'); process_inventory_dataset('Canada',std.traits='global');process_inventory_dataset('Canada',std.traits='no');"
Rscript -e "source('$<'); process_inventory_dataset('Canada',std.traits='local'); process_inventory_dataset('Canada',std.traits='global');"
$(D2)/Canada/traits.csv: R/find.trait/Canada.R R/find.trait/trait-fun.R $(D2)/Canada/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $<
......@@ -90,7 +90,7 @@ $(D2)/Canada/tree.csv: R/format.data/Canada.R $(shell find $(D1)/Canada -type f)
France: $(D3)/France/Done.txt
$(D3)/France/Done.txt: R/process.data/process-fun.R $(D2)/France/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('France',std.traits='local'); process_inventory_dataset('France',std.traits='global');process_inventory_dataset('France',std.traits='no');"
Rscript -e "source('$<'); process_inventory_dataset('France',std.traits='local'); process_inventory_dataset('France',std.traits='global');"
$(D2)/France/traits.csv: R/find.trait/France.R R/find.trait/trait-fun.R $(D2)/France/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $<
......@@ -103,7 +103,7 @@ $(D2)/France/tree.csv: R/format.data/France.R $(shell find $(D1)/France -type f)
Fushan: $(D3)/Fushan/Done.txt
$(D3)/Fushan/Done.txt: R/process.data/process-fun.R $(D2)/Fushan/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Fushan', Rlim=15,std.traits='local'); process_bigplot_dataset('Fushan', Rlim=15,std.traits='global');process_bigplot_dataset('Fushan', Rlim=15,std.traits='no');"
Rscript -e "source('$<'); process_bigplot_dataset('Fushan', Rlim=15,std.traits='local'); process_bigplot_dataset('Fushan', Rlim=15,std.traits='global');"
$(D2)/Fushan/traits.csv: R/find.trait/Fushan.R R/find.trait/trait-fun.R $(D2)/Fushan/tree.csv
......@@ -117,7 +117,7 @@ $(D2)/Fushan/tree.csv: R/format.data/Fushan.R $(shell find $(D1)/Fushan -type f)
NSW: $(D3)/NSW/Done.txt
$(D3)/NSW/Done.txt: R/process.data/process-fun.R $(D2)/NSW/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('NSW',std.traits='local'); process_inventory_dataset('NSW',std.traits='global');process_inventory_dataset('NSW',std.traits='no');"
Rscript -e "source('$<'); process_inventory_dataset('NSW',std.traits='local'); process_inventory_dataset('NSW',std.traits='global');"
$(D2)/NSW/traits.csv: R/find.trait/NSW.R R/find.trait/trait-fun.R $(D2)/NSW/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $<
......@@ -130,7 +130,7 @@ $(D2)/NSW/tree.csv: R/format.data/NSW.R $(shell find $(D1)/NSW -type f)
NVS: $(D3)/NVS/Done.txt
$(D3)/NVS/Done.txt: R/process.data/process-fun.R $(D2)/NVS/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('NVS',std.traits='local'); process_inventory_dataset('NVS',std.traits='global'); process_inventory_dataset('NVS',std.traits='no');"
Rscript -e "source('$<'); process_inventory_dataset('NVS',std.traits='local'); process_inventory_dataset('NVS',std.traits='global');"
$(D2)/NVS/traits.csv: R/find.trait/NVS.R R/find.trait/trait-fun.R $(D2)/NVS/tree.csv
Rscript $<
......@@ -143,7 +143,7 @@ $(D2)/NVS/tree.csv: R/format.data/NVS.R $(shell find $(D1)/NVS -type f)
Paracou: $(D3)/Paracou/Done.txt
$(D3)/Paracou/Done.txt: R/process.data/process-fun.R $(D2)/Paracou/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Paracou', Rlim=15,std.traits='local'); process_bigplot_dataset('Paracou', Rlim=15,std.traits='global'); process_bigplot_dataset('Paracou', Rlim=15,std.traits='no');"
Rscript -e "source('$<'); process_bigplot_dataset('Paracou', Rlim=15,std.traits='local'); process_bigplot_dataset('Paracou', Rlim=15,std.traits='global');"
$(D2)/Paracou/traits.csv: R/find.trait/Paracou.R R/find.trait/trait-fun.R $(D2)/Paracou/tree.csv
Rscript $<
......@@ -156,7 +156,7 @@ $(D2)/Paracou/tree.csv: R/format.data/Paracou.R $(shell find $(D1)/Paracou -type
Spain: $(D3)/Spain/Done.txt
$(D3)/Spain/Done.txt: R/process.data/process-fun.R $(D2)/Spain/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('Spain',std.traits='local'); process_inventory_dataset('Spain',std.traits='global'); process_inventory_dataset('Spain',std.traits='no');"
Rscript -e "source('$<'); process_inventory_dataset('Spain',std.traits='local'); process_inventory_dataset('Spain',std.traits='global');"
$(D2)/Spain/traits.csv: R/find.trait/Spain.R R/find.trait/trait-fun.R $(D2)/Spain/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $<
......@@ -169,7 +169,7 @@ $(D2)/Spain/tree.csv: R/format.data/Spain.R $(shell find $(D1)/Spain -type f)
Sweden: $(D3)/Sweden/Done.txt
$(D3)/Sweden/Done.txt: R/process.data/process-fun.R $(D2)/Sweden/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('Sweden',std.traits='local'); process_inventory_dataset('Sweden',std.traits='global'); process_inventory_dataset('Sweden',std.traits='no');"
Rscript -e "source('$<'); process_inventory_dataset('Sweden',std.traits='local'); process_inventory_dataset('Sweden',std.traits='global');"
$(D2)/Sweden/traits.csv: R/find.trait/Sweden.R R/find.trait/trait-fun.R $(D2)/Sweden/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $<
......@@ -182,7 +182,7 @@ $(D2)/Sweden/tree.csv: R/format.data/Sweden.R $(shell find $(D1)/Sweden -type f)
Swiss: $(D3)/Swiss/Done.txt
$(D3)/Swiss/Done.txt: R/process.data/process-fun.R $(D2)/Swiss/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('Swiss',std.traits='local'); process_inventory_dataset('Swiss',std.traits='global'); process_inventory_dataset('Swiss',std.traits='no');"
Rscript -e "source('$<'); process_inventory_dataset('Swiss',std.traits='local'); process_inventory_dataset('Swiss',std.traits='global');"
$(D2)/Swiss/traits.csv: R/find.trait/Swiss.R R/find.trait/trait-fun.R $(D2)/Swiss/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $<
......@@ -195,7 +195,7 @@ $(D2)/Swiss/tree.csv: R/format.data/Swiss.R $(shell find $(D1)/Swiss -type f)
US: $(D3)/US/Done.txt
$(D3)/US/Done.txt: R/process.data/process-fun.R $(D2)/US/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('US',std.traits='local'); process_inventory_dataset('US',std.traits='global'); process_inventory_dataset('US',std.traits='no');"
Rscript -e "source('$<'); process_inventory_dataset('US',std.traits='local'); process_inventory_dataset('US',std.traits='global');"
$(D2)/US/traits.csv: R/find.trait/US.R R/find.trait/trait-fun.R $(D2)/US/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $<
......
R/analysis/lmer.nolog.all.output.R 0 → 100644
View file @ f1d3e6cc
#!/usr/bin/env Rscript
source("R/analysis/lmer.output-fun.R")
source("R/analysis/lmer.run.nolog.all.R")
## LOOP OVER FILES AND NOT SCAN ALL FILES
traits <- c("SLA", "Wood.density","Max.height","Leaf.N","Seed.mass")
type.filling <- 'species'
files <- c()
for (trait in trait){
for(model in c(model.files.lmer.Tf.1,model.files.lmer.Tf.2)){
source(model, local = TRUE)
model.obj <- load.model()
pathout <- output.dir.lmer(model.obj$name, trait, 'all',
type.filling=type.filling)
files <- c(files,file.path(pathout,"results.nolog.all.rds"))
}
}
out <- lapply(files, summarise.lmer.output.all.list)
names(out) <- lapply(lapply(files,files.details.all),
function(x) paste(as.vector(x[names(x) != 'file']),
collapse="_"))
### remove missing
out <- out[!unlist(lapply(out,FUN=function(x) is.null(x$lmer.summary)))]
saveRDS(out,file='output/list.lmer.out.nolog.all.rds')
names.param <- unique(unlist(lapply(out,function(list.res)
names(list.res$lmer.summary$fixed.coeff.E))))
DF.results <- do.call("rbind",lapply(out, fun.format.in.data.frame,
names.param = names.param))
DF.results$id <- paste(DF.results$set, DF.results$ecocode, sep=".")
saveRDS(DF.results,file='output/DF.lmer.nolog.all.rds')
R/analysis/lmer.nolog.output.R
View file @ f1d3e6cc
......@@ -3,7 +3,7 @@
source("R/analysis/lmer.output-fun.R")
source("R/analysis/lmer.run.nolog.R")
## TODO NEED TO CHANGE THAT TO LOOP OVER FILES AND NOT SCAN ALL FILES
## LOOP OVER FILES AND NOT SCAN ALL FILES
sets <- c('BCI','Canada','France','Fushan','NVS','Paracou',
'Spain','US','Swiss','Sweden','NSW','Mbaiki','Luquillo','Japan')
traits <- c("SLA", "Wood.density","Max.height","Leaf.N","Seed.mass")
......@@ -17,7 +17,8 @@ ecoregions <- get.ecoregions.for.set(set)
source(model, local = TRUE)
model.obj <- load.model()
pathout <- output.dir.lmer(model=model.obj$name, trait, set, ecoregion,type.filling)
pathout <- output.dir.lmer(model=model.obj$name, trait,
set, ecoregion,type.filling)
files <- c(files,file.path(pathout,"results.nolog.rds"))
}
}
......@@ -26,11 +27,16 @@ ecoregions <- get.ecoregions.for.set(set)
out <- lapply(files, summarise.lmer.output.list)
names(out) <- lapply(lapply(files,files.details),function(x) paste(as.vector(x[-6]),collapse="_"))
names(out) <- lapply(lapply(files,files.details),
function(x) paste(as.vector(x[names(x) != 'file']),
collapse="_"))
### remove missing
out <- out[!unlist(lapply(out,FUN=function(x) is.null(x$lmer.summary)))]
saveRDS(out,file='output/lmer.list.out.nolog.rds')
names.param <- unique(unlist(lapply(out,function(list.res) names(list.res$lmer.summary$fixed.coeff.E))))
DF.results <- do.call("rbind",lapply(out,fun.format.in.data.frame,names.param=names.param))
DF.results$id <- paste(DF.results$set,DF.results$ecocode,sep=".")
names.param <- unique(unlist(lapply(out,
function(list.res) names(list.res$lmer.summary$fixed.coeff.E))))
DF.results <- do.call("rbind",lapply(out, fun.format.in.data.frame,
names.param = names.param))
DF.results$id <- paste(DF.results$set, DF.results$ecocode, sep=".")
saveRDS(DF.results,file='output/lmer.nolog.DF.rds')
R/analysis/lmer.output-fun.R
View file @ f1d3e6cc
......@@ -36,6 +36,18 @@ summarise.lmer.output.list <- function(f ){
}
summarise.lmer.output.all.list <- function(f ){
output.lmer <- read.lmer.output(f)
if(!is.null(output.lmer)){
res <- list(files.details=files.details.all(f),
lmer.summary=summarise.lmer.output( output.lmer))
}else{
res <- NULL
}
return(res)
}
files.details <- function(x){
s <- data.frame(t(strsplit(x, "/", fixed=TRUE)[[1]]), stringsAsFactors= FALSE)
......@@ -44,8 +56,14 @@ files.details <- function(x){
}
files.details.all <- function(x){
s <- data.frame(t(strsplit(x, "/", fixed=TRUE)[[1]]), stringsAsFactors= FALSE)
names(s) <- c("d1", "d2", "set", "trait", "filling", "model", "file" )
s[-(1:2)]
}
#### R squred functions
#### R squared functions
# Function rsquared.glmm requires models to be input as a list (can include fixed-
# effects only models,but not a good idea to mix models of class "mer" with models
......
R/format.data/US.R
View file @ f1d3e6cc
......@@ -7,26 +7,34 @@ dir.create("output/formatted/US", recursive=TRUE,showWarnings=FALSE)
### read species names
species.clean <- read.csv("data/raw/US/REF_SPECIES.CSV", stringsAsFactors = FALSE)
species.clean <- read.csv("data/raw/US/REF_SPECIES.CSV",
stringsAsFactors = FALSE)
## select column to keep
species.clean <- subset(species.clean, select = c("SPCD", "GENUS", "SPECIES", "VARIETY",
"SUBSPECIES", "SPECIES_SYMBOL"))
species.clean$Latin_name <- paste(species.clean[["GENUS"]], species.clean[["SPECIES"]],
sep = " ")
species.clean$Latin_name_syn <- paste(species.clean[["GENUS"]], species.clean[["SPECIES"]],
species.clean <- subset(species.clean,
select = c("SPCD", "GENUS", "SPECIES", "VARIETY",
"SUBSPECIES", "SPECIES_SYMBOL"))
species.clean$Latin_name <- paste(species.clean[["GENUS"]],
species.clean[["SPECIES"]],
sep = " ")
species.clean$Latin_name_syn <- paste(species.clean[["GENUS"]],
species.clean[["SPECIES"]],
sep = " ")
names(species.clean)[1] <- "sp"
species.clean[["sp"]] <- paste("sp", species.clean[["sp"]], sep = ".")
## LOAD QUANTILE OF HEIGHT AND FORMAT AS OTHER
data.max.height <- read.csv("data/raw/US/FiaSpMaxHt.csv", stringsAsFactors = FALSE)
data.max.height <- data.frame(sp=data.max.height$SpCd,Max.height.mean=data.max.height$Ht99,
Max.height.sd=rep(NA,nrow(data.max.height)), Max.height.nobs=rep(NA,nrow(data.max.height)))
write.csv(data.max.height, file = "output/formatted/US/max.height.csv")
data.max.height <- read.csv("data/raw/US/FiaSpMaxHt.csv",
stringsAsFactors = FALSE)
data.max.height <- data.frame(sp=data.max.height$SpCd,
Max.height.mean=data.max.height$Ht99,
Max.height.sd=rep(NA,nrow(data.max.height)),
Max.height.nobs=rep(NA,nrow(data.max.height)))
write.csv(data.max.height, file = "output/formatted/US/max.height.csv")
#### READ DATA US
data.us <- read.csv("data/raw/US/FIA51_trees_w_supp.csv",header=TRUE,stringsAsFactors =FALSE)
data.us <- read.csv("data/raw/US/FIA51_trees_w_supp.csv",
header = TRUE, stringsAsFactors = FALSE)
#MCV remove trees <10 cm dbh
data.us <- subset(data.us,subset=(data.us$InitDbh>=10) )
......@@ -37,63 +45,80 @@ data.us <- subset(data.us,subset=(data.us$InitDbh>=10) )
############## WHEN WE ANALYZE THAT DATASET LATER ON
############## FORMAT INDIVIDUAL TREE DATA
## change unit and names of variables to be the same in all data for the tree
data.us$G <- 10 * (data.us$FinalDbh - data.us$InitDbh)/data.us$Interval ## diameter growth in mm per year
data.us$BA.G <- (pi*(data.us$FinalDbh/2)^2-pi*(data.us$InitDbh/2)^2)/data.us$Interval ## ba growth in cm^2/year
data.us$G <- 10 * (data.us$FinalDbh - data.us$InitDbh)/data.us$Interval
## diameter growth in mm per year
data.us$BA.G <- (pi*(data.us$FinalDbh/2)^2-pi*(data.us$InitDbh/2)^2)/
data.us$Interval ## ba growth in cm^2/year
data.us$G[which(data.us$InitDbh == 0 | data.us$FinalDbh == -999)] <- NA
data.us$BA.G[which(data.us$InitDbh == 0 | data.us$FinalDbh == -999)] <- NA
## ## remove too big error NOT AT THIS STAGE LATER ON IN SAME WAY FOR ALL DATA
## data.us$G[which((data.us$G < -15 | data.us$G>100) & !is.na(data.us$G) )] <- NA
## data.us$BA.G[which((data.us$G < -15 | data.us$G>100) & !is.na(data.us$G) )] <- NA
data.us$year <- data.us$Interval ## number of year between measuremen
data.us$D <- data.us[["InitDbh"]]
data.us$D[data.us$D == 0] <- NA
## diameter in cm
data.us$dead <- as.numeric(data.us$FinalDbh > 0) ## dummy variable for dead tree 0 alive 1 dead
data.us$dead <- as.numeric(data.us$FinalDbh > 0)
## dummy variable for dead tree 0 alive 1 dead
data.us$sp <- as.character(data.us[["Species"]]) ## species code
data.us[["sp"]] <- paste("sp", data.us[["sp"]], sep = ".")
data.us$cluster <- as.character(data.us[["PlotID"]]) ## cluster code
data.us$plot <- paste(as.character(data.us[["PlotID"]]), as.character(data.us[["SubplotNumber"]]),
sep = ".") ## plot code
data.us$htot <- rep(NA, length(data.us[["Species"]])) ## height of tree in m - MISSING
data.us$cluster <- rep(NA,nrow(data.us)) # ## cluster code
data.us$plot <- as.character(data.us[["PlotID"]]) ## plot code
data.us$htot <- rep(NA, length(data.us[["Species"]]))
## height of tree in m - MISSING
data.us$tree.id <- as.character(data.us$TreeID)
data.us$obs.id <- 1:nrow(data.us) ## tree unique id
data.us$sp.name <- species.clean[match(data.us$sp,species.clean$sp),"Latin_name"] ## GET THE SPECIES LATIN NAME species.clean[match(data.us$sp,res.mat[["obs.id"]]),"Latin_name"])
data.us$sp.name <- species.clean[match(data.us$sp,species.clean$sp),
"Latin_name"]
## GET THE SPECIES LATIN NAME
## census is missing use as only one census and check with mark
data.us$census <- rep(1,nrow(data.us))
### add plot weights for computation of competition index (in 1/m^2)
data.us$weights <- 1/(10000 * data.us[["PlotSize"]])
###################### ECOREGION merge greco to have no ecoregion with low number of observation
#### UPDATE NEED TO USE PLOT AND NOT SUBPLOT for competition index based on Mark
## ECOREGION merge greco to have no ecoregion
## with low number of observation
greco <- read.csv(file = "data/raw/US/EcoregionCodes.csv", header = T)
colnames(greco)[1] <- "Ecocode"
table(data.us$Ecocode)
data.us <- merge(data.us, greco[, -4], by = "Ecocode")
data.us$DIVISION <- factor(data.us$DIVISION)
## Some ecoregions still have small # of individuals, so create a variable which
## does division if # ind < 10000 else it reads Domain
data.us$DIVISION <- as.character(data.us$DIVISION)
data.us$DIVISION[data.us$DIVISION %in%
c('Marine Division')] <-
'Warm Continental Division'
data.us$DIVISION[data.us$DIVISION %in%
c('Marine Mountains')] <-
'Warm Continental Mountains'
data.us$DIVISION[data.us$DIVISION %in%
c('Mediterranean Division','Mediterranean Mountains')] <-
'Mediterranean'
data.us$DIVISION[data.us$DIVISION %in%
c('Temperate Desert Division','Temperate Desert Mountains',
'Tropical/Subtropical Desert Division')] <-
'Desert'
data.us$DIVISION[data.us$DIVISION %in%
c('Temperate Steppe Division','Temperate Steppe Mountains',
'Tropical/Subtropical Steppe Division'
,'Tropical/Subtropical Steppe Mountains')] <-
'Steppe'
data.us$DIVISION[data.us$DIVISION %in%
c('Savanna Division','Savanna Mountains')] <-'Savanna'
## Some ecoregions still have small # of individuals, BUT OVERALL FINE
data.us$eco_codemerged <- as.character(data.us$DIVISION)
tab.small.div <- table(data.us$eco_codemerged)
sel.small.div <- which(table(data.us$eco_codemerged) < 10000)
for (i in 1:length(sel.small.div)) {
find.ind <- which(data.us$eco_codemerged == names(tab.small.div)[sel.small.div[i]])
print(length(find.ind))
data.us$eco_codemerged[find.ind] <- as.character(data.us$DOMAIN)[find.ind]
}
#### create good data format to be run per ecoregion
data.us$ecocode <- unlist(lapply(lapply(strsplit(data.us[["eco_codemerged"]], " "),
data.us$ecocode <- unlist(lapply(lapply(strsplit(data.us[["eco_codemerged"]],
" "),
FUN = substr, 1, 2), FUN = paste, collapse = "."))
## ### TEST Olson ecocoregion
## source("R/utils/ecoregions.R")
## data.us$ecocode2 <- GetEcoregions(data.us$Lon,data.us$Lat)
## data.us$ecocode2 <- as.character(data.us$ecocode2)
## ### NEED TO VERIFY THAT ECOREGIONS ARE ASSIGNED PROPERLY
## table(as.character(data.us$ecocode2))
## table(substr(as.character(data.us$ecocode2),1,4))
## table(as.character(data.us$ecocode))
## ## ## plot map to check coordinates syste
## library(RColorBrewer)
......@@ -101,27 +126,31 @@ data.us$ecocode <- unlist(lapply(lapply(strsplit(data.us[["eco_codemerged"]], "
## newmap <- getMap(resolution = 'coarse')
## # different resolutions available
## plot(newmap,xlim=c(-160,-50),ylim=c(45,50))
## ecoreg <- factor(data.us$ecocode);# levels(ecoreg) <- mycols <- brewer.pal(length(levels(ecoreg)), "Set3")
## ecoreg <- factor(data.us$ecocode);# levels(ecoreg) <- mycols <-
## brewer.pal(length(levels(ecoreg)), "Set3")
## points(data.us[['Lon']],data.us[['Lat']],pty='.',cex=.2,col = (ecoreg))
## legend('bottomleft', col = mycols, legend = names(table(data.us$ecocode)), pch = rep(19,length(levels(ecoreg))))
## legend('bottomleft', col = mycols,
## legend = names(table(data.us$ecocode)), pch = rep(19,length(levels(ecoreg))))
###### 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.us[["dead"]], INDEX = data.us[["plot"]], FUN = function.perc.dead)
perc.dead <- tapply(data.us[["dead"]], INDEX = data.us[["plot"]],
FUN = function.perc.dead)
# ## VARIABLE TO SELECT PLOT WITH NOT BIG DISTURBANCE KEEP OFTHER VARIABLES IF
# AVAILABLE (disturbance record)
data.us <- merge(data.us, data.frame(plot = names(perc.dead), perc.dead = perc.dead, stringsAsFactors = FALSE),
by = "plot", sort = FALSE)
data.us <- merge(data.us, data.frame(plot = names(perc.dead),
perc.dead = perc.dead,
stringsAsFactors = FALSE),
by = "plot", sort = FALSE)
##### PLOT SELECTION FOR THE ANALYSIS
## remove everything from memory not need before computation
rm(greco, perc.dead, tab.small.div, sel.small.div)
## variables to keep
vec.abio.var.names <- c("MAT", "MAP")
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","weights","census")
data.tree <- subset(data.us, select = c(vec.basic.var, vec.abio.var.names))
## select only tree above 10cm DBH
......@@ -130,10 +159,11 @@ data.tree <- subset(data.tree, subset =data.tree$D>10 | is.na(data.tree$D))
data.tree <- fun.convert.type.I(data.tree)
write.csv(data.tree,file="output/formatted/US/tree.csv",row.names = FALSE)
### write data plot with variables only at the plot level. I should delete them from the tree table but so far I keep them to not distroy later code
### write data plot with variables only at the plot level.
vec.abio.var.names <- c("MAT", "MAP")
vec.basic.var.p <- c( "cluster", "plot", "ecocode", "Lon", "Lat", "perc.dead")
data.plot <- subset(data.us, subset=!duplicated(data.us$plot),select = c(vec.basic.var.p, vec.abio.var.names))
data.plot <- subset(data.us, subset=!duplicated(data.us$plot),
select = c(vec.basic.var.p, vec.abio.var.names))
write.csv(data.plot,file="output/formatted/US/plot.csv",row.names = FALSE)
rm(data.us)
......
R/process.data/explore.processed.data.R
View file @ f1d3e6cc
......@@ -24,22 +24,26 @@ if(dim(data.all)[1] != sum(mat.perc[['num.obs']]))
fun.plot.hist.trait.per.set(data.all)
to.pdf(fun.hist.var.set(data.all,var='BATOT',cex=0.6),
filename="figs/test.processed/fig.BATOT.set.pdf")
to.pdf(fun.hist.var.set(data.all,var='G',cex=0.6),
filename="figs/test.processed/fig.G.set.pdf")
to.pdf(fun.hist.var.set(data.all,var='BA.G',cex=0.6),
filename="figs/test.processed/fig.BA.G.set.pdf")
to.pdf(fun.hist.var.set(data.all,var='D',cex=0.6),
filename="figs/test.processed/fig.D.set.pdf")
to.dev(fun.plot.xy.set(data.all,var.x='BATOT',var.y='BA.G',cex=0.6),dev=png,
filename="figs/test.processed/fig.xy.BATOT.BA.G.set.png")
to.dev(fun.plot.xy.set(data.all,var.x='D',var.y='BA.G',cex=0.6),dev=png,
filename="figs/test.processed/fig.xy.D.BA.G.set.png")
to.dev(fun.plot.xy.set(data.all,var.x='D',var.y='G',cex=0.6,col=col.vec[data.all$set]),dev=png,
filename="figs/test.processed/fig.xy.D.G.set.png")
to.dev(fun.plot.xy.set(data.all,var.x='BATOT',var.y='G',cex=0.6),dev=png,
filename="figs/test.processed/fig.xy.BATOT.G.set.png")
to.pdf(fun.hist.var.set(data.all, var='BATOT', cex=0.6),
filename = "figs/test.processed/fig.BATOT.set.pdf")
to.pdf(fun.hist.var.set(data.all, var='G', cex=0.6),
filename = "figs/test.processed/fig.G.set.pdf")
to.pdf(fun.hist.var.set(data.all, var='BA.G', cex=0.6),
filename = "figs/test.processed/fig.BA.G.set.pdf")
to.pdf(fun.hist.var.set(data.all, var='D', cex=0.6),
filename = "figs/test.processed/fig.D.set.pdf")
to.dev(fun.plot.xy.set(data.all, var.x = 'BATOT', var.y = 'BA.G',
cex = 0.6), dev = png,
filename = "figs/test.processed/fig.xy.BATOT.BA.G.set.png")
to.dev(fun.plot.xy.set(data.all, var.x = 'D', var.y = 'BA.G', cex = 0.6),
dev = png,
filename = "figs/test.processed/fig.xy.D.BA.G.set.png")
to.dev(fun.plot.xy.set(data.all, var.x = 'D', var.y = 'G',
cex = 0.6, col = col.vec[data.all$set]), dev = png,
filename = "figs/test.processed/fig.xy.D.G.set.png")
to.dev(fun.plot.xy.set(data.all, var.x = 'BATOT', var.y = 'G', cex = 0.6),
dev = png,
filename = "figs/test.processed/fig.xy.BATOT.G.set.png")
########
......@@ -49,14 +53,16 @@ to.dev(fun.plot.xy.set(data.all,var.x='BATOT',var.y='G',cex=0.6),dev=png,
#- pattern of CWM
#- pattern ED angio/conif
## remove outlier following Condit approach
data.all[['G']][!(trim.positive.growth(data.all[['G']]) &
trim.negative.growth(dbh1=data.all[['D']]*10,
dbh2=data.all[['D']]*10 +data.all[['year']]*data.all[['G']]))] <- NA
data.all[['BA.G']][!(trim.positive.growth(data.all[['G']]) &
trim.negative.growth(data.all[['D']]*10,dbh2=data.all[['D']]*10 +data.all[['year']]*data.all[['G']]))] <- NA
## 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)
boxplot(data.all$BATOT ~ data.all$ecocode, las = 3,
at = log(MAP.ECOCODE), boxwex = 0.05, outline = FALSE,
col=col.vec[SET.ECOCODE])
boxplot(data.all$BATOT ~ data.all$set, las = 3)
### compute mean BATOT, number of species, traits and VAR OF TRAITS per cluster
......
R/process.data/plot.shape.var.R
View file @ f1d3e6cc
#######
## SCRIPT TO EXPLORE THE SHAPE OF THE FUNCTION FO RTHE ANALYSIS
## SCRIPT TO EXPLORE THE SHAPE OF THE FUNCTION FOR THE ANALYSIS
source("R/process.data/process.fun.R")
source("R/process.data/test.tree.CWM-fun.R")
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment