#!/usr/bin/env Rscript

### MERGE NSW DATA
rm(list = ls())
source("R/format.data/format.fun.R")
dir.create("output/formatted/NSW", recursive=TRUE,showWarnings=FALSE)

library(reshape)
######################### READ DATA read individuals tree data
data.nswbrc <- read.csv("data/raw/NSW/NSW_data_BRcontrols.csv", header = TRUE, 
    stringsAsFactors = FALSE)
data.nswbrc$Date.of.measure <- as.vector(sapply(data.nswbrc$Date.of.measure, function(x) {
    unlist(strsplit(x, "/"))[3]}))  ## Extract the years only


## write a function that create a 10*10m quadrat per plot
fun.quadrat <- function(x,y,max.x,max.y){
vec.x <- seq(0,max.x,by=10)
vec.y <- seq(0,max.y,by=10)
x.cut <- cut(x,breaks=vec.x)
y.cut <- cut(y,breaks=vec.y)

return(unclass(factor(paste(x.cut,y.cut))))

}



data.nswbrt <- read.csv("data/raw/NSW/NSW_data_BRtreatments.csv", header = TRUE, 
    stringsAsFactors = FALSE)
data.nswbrt$Date.of.measure <- as.vector(sapply(data.nswbrt$Date.of.measure, function(x) {
    unlist(strsplit(x, "/"))[3]}))  ## Extract the years only
data.nswbs1 <- read.csv("data/raw/NSW/NSW_data_BS1.csv", header = TRUE, stringsAsFactors = FALSE)

data.nswbs1$Date.of.measure <- as.character(format(as.Date(data.nswbs1$Date.of.measure, 
    format = "%d-%b-%y"), format = "%d/%b/%Y"))
data.nswbs1$Date.of.measure <- as.vector(sapply(data.nswbs1$Date.of.measure, function(x) {
    unlist(strsplit(x, "/"))[3]}))  ## Extract the years only
## data.nswbs2 has a different format to the other datasets, so format to match
## the above
data.nswbs2 <- read.csv("data/raw/NSW/NSW_data_BS2.csv", header = TRUE, stringsAsFactors = FALSE)
data.nswbs2$Plot <- apply(data.nswbs2[, 1:2], 1, paste, collapse = "")
data.nswbs2$Subplot <- NULL
data.nswbs2$Family <- NULL
colnames(data.nswbs2)[3] <- "species"
data.nswbs2$x <- data.nswbs2$y <- rep(NA, nrow(data.nswbs2))
data.nswbs22 <- data.nswbs2[rep(1:nrow(data.nswbs2), each = 2), ]
data.nswbs22$Date.of.measure <- rep(c(1988, 2000), nrow(data.nswbs2))
data.nswbs22$Dbh <- c(rbind(data.nswbs2[["DBH.cm..1988."]], data.nswbs2[["DBH.cm..2000."]]))
data.nswbs22[["DBH.cm..1988."]] <- data.nswbs22[["DBH.cm..2000."]] <- NULL
data.nswbs2 <- data.nswbs22[, c(1:5, 8:9, 7, 6)]
rm(data.nswbs22)
data.nswtnd <- read.csv("data/raw/NSW/NSW_data_TND.csv", header = TRUE, stringsAsFactors = FALSE)
data.nswtnd$Date.of.measure <- as.character(format(as.Date(data.nswtnd$Date.of.measure, 
    format = "%d-%b-%y"), format = "%d/%b/%Y"))
data.nswtnd$Date.of.measure <- as.vector(sapply(data.nswtnd$Date.of.measure, function(x) {
    unlist(strsplit(x, "/"))[3]
}))  ## Extract the years only
data.nsw <- rbind(data.nswbrc, data.nswbrt, data.nswbs1, data.nswbs2, data.nswtnd)


######
## NEED TO CREATE plot of 10x10m for all tree with either x y coordinates or subplot in BS2



###################################### MASSAGE TRAIT DATA
data.trait <- read.csv("data/raw/NSW/NSW_traits.csv", header = TRUE, stringsAsFactors = FALSE)
data.trait$sp <- data.trait[["Species.all"]]; data.trait[["Species.all"]] <- NULL ## There is not sp.code in data.nsw; using spp name as code
data.trait$Latin_name <- data.trait$sp
data.trait$Leaf.N.mean <- NA
data.trait$Leaf.N.sd <- NA
data.trait$Seed.mass.mean <- (10^data.trait$SDM_log10_g)*1000; data.trait$SDM_log10_g <- NULL ## conversion from log10 g to mg
data.trait$Seed.mass.sd <- NA
data.trait$SLA.mean <- NA
data.trait$SLA.sd <- NA
data.trait$Wood.density.mean <- data.trait[["WD_basic_kg.m3"]]/1000; data.trait[["WD_basic_kg.m3"]] <- NULL ## conversion from kg/m3 to mg/mm3
data.trait$Wood.density.sd <- NA
data.trait$Max.height.mean <- data.trait$Log10_Hmax_m; data.trait$Log10_Hmax_m <- NULL
data.trait$Max.height.sd <- NA


