diff --git a/R/format.function.R b/R/format.function.R index 90f7f4c3bec476ceea91ad8603ef15ceac20c03b..c87fcac3729d25dc42a8fd8167f4a8bd1d5408cb 100644 --- a/R/format.function.R +++ b/R/format.function.R @@ -3,6 +3,20 @@ ################################################### ##### FUNCTION TO FORMAT DATA FOR THE WORKSHOP ANALYSIS + +###### +###### +## FUNCTION TO PLOT MAP OF TREE +fun.circles.plot <- function(plot.select,x,y,plot,D,inches){ +x.t <- x[plot==plot.select] +y.t <- y[plot==plot.select] +D.t <- D[plot==plot.select] +D.t[is.na(D.t)] <- 0 +symbols(x.t,y.t,circles=D.t ,main=plot.select,inches=inches) +} + + + ######################### ## ##' .. Compute the basal area of competitor in a plot.. @@ -19,6 +33,7 @@ BA.fun <- function(diam,weights){ + #### ##' .. function to compute the basal area of neighborood tree in plots .. ##' @@ -32,34 +47,73 @@ BA.fun <- function(diam,weights){ ##' @param weights.full.plot weights for the whole plot to compute basal area in cm^2/m^2 or if NA use weights of the individuals (for simple plots) ##' @return data frame with tree.id and one column per species with basal area of the species (without the target tree) ##' @author Kunstler -BA.SP.FUN <- function(id.tree,diam,sp,id.plot,weights,weight.full.plot){ +BA.SP.FUN <- function(id.tree,diam,sp,id.plot,weights,weight.full.plot,name,num){ +require(data.table) +id.plot <- as.character(id.plot) +id.rtee <- as.character(id.tree) # compute BA tot per species per plot if(!(length(id.tree)==length(diam) & length(id.tree)==length(sp) & length(id.tree)==length(id.plot) & length(id.tree)==length(weights))) stop("length of id.tree diam,sp id.plot & weights need to be the same") -sp <- as.character(sp) +sp <- paste("sp",sp,sep=".") BASP <- tapply(BA.fun(diam,weights),INDEX=list(id.plot,sp),FUN=sum,na.rm=T) print(dim(BASP)) -DATA.BASP <- data.frame(id.plot= rownames(BASP),BASP) -names( DATA.BASP) <- c("id.plot",colnames(BASP)) - +DATA.BASP <- data.table(id.plot= rownames(BASP),BASP) +setnames( DATA.BASP,old=1:ncol(DATA.BASP), c("id.plot",colnames(BASP))) +setkeyv(DATA.BASP,c("id.plot")) +sp.name <- colnames(BASP) +rm(BASP) +print("first table created") #### MERGE with indivudal tree -data.indiv <- data.frame(id.tree,sp,id.plot,diam,weights) -data.merge <- merge(data.indiv,DATA.BASP,by="id.plot", sort = FALSE) - -## substracte the BA of target species in the column of the good species (with transpose of the matrix ) -t.arbre.BASP <- t(as.matrix(data.merge[,-(1:5)])) +## use library(data.table) if(!is.na(weight.full.plot)){ - t.arbre.BASP[t(outer(data.merge$sp,names(data.merge[,-(1:5)]),FUN="=="))] <- - t.arbre.BASP[t(outer(data.merge$sp,names(data.merge[,-(1:5)]),FUN="=="))] -BA.fun(data.merge$diam,rep(weight.full.plot,length(diam))) - }else{ - t.arbre.BASP[t(outer(data.merge$sp,names(data.merge[,-(1:5)]),FUN="=="))] <- - t.arbre.BASP[t(outer(data.merge$sp,names(data.merge[,-(1:5)]),FUN="=="))] -BA.fun(data.merge$diam,data.merge$weights) - } -data.res <- data.frame(data.merge$id.tree,t(t.arbre.BASP)) -names(data.res) <- c("id.tree",colnames(BASP)) -return( data.res) -} + data.indiv <- data.table(id.tree=id.tree,sp=sp, + id.plot=id.plot,diam=diam, + BA.indiv=BA.fun(diam,rep(weight.full.plot,length(diam)))) + setkeyv(data.indiv,"id.plot") + print("second table created") + data.merge <- merge(data.indiv,DATA.BASP) + print("merge done") + ## mat.arbre.BASP <- as.matrix(subset(data.merge,select=sp.name)) + ## rownames(mat.arbre.BASP) <- data.merge[["id.tree"]] + ## mat.arbre.BASP[cbind(1:nrow(data.merge),match(data.merge[["sp"]],sp.name))] <- mat.arbre.BASP[cbind(1:nrow(data.merge),match(data.merge[["sp"]],sp.name))] - data.merge[["BA.indiv"]] + + for (i in (sp.name)){ + eval(parse(text=paste("data.merge[sp==\'",i,"\',",i,":=",i,"-BA.indiv]",sep=""))) + } + +}else{ + data.indiv <- data.table(id.tree=id.tree,sp=sp,id.plot=id.plot, + diam=diam,weights=weights, + BA.indiv=BA.fun(diam,weights)) + setkeyv(data.indiv,"id.plot") + print("second table created") + data.merge <- merge(data.indiv,DATA.BASP) + print("merge done") + ## mat.arbre.BASP <- as.matrix(subset(data.merge,select=sp.name)) + ## rownames(mat.arbre.BASP) <- data.merge[["id.tree"]] + ## mat.arbre.BASP[cbind(1:nrow(data.merge),match(data.merge[["sp"]],sp.name))] <- mat.arbre.BASP[cbind(1:nrow(data.merge),match(data.merge[["sp"]],sp.name))] - data.merge[["BA.indiv"]] + + for (i in (sp.name)){ + eval(parse(text=paste("data.merge[sp==\'",i,"\',",i,":=",i,"-BA.indiv]",sep=""))) + } +} +## print("replacment done") +## data.res <- data.table(id.tree=data.merge$id.tree,t(t.arbre.BASP)) +## print("data.table created") +## setnames(data.res,old=1:ncol(data.res),c("id.tree",sp.name)) + +print("replacment done") +data.merge[,BA.indiv:=NULL] +print("first column removed") +data.merge[,sp:=NULL] +data.merge[,diam:=NULL] +data.merge[,id.plot:=NULL] +data.merge[,weights:=NULL] +print("columns removed") +## save(data.merge,file=paste("./data/process/data.BA.SP",name,num,".Rdata",sep="")) + return( data.merge) +} #### function with X Y coordinates based on a neighborhood of radius R diff --git a/merge.data.PARACOU.R b/merge.data.PARACOU.R index fa50b4789b5fc93936a6c91641a5e4cc1b8c9814..45f690a7cd9fe6ad8cba475279435ea236f5e54b 100644 --- a/merge.data.PARACOU.R +++ b/merge.data.PARACOU.R @@ -21,23 +21,19 @@ data.paracou$treeid <- apply(data.paracou[,1:4],1,paste,collapse="."); ## Create data.paracou <- data.paracou[,c(ncol(data.paracou),1:(ncol(data.paracou)-1))] ## REMOVE ALL TREES WITH X OR Y >250 m or NA ASK GHISLAIN!!! - data.paracou <- subset(data.paracou,subset=(!is.na(data.paracou[["x"]])) & data.paracou[["x"]]<251 & data.paracou[["y"]]<251) -fun.circles.plot <- function(plot.select,x,y,plot,D){ - -x.t <- x[plot==plot.select] -y.t <- y[plot==plot.select] -D.t <- D[plot==plot.select] -D.t[is.na(D.t)] <- 0 -symbols(x.t,y.t,circles=D.t ,main=plot.select) -} -fun.circles.plot(unique(data.paracou[["plot"]])[1],data.paracou[['x']],data.paracou[['y']],data.paracou[["plot"]],data.paracou[["circum2009"]]) +## plot each +pdf("./figs/plots.paracou.pdf") +lapply(unique(data.paracou[["plot"]]),FUN=fun.circles.plot,data.paracou[['x']],data.paracou[['y']],data.paracou[["plot"]],data.paracou[["circum2009"]],inches=0.2) +dev.off() +table(data.paracou[["plot"]],data.paracou[["subplot"]]) +#################### +#### PLOT 16 17 18 ARE STRANGE ASK GHSILAIN +##################### -sum(is.na(data.paracou[["x"]])) -sum(is.na(data.paracou[["y"]])) ##draw map of each plots plot(data.paracou[["x"]],data.paracou[["y"]]) diff --git a/merge.data.US.R b/merge.data.US.R index 5acd84e22694e440027f440de0bb48d89c9d3e4d..701fea11ff89e8bab7e19c690ef734afa7607875 100644 --- a/merge.data.US.R +++ b/merge.data.US.R @@ -28,12 +28,13 @@ data.us$year <- data.us$Interval ## number of year between measuremen data.us$D <- data.us[["InitDbh"]]; data.us$D[data.us$D == 0] <- NA ;## diameter in cm data.us$dead <- as.numeric(data.us$FinalDbh > 0) ## dummy variable for dead tree 0 alive 1 dead data.us$sp <- as.character(data.us[["Species"]]) ## species code -data.us$plot <- (data.us[["PlotID"]]) ## plot code +data.us$plot <- as.character(data.us[["PlotID"]]) ## plot code +data.us$subplot <- paste(as.character(data.us[["PlotID"]]),as.character(data.us[["SubplotNumber"]]),sep=".") ## plot code data.us$htot <- rep(NA,length(data.us[["Species"]])) ## height of tree in m - MISSING -data.us$tree.id <- data.us$TreeID; ## tree unique id +data.us$tree.id <- as.character(data.us$TreeID); ## tree unique id data.us$sp.name <- NA; -v <- species.clean$SPCD; for(i in 1:length(unique(data.us$sp))) { - data.us$sp.name[which(data.us$sp == unique(data.us$sp)[i])] <- species.clean$COMMON_NAME[which(v == unique(data.us$sp)[i])] } +## v <- species.clean$SPCD; for(i in 1:length(unique(data.us$sp))) { +## data.us$sp.name[which(data.us$sp == unique(data.us$sp)[i])] <- species.clean$COMMON_NAME[which(v == unique(data.us$sp)[i])] } ###################### ## ECOREGION @@ -67,18 +68,52 @@ data.us <- merge(data.us,data.frame(plot=names(perc.dead),perc.dead=perc.dead), ################### ## Remove data with dead == 1 table(data.us$dead) -data.us <- data.us[data.us$dead == 1,] +## data.us <- data.us[data.us$dead == 1,] vec.abio.var.names <- c("MAT","MAP") vec.basic.var <- c("tree.id","sp","sp.name","plot","eco_codemerged","D","G","dead","year","htot","Lon","Lat","perc.dead") -data.tree <- subset(data.us,select=c(vec.basic.var,vec.abio.var.names)) +## data.tree <- subset(data.us,select=c(vec.basic.var,vec.abio.var.names)) + +## remove everything from memory not ned before computation +rm(greco,perc.dead,species.clean,tab.small.div,sel.small.div) +library(data.table) +data.us <- data.table(data.us) +data.us[,Species:=NULL] +data.us[,FinalDbh:=NULL] +data.us[,InitDbh:=NULL] +data.us[,IndWeight:=NULL] +data.us[,SubplotNumber:=NULL] +data.us[,IntervalYears:=NULL] + +gc() ############################################## ## 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.us[["tree.id"]]), diam=as.vector(data.us[["D"]]), - sp=as.vector(data.us[["sp"]]), id.plot=as.vector(data.us[["plot"]]), - weights=1/(10000*data.us[["PlotSize"]]), weight.full.plot=NA) +ecoregion.name <- unique(data.us[["eco_codemerged"]]) +subplot.1 <- unique(data.us[eco_codemerged %in% ecoregion.name[1:5],subplot]) +subplot.2 <- unique(data.us[eco_codemerged %in% ecoregion.name[6:11],subplot]) + +system.time(data.BA.SP.1 <- BA.SP.FUN(id.tree=as.vector(data.us[subplot %in% subplot.1,tree.id]), + diam=as.vector(data.us[subplot %in% subplot.1,D]), + sp=as.vector(data.us[subplot %in% subplot.1,sp]), + id.plot=as.vector(data.us[subplot %in% subplot.1,subplot]), + weights=1/(10000*data.us[subplot %in% subplot.1,PlotSize]), weight.full.plot=NA, + name="France",num=1)) + +system.time(data.BA.SP.2 <- BA.SP.FUN(id.tree=as.vector(data.us[subplot %in% subplot.2,tree.id]), + diam=as.vector(data.us[subplot %in% subplot.2,D]), + sp=as.vector(data.us[subplot %in% subplot.2,sp]), + id.plot=as.vector(data.us[subplot %in% subplot.2,subplot]), + weights=1/(10000*data.us[subplot %in% subplot.2,PlotSize]), weight.full.plot=NA, + name="France",num=2)) +save(data.BA.SP.1,file="./data/process/data.BA.SP.US.1.Rdata") +save(data.BA.SP.2,file="./data/process/data.BA.SP.US.2.Rdata") + +Possibility to merge that ? +# +data.BA.SP <- rbind(data.BA.SP.1,data.BA.SP.2) + ## 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