-
Daniel Falster authored7dbee2ef
### MERGE NZ DATA
rm(list = ls())
source("./R/format.function.R")
library(reshape)
######################### READ DATA read individuals tree data
data.nz <- read.csv("./data/raw/DataNVS/nz_treedata_growth_130801.csv", header = TRUE,
stringsAsFactors = FALSE, skip = 9)
data.nz <- data.nz[, -1]
data.trait <- read.csv("./data/raw/DataNVS/nz_traits_130801.csv", , header = TRUE,
stringsAsFactors = FALSE)
data.trait <- data.trait[, -1]
data.plot <- read.csv("./data/raw/DataNVS/nz_plotinfo_130801.csv", , header = TRUE,
stringsAsFactors = FALSE)
data.plot <- data.plot[, -1]
colnames(data.trait)[1] <- "sp"
data.plot$plid <- gsub("__", "_", data.plot$plid)
data.plot$plid <- gsub("_", ".", data.plot$plid) ## Replace all underscores with a single dot
data.nz$plid <- gsub("__", "_", data.nz$plid)
data.nz$plid <- gsub("_", ".", data.nz$plid) ## Replace all underscores with a single dot
###################################### MASSAGE TRAIT DATA Maximum height per species is already available from
###################################### data.trait (in m); so no sd's and only obs per spp
data.max.height <- data.frame(code = data.trait[["sp"]], Max.height.mean = log10(data.trait[["height.m"]]),
Max.height.sd = NA, Max.height.nobs = 1)
# write.csv(data.max.height,file='./data/process/data.max.height.nz.csv') ## I
# was planning to save processed data in that folder
################################################################ FORMAT INDIVIDUAL TREE DATA
## change unit and names of variables to be the same in all data for the tree
data.nz$G <- 10 * (data.nz[["D1"]] - data.nz[["D0"]])/(data.nz[["t1"]] - data.nz[["t0"]]) ## diameter growth in mm per year
data.nz$year <- (data.nz[["t1"]] - data.nz[["t0"]]) ## number of year between measurement
data.nz$D <- data.nz[["D0"]] ## diameter in mm convert to cm
data.nz$dead <- as.numeric(is.na(data.nz[["D1"]])) ## dummy variable for dead tree 0 alive 1 dead
data.nz$sp <- data.nz$sp
data.nz$plot <- data.nz$plid
data.nz$htot <- rep(NA, length(data.nz[["sp"]])) ## Max height is already available so have as missing
data.nz$tree.id <- data.nz[["tag"]]
data.nz$tree.id <- gsub("__", "_", data.nz$tree.id)
data.nz$tree.id <- gsub("_", ".", data.nz$tree.id) ## tree unique id
########################################## CHANGE COORDINATE SYSTEM DON'T KNOW THE EPSG CODE HERE change coordinates
########################################## system of Easting Northing to be in lat long WGS84
data.nz <- merge(data.plot[, c(1, 3, 4)], data.nz, by = "plid")
library(sp)
# library(dismo); library(rgdal);
data.sp <- data.nz[, c("tree.id", "Easting", "Northing")]
coordinates(data.sp) <- c("Easting", "Northing") # define x y
proj4string(data.sp) <- CRS("+init=epsg:23030") # define projection system of our data ## EPSG CODE 23030 ED50 / UTM zone 30N
summary(data.sp)
detach(package:rgdal)
data.sp2 <- spTransform(data.sp, CRS("+init=epsg:4326")) ## change projection in WGS84 lat lon
data.nz$Lon <- coordinates(data.sp2)[, "Easting"]
data.nz$Lat <- coordinates(data.sp2)[, "Northing"]
## ## plot on world map library(rworldmap) newmap <- getMap(resolution = 'coarse')
## # different resolutions available plot(newmap)
## points(data.sp2,cex=0.2,col='red')
rm(data.sp, data.sp2)
###################### ECOREGION merge greco to have no ecoregion with low number of observation
data.nz <- merge(data.nz, data.frame(plot = data.plot$plid, data.plot[, 11:12]),
by = "plot")
table(data.nz$Broad)
###################### PERCENT DEAD
perc.dead <- tapply(data.nz[["dead"]], INDEX = data.nz[["plot"]], FUN = function.perc.dead)
# ## VARIABLE TO SELECT PLOT WITH NOT BIG DISTURBANCE KEEP OFTHER VARIABLES IF
# AVAILABLE (disturbance record)
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109
data.nz <- merge(data.nz, data.frame(plot = names(perc.dead), perc.dead = perc.dead),
by = "plot", sort = FALSE)
########################################################### PLOT SELECTION FOR THE ANALYSIS
data.nz <- merge(data.nz, data.plot[, c(1, 8:10)], by = "plid")
colnames(data.nz)[colnames(data.nz) %in% c("mat", "map")] <- c("MAT", "MAP")
# colnames(data.nz)[names(data.nz) =='eco_codemerged' ] <- c('ecocode')
vec.abio.var.names <- c("MAT", "MAP")
vec.basic.var <- c("tree.id", "sp", "spname", "plot", "ecocode", "D", "G", "dead",
"year", "htot", "Lon", "Lat", "perc.dead")
data.tree <- subset(data.nz, 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
data.BA.SP <- BA.SP.FUN(id.tree = as.vector(data.nz[["tree.id"]]), diam = as.vector(data.nz[["D"]]),
sp = as.vector(data.nz[["sp"]]), id.plot = as.vector(data.nz[["plot"]]), weights = as.vector(1/(pi *
(0.5 * data.nz[["D0"]]/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.nz[["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], MARGIN = 1, FUN = sum, na.rm = TRUE)
data.BA.SP$BATOT.COMPET <- BATOT.COMPET
### create data frame
names(data.BA.SP) <- c("tree.id", names(data.BA.SP)[-1])
data.BA.sp <- merge(data.frame(tree.id = dataIFN.spain[["tree.id"]], ecocode = dataIFN.spain[["ecocode"]]),
data.BA.SP, by = "tree.id", sort = FALSE)
## test
if (sum(!data.BA.sp[["tree.id"]] == data.tree[["tree.id"]]) > 0) stop("competition index not in the same order than data.tree")
## save everything as a list
list.nz <- list(data.tree = data.tree, data.BA.SP = data.BA.sp, data.traits = data.trait)
save(list.nz, file = "./data/process/list.nz.Rdata")