########################################## FORMAT INDIVIDUAL TREE DATA Each tree has at most 3 observations (from prelim
########################################## checks of the data)
data.nsw$tree.id <- apply(data.nsw[, 1:2], 1, paste, collapse = "_")
data.nsw2 <- data.frame(data.nsw[1, ], year1 = NA, year2 = NA, dbh1 = NA, dbh2 = NA)
for (k in 1:length(unique(data.nsw$tree.id))) {
    sub.datansw <- as.data.frame(data.nsw[which(data.nsw$tree.id == unique(data.nsw$tree.id)[k]), 
        ])
    if (nrow(sub.datansw) == 1) {
        data.nsw2 <- rbind(data.nsw2, data.frame(sub.datansw, year1 = sub.datansw$Date.of.measure[1], 
            year2 = NA, dbh1 = sub.datansw$Dbh[1], dbh2 = NA, stringsAsFactors = F))
    }
    if (nrow(sub.datansw) == 2) {
        data.nsw2 <- rbind(data.nsw2, data.frame(sub.datansw[1, ], year1 = sub.datansw$Date.of.measure[1], 
            year2 = sub.datansw$Date.of.measure[2], dbh1 = sub.datansw$Dbh[1], dbh2 = sub.datansw$Dbh[2], stringsAsFactors = F))
    }
    if (nrow(sub.datansw) == 3) {
        data.nsw2 <- rbind(data.nsw2, data.frame(sub.datansw[1, ], year1 = sub.datansw$Date.of.measure[1], 
            year2 = sub.datansw$Date.of.measure[2], dbh1 = sub.datansw$Dbh[1], dbh2 = sub.datansw$Dbh[2], stringsAsFactors = F))
        data.nsw2 <- rbind(data.nsw2, data.frame(sub.datansw[1, ], year1 = sub.datansw$Date.of.measure[2], 
            year2 = sub.datansw$Date.of.measure[3], dbh1 = sub.datansw$Dbh[2], dbh2 = sub.datansw$Dbh[3], stringsAsFactors = F))
    }
}
data.nsw2 <- data.nsw2[-1, ]
data.nsw2$Date.of.measure <- data.nsw2$Dbh <- NULL
data.nsw <- data.nsw2
for (k in 9:12) data.nsw[, k] <- as.numeric(data.nsw[, k])
## change unit and names of variables to be the same in all data for the tree
data.nsw$year <- (data.nsw$year2 - data.nsw$year1)  ## number of year between measurements
data.nsw$G <- 10 * (data.nsw$dbh2 - data.nsw$dbh1)/(data.nsw$year)  ## diameter growth in mm per year
#head(data.nsw[order(data.nsw$G), ])
data.nsw$D <- data.nsw[["dbh1"]] ## diameter in cm
data.nsw$dead <- rep(NA, nrow(data.nsw))  ## dummy variable for dead tree 0 alive 1 dead - MISSING USE M
data.nsw$sp <- as.numeric(factor(data.nsw[["species"]])); ## species code - use the spp name as code
data.nsw$sp.name <- data.nsw[["species"]]; data.nsw$species <- NULL
data.nsw$plot <- as.character(data.nsw[["Plot"]]); data.nsw$Plot <- NULL ## plot code
data.nsw$subplot <- rep(NA,nrow(data.nsw)) 
data.nsw$census <- rep(NA,nrow(data.nsw)) 
data.nsw$htot <- rep(NA, nrow(data.nsw))  ## height of tree in m
### add plot weights for computation of competition index (in 1/m^2) - from the
### original excel file
data.nsw$weights <- rep(NA, nrow(data.nsw))
data.nsw$weights[grep("AA", data.nsw$Plot)] <- 1/(20 * 80)
data.nsw$weights[grep("BB", data.nsw$Plot)] <- 1/(20 * 80)
data.nsw$weights[grep("CC", data.nsw$Plot)] <- 1/(20 * 80)
data.nsw$weights[grep("DD", data.nsw$Plot)] <- 1/(20 * 80)
data.nsw$weights[grep("BS", data.nsw$Plot)] <- 1/(25 * 30)
data.nsw$weights[grep("BR", data.nsw$Plot)] <- 1/(60.4 * 60.4)
data.nsw$weights[grep("END", data.nsw$Plot)] <- 1/(40 * 50)
data.nsw$weights[grep("TND", data.nsw$Plot)] <- 1/(40 * 50)
data.nsw$cenus <- data.nsw$year1
data.nsw$obs.id <- 1:nrow(data.nsw)
###################### ECOREGION nsw has only 1 eco-region
###################### PERCENT DEAD NO DATA ON MORTALITY
perc.dead <- tapply(data.nsw[["dead"]], INDEX = data.nsw[["plot"]], FUN = function.perc.dead2)
data.nsw <- merge(data.nsw, data.frame(plot = names(perc.dead), perc.dead = perc.dead), 
    by = "plot", sort = FALSE)
