diff --git a/merge.data.CANADA.R b/merge.data.CANADA.R index bb34d7932bc46fe2d24c830ee7c81680c65fe512..46571a4b023293bba6296a75f4687e9902ea9450 100644 --- a/merge.data.CANADA.R +++ b/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")