-
Daniel Falster authored0ed8f24b
### MERGE NVS DATA
rm(list = ls())
source("format.fun.R")
library(reshape)
dir.create("../../output/formatted/NVS", recursive=TRUE,showWarnings=FALSE)
###################################### 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.trait <- read.csv("../../data/raw/NVS/nz_traits_130801.csv", , header = TRUE,
stringsAsFactors = FALSE)
data.trait <- data.trait[, -1]
colnames(data.trait)[1] <- "sp"
data.trait$Latin_name <- data.trait$sp
data.trait$Leaf.N.mean <- data.trait$leafn*10; data.trait$leafn <- NULL ## conversion from % to mg/g
data.trait$Leaf.N.sd <- NA ## conversion from % to mg/g
data.trait$Seed.mass.mean <- data.trait$seed.mg; data.trait$seed.mg <- NULL
data.trait$Seed.mass.sd <- NA
data.trait$SLA.mean <- 1/data.trait$lma.gm2; data.trait$SLA.mean <- data.trait$SLA.mean*1000;data.trait$lma.gm2 <- NULL ## conversion of g/m2 to mm2/g1
data.trait$SLA.sd <- NA
data.trait$Wood.density.mean <- data.trait$wood; data.trait$wood <- NULL
data.trait$Wood.density.sd <- NA
data.trait$Max.height.mean <- log10(data.trait$height.m); data.trait$height.m <- NULL
data.trait$Max.height.sd <- NA
write.csv(data.trait,file='../../output/formatted/NVS/trait.csv')
######################### READ DATA read individuals tree and environmental data
data.nz <- read.csv("../../data/raw/NVS/nz_treedata_growth_130801.csv", header = TRUE,
stringsAsFactors = FALSE, skip = 9)
data.nz <- data.nz[, -1]
data.nz$plid <- gsub("__", "_", data.nz$plid)
data.nz$plid <- gsub("_", ".", data.nz$plid) ## Replace all underscores with a single dot
data.plot <- read.csv("../../data/raw/NVS/nz_plotinfo_130801.csv", , header = TRUE,
stringsAsFactors = FALSE)
data.plot <- data.plot[, -1]
data.plot$plid <- gsub("__", "_", data.plot$plid)
data.plot$plid <- gsub("_", ".", data.plot$plid) ## Replace all underscores with a single dot
data.nz <- merge(data.nz, data.frame(plid = data.plot$plid, Easting = data.plot$Easting,
Northing = data.plot$Northing, Locality = data.plot$Locality, MAP = data.plot$map, MAT = data.plot$mat, Broadclass = data.plot$BroadClass,
Specificclass = data.plot$SpecificClass, stringsAsFactors = F), sort = T, by = "plid")
rm(data.plot)
################################################################ 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$sp.name <- data.nz$sp ## no latin name so use that as not needed to get traits
data.nz$year <- (data.nz[["t1"]] - data.nz[["t0"]]) ## number of year between measurement
data.nz$D <- data.nz[["D0"]] ## diameter in cm
data.nz$dead <- as.numeric(is.na(data.nz[["D1"]])) ## dummy variable for dead tree 0 alive 1 dead
data.nz$plot <- data.nz$plid; data.nz$plid <- NULL
data.nz$cluster <- rep(NA, nrow(data.nz)) ## no plot cluster
data.nz$htot <- rep(NA, nrow(data.nz)) ## Max height is already available so have as missing
data.nz$tree.id <- data.nz$tag
head(subset(data.nz,subset=data.nz$tag %in% names(table(data.nz$tag))[table(data.nz$tag)==2])) ### 3 individuals with two observation with strange value check with Daniel and David add an issue
data.nz$obs.id <- as.character(1:nrow(data.nz))
data.nz$census <- rep(1,nrow(data.nz)) # only one census NEED TO CHECK WITH DANIEL AND DAVID add issue
data.nz$weights <- 1/(20*20) ## need to check with david OK
###################### ECOREGION merge to have no ecoregion with low number of observation
table(data.nz$Broad)
## merging two ecoregion based on Daniel Laughlin advice.
data.nz$Broad[data.nz$Broad=="Wetland"] <- "BeechHumid"
data.nz$Broad[data.nz$Broad=="ConiferDry"] <- "BeechHumid"
data.nz$ecocode <- data.nz$Broad
########################################## CHANGE COORDINATE SYSTEM
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112
########################################## system of Easting Northing to be in lat long WGS84
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:27200") # define projection system of our data ## EPSG CODE 27200 NZGD49 / New Zealand Map Grid
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"]
rm(data.sp, data.sp2)
## ## plot on world map
# pdf("./figs/map.NZ.pdf")
# library(RColorBrewer)
# library(rworldmap)
# newmap <- getMap(resolution = 'coarse')
# ## # different resolutions available
# plot(newmap,xlim=c(166,177),ylim=c(-34,-47))
# ecoreg <- factor(data.nz$ecocode); levels(ecoreg) <- mycols <- brewer.pal(length(levels(ecoreg)), "Set1")
# points(data.nz[['Lon']],data.nz[['Lat']],pty='.',cex=.2,col = as.character(ecoreg))
# legend('bottomleft', col = mycols, legend = names(table(data.nz$ecocode)), pch = rep(19,length(levels(ecoreg))))
# points(data.sp2,cex=0.2,col='red')
# dev.off()
###################### PERCENT DEAD
perc.dead <- tapply(data.nz[["dead"]], INDEX = data.nz[["plot"]], FUN = function.perc.dead2)
data.nz <- merge(data.nz, data.frame(plot = names(perc.dead), perc.dead = perc.dead),
by = "plot", sort = FALSE)
########################################################### PLOT SELECTION FOR THE ANALYSIS
#####################
###### REMOVE DEAD TREE at first census
data.nz <- subset(data.nz,subset=!is.na(data.nz[["D"]]))
vec.abio.var.names <- c("MAT", "MAP")
vec.basic.var <- c("obs.id","tree.id", "sp", "sp.name","cluster","plot", "ecocode", "D", "G", "dead",
"year", "htot", "Lon", "Lat", "perc.dead","weights","census")
data.tree <- subset(data.nz, select = c(vec.basic.var, vec.abio.var.names))
write.csv(data.tree,file="../../output/formatted/NVS/tree.csv")