########################################### VARIABLES SELECTION FOR THE ANALYSIS
vec.basic.var <- c("obs.id","tree.id", "sp", "sp.name","plot", "subplot", "D", "G", "dead", 
    "year", "htot", "x", "y", "census")
data.tree <- subset(data.nsw, select = c(vec.basic.var))
write.csv(data.tree,file="output/formatted/NSW/tree.csv",row.names = FALSE)


## ############################################## COMPUTE MATRIX OF COMPETITION INDEX WITH SUM OF BA PER SPECIES IN EACH PLOT in
## ############################################## m^2/ha without the target species NEED TO COMPUTE BASED ON RADIUS AROUND TARGET
## ############################################## TREE species as factor because number
## data.nsw[["species"]] <- factor(data.nsw[["species"]])
## data.nsw$spcode <- data.nsw[["species"]]
## levels(data.nsw$spcode) <- 1:length(levels(data.nsw$spcode))
## data.BA.SP <- BA.SP.FUN(obs.id = as.vector(data.nsw[["tree.id"]]), diam = as.vector(data.nsw[["D"]]), 
##     sp = as.vector(data.nsw[["spcode"]]), id.plot = as.vector(data.nsw[["Plot"]]), 
##     weights = data.nsw[["weights"]], weight.full.plot = NA)
## ## change NA and <0 data for 0
## data.BA.SP[is.na(data.BA.SP)] <- 0
## data.BA.SP2 <- data.frame(data.BA.SP)
## colnames(data.BA.SP2) <- colnames(data.BA.SP)
## ### CHECK IF sp and sp name for column are the same
## if (sum(!(names(data.BA.SP2)[-1] %in% unique(data.nsw[["spcode"]]))) > 0) stop("competition index sp name not the same as in data.tree")
## #### compute BA tot for all competitors
## BATOT.COMPET <- apply(data.BA.SP2[, -1], 1, sum, na.rm = TRUE)
## data.BA.SP2$BATOT.COMPET <- BATOT.COMPET
## rm(BATOT.COMPET)
## data.BA.SP <- data.BA.SP2
## # Rlim <- 15 # set size of neighborhood for competition index system.time(test <-
## # fun.compute.BA.SP.XY.per.plot(1,data.tree=data.nsw,Rlim=15,parallel=TRUE,rpuDist=FALSE))
## # list.BA.SP.data <-
## # mclapply(unique(data.nsw[['plot']]),FUN=fun.compute.BA.SP.XY.per.plot,data.tree=data.nsw,Rlim=Rlim,mc.cores=4)
## # data.BA.sp <- rbind.fill(list.BA.SP.data) dim(data.BA.SP) ### TEST DATA FORMAT
## # if(sum(! rownames(BA.SP.temp)==data.tree[['obs.id']]) >0) stop('rows not in the
## # good order')
## # if(sum(!colnames(BA.SP.temp)==as.character((levels(data.tree[['species']]))))>0)
## # stop('colnames does mot match species name') ## test same order as data.nsw
## # if(sum(!data.BA.SP[['obs.id']] == data.nsw[['obs.id']]) >0) stop('competition
## # index not in the same order than data.nsw') ## REMOVE TREE IN BUFFER ZONE
## # not.in.buffer.zone <- (data.nsw[['x']]<(250-Rlim) & data.nsw[['x']]>(0+Rlim) &
## # data.nsw[['y']]<(250-Rlim) & data.nsw[['y']]>(0+Rlim)) # remove subset data.nsw
## # <- subset(data.nsw,subset=not.in.buffer.zone) data.BA.sp <-
## # subset(data.BA.sp,subset=not.in.buffer.zone) ## plot each plot
## # pdf('./figs/plots.tree.pdf')
## # lapply(unique(data.nsw[['plot']]),FUN=fun.circles.plot,data.nsw[['x']],data.nsw[['y']],data.nsw[['plot']],data.nsw[['D']],inches=0.2,xlim=c(0,250),ylim=c(0,250))
## # dev.off()
## ## save everything as a list
## list.nsw <- list(data.tree = data.nsw, data.BA.SP = data.BA.sp, data.traits = data.traits)
## save(list.nsw, file = "output/formatted/NSW/list.nsw.Rdata")