merge.data.BCI.R 5.66 KiB
### MERGE BCI DATA Edited by FH
rm(list = ls())
source("./R/format.function.R")
library(reshape)
######################### READ DATA read individuals tree data Requires careful formatting of 7 census
######################### datasets The raw data is such that, once a tree dies in census X, then it no
######################### longer exists in census X+1, X+2 etc...
data.bci1 <- read.table("./data/raw/DataBCI/census1/PlotsDataReport.txt", header = TRUE, 
    stringsAsFactors = FALSE, sep = "\t")
big.bci <- NULL
for (k in 2:7) {
    new.directory <- paste("./data/raw/DataBCI/census", k, "/PlotsDataReport.txt", 
        sep = "")
    data.bci2 <- read.table(new.directory, header = TRUE, stringsAsFactors = FALSE, 
        sep = "\t")
    if (!is.null(big.bci)) {
        sub.bci <- merge(data.bci1[, c(2:7, 11, 13)], data.frame(TreeID = data.bci2[["TreeID"]], 
            DBH2 = data.bci2[["DBH"]], Date2 = data.bci2[["Date"]], dead = as.numeric(data.bci2[["Status"]] == 
                "dead")), sort = T, by = "TreeID")
        big.bci <- rbind(big.bci, sub.bci)
    if (is.null(big.bci)) {
        big.bci <- merge(data.bci1[, c(2:7, 11, 13)], data.frame(TreeID = data.bci2[["TreeID"]], 
            DBH2 = data.bci2[["DBH"]], Date2 = data.bci2[["Date"]], dead = as.numeric(data.bci2[["Status"]] == 
                "dead")), sort = T, by = "TreeID")
    data.bci1 <- data.bci2
    cat("Census", k, "now included\n")
rm(data.bci1, data.bci2, sub.bci)
big.bci <- big.bci[order(big.bci$TreeID), ]
colnames(big.bci)[c(7:8)] <- c("DBH1", "Date1")
data.bci <- big.bci
rm(big.bci)
### read species names
species.clean <- read.table("./data/raw/DataBCI/TaxonomyDataReport.txt", stringsAsFactors = FALSE, 
    header = T, sep = "\t")
## Try to relate SpeciesID in species.clean species names in data.bci
data.bci$sp2 = species.clean$SpeciesID[match(data.bci$Latin, paste(species.clean[["Genus"]], 
    species.clean[["species"]]))]
length(unique(data.bci$sp))
###################################### MASSAGE TRAIT DATA Use HEIGHT_AVG, LMALAM_AVD, SEED_DRY, BUT I DO NOT KNOW
###################################### WHICH WOOD DENSITY VARIABLE TO USE
data.trait <- read.csv("./data/raw/DataBCI/BCITRAITS_20101220.csv", stringsAsFactors = FALSE, 
    header = T)
data.trait$Latin <- apply(data.trait[, 1:2], 1, paste, collapse = " ")
data.bci <- merge(data.bci, data.trait[, c(ncol(data.trait), 3, 7:10, 13, 15, 18, 
    20:21)], by = "Latin", all.x = T)
data.bci <- data.bci[order(data.bci$TreeID), ]
########################################## FORMAT INDIVIDUAL TREE DATA
data.bci$Date1 <- as.Date(data.bci$Date1)
data.bci$Date2 <- as.Date(data.bci$Date2)
# data.bci$yr1 <- format(strptime(data.bci$Date1, format = '%Y-%m-%d'),'%Y')
# data.bci$yr2 <- format(strptime(data.bci$Date2, format = '%Y-%m-%d'),'%Y')
data.bci$year <- as.numeric(difftime(data.bci$Date2, data.bci$Date1, units = "weeks")/52)  ## Not rounded 
## change unit and names of variables to be the same in all data for the tree
data.bci$G <- 10 * (data.bci$DBH1 - data.bci$DBH1)/data.bci$year  ## diameter growth in mm per year - BASED ON UNROUNDED YEARS
data.bci$D <- data.bci[["DBH1"]]
data.bci$plot <- data.bci[["Quadrat"]]  ## plot code?
data.bci$htot <- data.bci$HEIGHT_AVG
data.bci$sp.name <- data.bci$Latin
###################### ECOREGION bci has only 1 eco-region
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120
###################### PERCENT DEAD variable percent dead/cannot do with since dead variable is ###################### missing compute numer of dead per plot to remove plot with disturbance function.perc.dead2 <- function(dead) { out <- sum(dead, na.rm = T)/length(dead[!is.na(dead)]) if (!is.finite(out)) out <- NA return(out) } perc.dead <- tapply(data.bci[["dead"]], INDEX = data.bci[["plot"]], FUN = function.perc.dead2) data.bci <- merge(data.bci, data.frame(plot = names(perc.dead), perc.dead = perc.dead), by = "plot", sort = FALSE) ########################################################### PLOT SELECTION FOR THE ANALYSIS Remove data with dead == 1 table(data.bci$dead) vec.abio.var.names <- c("MAT", "MAP") ## MISSING vec.basic.var <- c("treeid", "sp", "sp.name", "plot", "D", "G", "dead", "year", "htot", "x", "y", "perc.dead") data.tree <- subset(data.bci, select = c(vec.basic.var, vec.abio.var.names)) ############################################## COMPUTE MATRIX OF COMPETITION INDEX WITH SUM OF BA PER SPECIES IN EACH PLOT in ############################################## m^2/ha without the target species DON'T KNOW SUBPLOT SIZE! data.BA.SP <- BA.SP.FUN(id.tree = as.vector(data.bci[["TreeID"]]), diam = as.vector(data.bci[["D"]]), sp = as.vector(data.bci[["sp"]]), id.plot = as.vector(data.bci[["plot"]]), weights = 1/(pi * (0.5 * data.bci$D/100)^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 ### CHECK IF sp and sp name for column are the same if (sum(!(names(data.BA.SP)[-1] %in% unique(data.bci[["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], 1, sum, na.rm = TRUE) data.BA.SP$BATOT.COMPET <- BATOT.COMPET rm(BATOT.COMPET) ### create data frame names(data.BA.SP) <- c("TreeID", names(data.BA.SP)[-1]) data.BA.sp <- merge(data.frame(TreeID = data.bci[["TreeID"]], ecocode = data.bci[["ecocode"]]), data.BA.SP, by = "TreeID", sort = FALSE) ## test if (sum(!data.BA.sp[["TreeID"]] == data.tree[["TreeID"]]) > 0) stop("competition index not in the same order than data.tree") ## save everything as a list list.bci <- list(data.tree = data.tree, data.BA.SP = data.BA.sp, data.traits = data.traits) save(list.bci, file = "./data/process/list.bci.Rdata")