updated script for CANADA; uses new dataset

merge.data.CANADA.R

@@ -5,65 +5,59 @@ source("./R/FUN.TRY.R")
 library(reshape)
 
 ######################### READ DATA read individuals tree data
-data.canada <- read.csv("./data/raw/DataCanada/Canada_Data2George_20130818.csv", 
+#data.canada <- read.csv("./data/raw/DataCanada/Canada_Data2George_20130818.csv", 
+#    header = TRUE, stringsAsFactors = FALSE)
+data.canada <- read.csv("./data/raw/DataCanada/Canada_Data2George_20130910.csv", 
     header = TRUE, stringsAsFactors = FALSE)
-data.canada <- data.canada[which(!is.na(data.canada$Species)),]
-colnames(data.canada)[2] <- "sp"
+sum(is.na(data.canada$species))
+data.canada$sp = data.canada$species_FIACode; data.canada$species_FIACode <- NULL
 
-### read species names
+### read species names and merge with data.canada
 species.clean <- read.csv("./data/raw/DataCanada/FIA_REF_SPECIES.csv", stringsAsFactors = FALSE)
+data.canada <- merge(data.canada, data.frame(sp = species.clean$SPCD, species.clean[,2:3], stringsAsFactors = F), by = "sp", all.x = T)
 
 ###################################### 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
 
+data.canada$FinalDBH[data.canada$FinalDBH < 0] <- NA 
 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
+data.canada$year <- data.canada$Interval; data.canada$IntervalYears <- NULL ## number of year between measurement
 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 < 0)  ## dummy variable for dead tree 0 alive 1 dead
-data.canada$plot <- (data.canada[["PLOT_ID"]]); data.canada[["PLOT_ID"]] <- NULL ## plot code
-data.canada$htot <- rep(NA, nrow(data.canada))  ## height of tree in m - MISSING
-data.canada$tree.id <- data.canada$PLOTTREE_I; data.canada$PLOTTREE_I <- NULL ## 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])]
-}
+data.canada$dead <- as.numeric(is.na(data.canada$FinalDBH))  ## dummy variable for dead tree 0 alive 1 dead
+data.canada$htot <- rep(NA, nrow(data.canada))  ## height of tree in m
+data.canada$plot.id <- data.canada$PlotID_InitYear; data.canada$PlotID_InitYear <- NULL ## plot id
+data.canada$tree.id <- data.canada[["PlotTree_ID"]]; data.canada[["PlotTree_ID"]] <- NULL ## tree unique id
+data.canada$census.id <- rep(1,nrow(data.canada))
+data.canada$obs.id <- data.canada$tree.id
 data.canada$weights <- 1/(pi*(0.5*data.canada$D/100)^2)
 
-
-###################### ECOREGION merge greco to have no ecoregion with low number of observation
+###################### ECOREGION merge to have no ecoregion with low number of observation
 greco <- read.csv(file = "./data/raw/DataCanada/EcoregionCodes.csv", header = T, 
     sep = "\t")
 
 table(data.canada$Ecocode)
 ## 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)
-
-# plot(data.canada[['CX']][order(ecoreg)],data.canada[['CY']][order(ecoreg)],pty='.',cex=.2,
-# col = rep(mycols,as.vector(table(ecoreg))))
+## analysis later on or wait for Quebec data to come in 
 
-# legend('bottomright', col = mycols, legend = levels(data.canada$eco_code), pch
-# = rep(19,length(levels(ecoreg))),cex=2) points(data.canada[['CX']][ecoreg ==
-# 9],data.canada[['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.canada$eco_codemerged <- combine_factor(data.canada$eco_code,
-# c(1:8,6,9))
+library(RColorBrewer)
+ecoreg <- factor(data.canada$Ecocode); levels(ecoreg) <- mycols <- brewer.pal(length(levels(ecoreg)), "Set1")
+plot(data.canada[['Lon']],data.canada[['Lat']],pty='.',cex=.2,col = as.character(ecoreg))
+legend('bottomleft', col = mycols, legend = names(table(data.canada$Ecocode)), pch = rep(19,length(levels(ecoreg)))) 
+points(data.canada[['Lon']][data.canada$Ecocode == "-331"], data.canada[['Lat']][data.canada$Ecocode == "-331"],pty=2,cex=1, col = 'black') 
+points(data.canada[['Lon']][data.canada$Ecocode == "-251"], data.canada[['Lat']][data.canada$Ecocode == "-251"],pty=2,cex=1, col = 'red') 
+#data.canada$eco_codemerged <- combine_factor(data.canada$eco_code,c(1:8,6,9))
 
 ###################### 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.canada[["dead"]], INDEX = data.canada[["plot"]], FUN = function.perc.dead2)
-data.canada <- merge(data.canada, data.frame(plot = names(perc.dead), perc.dead = perc.dead), 
-    by = "plot", sort = FALSE)
+perc.dead <- tapply(data.canada[["dead"]], INDEX = data.canada[["plot.id"]], FUN = function.perc.dead2)
+data.canada <- merge(data.canada, data.frame(plot.id = names(perc.dead), perc.dead = perc.dead), 
+    by = "plot.id", sort = FALSE)
 
-########################################################### PLOT SELECTION FOR THE ANALYSIS Remove data with dead == 1
+####################################### PLOT SELECTION FOR THE ANALYSIS Remove data with dead == 1
 table(data.canada$dead)
 data.canada <- data.canada[data.canada$dead == 0,]
 
@@ -98,4 +92,4 @@ if (sum(!data.BA.sp[["tree.id"]] == data.tree[["tree.id"]]) > 0) stop("competiti
 
 ## save everything as a list
 list.canada <- list(data.tree = data.tree, data.BA.SP = data.BA.sp, data.traits = data.traits)
-save(list.spain, file = "./data/process/list.canada.Rdata") 
+save(list.canada, file = "./data/process/list.canada.Rdata")