Newer
Older
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) {
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) {
## 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)
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
## ############################################## 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")