merge nz data R script

merge.data.NZ-fhv1.R

+### 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 <- rep(NA,length(data.nz[["sp"]])) ## dummy variable for dead tree 0 alive 1 dead - MISSING
+data.nz$sp <- data.nz$sp
+data.nz$plot <- data.nz$plid 
+data.nz$htot <- rep(NA,length(data.nz[["sp"]])) ## I ASSUME THE "height.m" IN data.trait IS A MAX HEIGHT, WHICH IS NOT WHAT WE WANT HERE?
+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 coordinates system of Easting Northing to be in lat long WGS84 - DON'T KNOW THE EPSG CODE HERE
+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)
+
+############################
+## DON'T KNOW WHICH VARIABLE TO USE FOR ECOREGION =(
+greco <- greco[,c("tree.id","BIOME","eco_code")]
+
+#######################
+## variable percent dead - MISSING
+data.nz$perc.dead <- NA
+
+###########################################################
+### 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","PP")
+#colnames(data.nz)[names(data.nz) =="eco_codemerged"  ] <- c("ecocode")
+vec.abio.var.names <-  c("MAT","PP")
+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"]])^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")
+