updated R scripts after discussion 190813

merge.data.CANADA.R

@@ -6,8 +6,8 @@ rm(list = ls()); source("./R/format.function.R"); library(reshape)
 ## READ DATA
 ####################
 ### read individuals tree data
-#data.canada <- read.csv("./data/raw/DataCanada/Canada_Data2George_20130815.csv",header=TRUE,stringsAsFactors =FALSE)
-data.canada <- read.csv("./data/raw/DataCanada/Canada_Data2George_20130816.csv",header=TRUE,stringsAsFactors =FALSE)
+#data.canada <- read.csv("./data/raw/DataCanada/Canada_Data2George_20130816.csv",header=TRUE,stringsAsFactors =FALSE)
+data.canada <- read.csv("./data/raw/DataCanada/Canada_Data2George_20130818.csv",header=TRUE,stringsAsFactors =FALSE)
 data.canada <- data.canada[which(!is.na(data.canada$Species)),]
 colnames(data.canada)[2] <- "Species"
 
@@ -18,6 +18,7 @@ species.clean <- read.csv("./data/raw/DataCanada/FIA_REF_SPECIES.csv",stringsAsF
 ## MASSAGE TRAIT DATA
 ############################
 ## HEIGHT DATA FOR TREE MISSING
+## BRING US DATA FOR HEIGHT OVER WHEN WE ANALYZE THAT DATASET LATER ON
 
 ##########################################
 ## FORMAT INDIVIDUAL TREE DATA
@@ -26,19 +27,18 @@ species.clean <- read.csv("./data/raw/DataCanada/FIA_REF_SPECIES.csv",stringsAsF
 ## change unit and names of variables to be the same in all data for the tree 
 data.canada$G <- 10*(data.canada$FinalDBH-data.canada$InitDBH)/data.canada$Interval ## diameter growth in mm per year
 data.canada$G[which(data.canada$InitDBH == 0 | data.canada$FinalDBH == -999)] <- NA
-data.canada$year <- data.canada$Interval ## number of year between measurement - MISSING
+data.canada$year <- data.canada$Interval ## number of year between measuremen
 data.canada$D <- data.canada[["InitDBH"]]; data.canada$D[data.canada$D == 0] <- NA ;## diameter in cm
-data.canada$dead <- as.numeric(data.canada$FinalDBH == -999) ## dummy variable for dead tree 0 alive 1 dead
+data.canada$dead <- as.numeric(data.canada$FinalDBH > 0) ## dummy variable for dead tree 0 alive 1 dead
 data.canada$sp <- as.character(data.canada[["Species"]]) ## species code
 data.canada$plot <- (data.canada[["PLOT_ID"]]) ## plot code
 data.canada$htot <- rep(NA,length(data.canada[["Species"]])) ## height of tree in m - MISSING
-data.canada$tree.id <- data.canada$PLOTTREE ## tree unique id
+data.canada$tree.id <- gsub("_",".",data.canada$PLOTTREE); ## tree unique id
 data.canada$sp.name <- NA; 
 for(i in 1:length(unique(data.canada$sp))) { 
 	v <- species.clean$SPCD
 	data.canada$sp.name[which(data.canada$sp == unique(data.canada$sp)[i])] <- species.clean$COMMON_NAME[which(v == unique(data.canada$sp)[i])] }
 
-
 ######################
 ## ECOREGION
 ###################
@@ -46,7 +46,7 @@ for(i in 1:length(unique(data.canada$sp))) {
 greco <- read.csv(file = "./data/raw/DataCanada/EcoregionCodes.csv", header = T, sep = "\t") 
 
 table(data.canada$Ecocode) 
-## There is only four ecoregions though; do you want to aggregate into divsion still?
+## Some ecoregions still have small # of individuals, so either drop off for analysis later on or wait for Quebec data to come in
 # 
 # library(RColorBrewer); mycols <- brewer.pal(10,"Set3"); 
 # ecoreg <- unclass(data.canada$eco_code); 

merge.data.NZ.R

@@ -17,7 +17,7 @@ data.nz$plid <- gsub("__", "_",data.nz$plid); data.nz$plid <- gsub("_", ".",data
 ############################
 ## Maximum height per species is already available from data.trait (in m); so no sd's and only obs per spp
 data.max.height <- data.frame(code=data.trait[["sp"]],Max.height.mean=log10(data.trait[["height.m"]]),Max.height.sd=NA,Max.height.nobs=1)
-write.csv(data.max.height,file="./data/process/data.max.height.nz.csv") ## I was planning to save processed data in that folder
+#write.csv(data.max.height,file="./data/process/data.max.height.nz.csv") ## I was planning to save processed data in that folder
 
 ################################################################
 ## FORMAT INDIVIDUAL TREE DATA
@@ -27,10 +27,10 @@ write.csv(data.max.height,file="./data/process/data.max.height.nz.csv") ## I was
 data.nz$G <- 10*(data.nz[["D1"]]-data.nz[["D0"]])/(data.nz[["t1"]]-data.nz[["t0"]]) ## diameter growth in mm per year
 data.nz$year <- (data.nz[["t1"]]-data.nz[["t0"]]) ## number of year between measurement
 data.nz$D <- data.nz[["D0"]] ## diameter in mm convert to cm
-data.nz$dead <- rep(NA,length(data.nz[["sp"]])) ## dummy variable for dead tree 0 alive 1 dead - MISSING
+data.nz$dead <- as.numeric(is.na(data.nz[["D1"]])) ## dummy variable for dead tree 0 alive 1 dead
 data.nz$sp <- data.nz$sp
-data.nz$plot <- data.nz$plid 
-data.nz$htot <- rep(NA,length(data.nz[["sp"]])) ## I ASSUME THE "height.m" IN data.trait IS A MAX HEIGHT, WHICH IS NOT WHAT WE WANT HERE?
+data.nz$plot <- data.nz$plid
+data.nz$htot <- rep(NA,length(data.nz[["sp"]])) ## Max height is already available so have as missing
 data.nz$tree.id <- data.nz[["tag"]]; data.nz$tree.id <- gsub("__", "_",data.nz$tree.id); data.nz$tree.id <- gsub("_", ".",data.nz$tree.id) ## tree unique id
 
 ##########################################
@@ -59,22 +59,24 @@ rm(data.sp,data.sp2)
 ######################
 ## ECOREGION
 ###################
-## DON'T KNOW WHICH VARIABLE TO USE FOR ECOREGION
-greco <- greco[,c("tree.id","BIOME","eco_code")]
+## merge greco to have no ecoregion with low number of observation
+data.nz <- merge(data.nz, data.frame(plot = data.plot$plid, data.plot[,11:12]), by = "plot")
+table(data.nz$Broad); 
 
 ######################
 ## PERCENT DEAD
 ###################
-## variable percent dead - MISSING
-data.nz$perc.dead <- NA
+perc.dead <- tapply(data.nz[["dead"]],INDEX=data.nz[["plot"]],FUN=function.perc.dead)
+# ## VARIABLE TO SELECT PLOT WITH NOT BIG DISTURBANCE KEEP OFTHER VARIABLES IF AVAILABLE (disturbance record)
+data.nz <- merge(data.nz,data.frame(plot=names(perc.dead),perc.dead=perc.dead), by = "plot", sort=FALSE)
 
 ###########################################################
 ### PLOT SELECTION FOR THE ANALYSIS
 ###################
 data.nz <- merge(data.nz, data.plot[,c(1,8:10)], by = "plid")
-colnames(data.nz)[colnames(data.nz) %in% c("mat","map")] <- c("MAT","PP")
+colnames(data.nz)[colnames(data.nz) %in% c("mat","map")] <- c("MAT","MAP")
 #colnames(data.nz)[names(data.nz) =="eco_codemerged"  ] <- c("ecocode")
-vec.abio.var.names <-  c("MAT","PP")
+vec.abio.var.names <-  c("MAT","MAP")
 vec.basic.var <-  c("tree.id","sp","spname","plot","ecocode","D","G","dead","year","htot","Lon","Lat","perc.dead")
 data.tree <- subset(data.nz,select=c(vec.basic.var,vec.abio.var.names))
 
@@ -86,8 +88,7 @@ data.BA.SP <- BA.SP.FUN(id.tree=as.vector(data.nz[["tree.id"]]), diam=as.vector(
 	weights=as.vector(1/(pi*(0.5*data.nz[["D0"]])^2)), weight.full.plot=NA)
 
 ## change NA and <0 data for 0
-data.BA.SP[is.na(data.BA.SP)] <- 0
-data.BA.SP[,-1][data.BA.SP[,-1]<0] <- 0
+data.BA.SP[is.na(data.BA.SP)] <- 0; data.BA.SP[,-1][data.BA.SP[,-1]<0] <- 0
 
 ### CHECK IF sp and sp name for column are the same
 if(sum(!(names(data.BA.SP)[-1] %in% unique(data.nz[["sp"]]))) >0) stop("competition index sp name not the same as in data.tree")

merge.data.SPAIN.R

@@ -86,6 +86,7 @@ rm(greco2)
 
 table(data.spain$eco_code)
 ## There's an eco-region with no code, and one with 55 sites
+## The former we could drop as they were on the border of Spain
 
 library(RColorBrewer); mycols <- brewer.pal(10,"Set3"); 
 ecoreg <- unclass(data.spain$eco_code); 
@@ -94,6 +95,7 @@ legend("bottomright", col = mycols, legend = levels(data.spain$eco_code), pch =
 points(data.spain[["CX"]][ecoreg == 9],data.spain[["CY"]][ecoreg == 9],pty=".",cex=.2, col = "black"); ## Highlight the region with 55 sites
 ## PA1219 looks to be similar to PA1209; merge them together
 data.spain$eco_codemerged <- combine_factor(data.spain$eco_code, c(1:8,6,9))
+data.spain <- data.spain[-which(data.spain$eco_codemerged == ""),]
 
 #######################
 ## variable percent dead/cannot do with since dead variable is missing