modified canada R script for Canada_Data2George_20130816

merge.data.CANADA-fhv1.R → merge.data.CANADA.R

@@ -6,7 +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_20130815.csv",header=TRUE,stringsAsFactors =FALSE)
+data.canada <- read.csv("./data/raw/DataCanada/Canada_Data2George_20130816.csv",header=TRUE,stringsAsFactors =FALSE)
 data.canada <- data.canada[which(!is.na(data.canada$Species)),]
 colnames(data.canada)[2] <- "Species"
 
@@ -31,13 +32,14 @@ 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 <- (data.canada[["FinalDBH"]]-data.canada[["InitDBH"]])/data.canada$Interval ## diameter growth in mm per year
-data.canada$year <- data.canada$Interval ## number of year between measurement/missing!
-data.canada$D <- data.canada[["InitDBH"]] ## diameter in mm
-data.canada$dead <- rep(NA,length(data.canada[["Species"]])) ## dummy variable for dead tree 0 alive 1 dead/missing!
+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$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$sp <- as.character(data.canada[["Species"]]) ## species code
-data.canada$plot <- (data.canada[["PlotID"]]) ## plot code
-data.canada$htot <- rep(NA,length(data.canada[["Species"]]))## height of tree in m / missing
+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$sp.name <- NA; 
 for(i in 1:length(unique(data.canada$sp))) { 
@@ -47,16 +49,10 @@ for(i in 1:length(unique(data.canada$sp))) {
 
 ############################
 ## merge greco to have no ecoregion with low number of observation
-# greco <- read.csv(file = "./data/raw/DataSpain/R_Ecoregion.csv", header = T) 
-# greco 	<- greco[,c("Plot_ID_SFI","BIOME","eco_code")]
-# greco2 <- greco[!duplicated(greco$Plot),]; 
-# rm(greco)
-# 
-# data.canada <- merge(data.canada, greco2, by = "Plot_ID_SFI")
-# rm(greco2)
-# 
-# table(data.canada$eco_code)
-# ## There's an eco-region with no code, and one with 55 sites
+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?
 # 
 # library(RColorBrewer); mycols <- brewer.pal(10,"Set3"); 
 # ecoreg <- unclass(data.canada$eco_code); 
@@ -67,22 +63,21 @@ for(i in 1:length(unique(data.canada$sp))) {
 # data.canada$eco_codemerged <- combine_factor(data.canada$eco_code, c(1:8,6,9))
 
 #######################
-# ## 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.canada[["dead"]],INDEX=data.canada[["idp"]],FUN=function.perc.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.canada[["dead"]],INDEX=data.canada[["plot"]],FUN=function.perc.dead)
 # ## VARIABLE TO SELECT PLOT WITH NOT BIG DISTURBANCE KEEP OFTHER VARIABLES IF AVAILABLE (disturbance record)
-# data.canada <- merge(data.canada,data.frame(idp=as.numeric(names(perc.dead)),perc.dead=perc.dead),sort=FALSE)
-data.canada$perc.dead <- NA
+data.canada <- merge(data.canada,data.frame(plot=names(perc.dead),perc.dead=perc.dead), by = "plot", sort=FALSE)
 
 ###########################################################
 ### PLOT SELECTION FOR THE ANALYSIS
 ###################
-# ## Remove data with mortality == 1 or 2
-# table(data.canada$Mortality_Cut)
-# data.canada <- subset(data.canada,subset= (data.canada[["Mortality_Cut"]] == 0 |  data.canada[["Mortality_Cut"]] == ""))
+## Remove data with dead == 1
+table(data.canada$dead)
+## Nothing to remove
 
-colnames(data.canada)[c(2,5,10,12,14)] <- c("sp","plot","w","ecocode","PP")
-vec.abio.var.names <-  c("MAT","PP")
+colnames(data.canada)[c(3,1,11,13)] <- c("sp","plot","w","ecocode")
+vec.abio.var.names <-  c("MAT","MAP")
 vec.basic.var <-  c("tree.id","sp","sp.name","plot","ecocode","D","G","dead","year","htot","Lon","Lat","perc.dead")
 data.tree <- subset(data.canada,select=c(vec.basic.var,vec.abio.var.names))
 
@@ -91,17 +86,16 @@ data.tree <- subset(data.canada,select=c(vec.basic.var,vec.abio.var.names))
 ###########################
 data.BA.SP <- BA.SP.FUN(id.tree=as.vector(data.canada[["tree.id"]]), diam=as.vector(data.canada[["D"]]),
 	sp=as.vector(data.canada[["sp"]]), id.plot=as.vector(data.canada[["plot"]]),
-	weights=1/(10000*data.canada[["SubPlotSize"]]), weight.full.plot=NA)
+	weights=1/(10000*data.canada[["SubPlot_Size"]]), 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.canada[["sp"]]))) >0) stop("competition index sp name not the same as in data.tree")
 
 #### compute BA tot for all competitors
-BATOT.COMPET <- apply(data.BA.SP[,-1],MARGIN=1,FUN=sum,na.rm=TRUE)
+BATOT.COMPET <- apply(data.BA.SP[,-1],1,sum,na.rm=TRUE)
 data.BA.SP$BATOT.COMPET <- BATOT.COMPET; rm(BATOT.COMPET)
 ### create data frame
 names(data.BA.SP) <- c("tree.id",names(data.BA.SP)[-1])