Skip to content
GitLab
Commit b5b1ce84 authored by Georges Kunstler's avatar Georges Kunstler
Browse files

new competition index function

parent 6e6c9f83
master BRIEFCOM ctfs package test.dplyr
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 126 additions and 41 deletions
+126 -41
R/format.function.R
+ 74
− 20
View file @ b5b1ce84
...@@ -3,6 +3,20 @@ ...@@ -3,6 +3,20 @@
################################################### ###################################################
##### FUNCTION TO FORMAT DATA FOR THE WORKSHOP ANALYSIS ##### 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.. ##' .. Compute the basal area of competitor in a plot..
...@@ -19,6 +33,7 @@ BA.fun <- function(diam,weights){ ...@@ -19,6 +33,7 @@ BA.fun <- function(diam,weights){
#### ####
##' .. function to compute the basal area of neighborood tree in plots .. ##' .. function to compute the basal area of neighborood tree in plots ..
##' ##'
...@@ -32,34 +47,73 @@ BA.fun <- function(diam,weights){ ...@@ -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) ##' @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) ##' @return data frame with tree.id and one column per species with basal area of the species (without the target tree)
##' @author Kunstler ##' @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 # 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))) 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") 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) BASP <- tapply(BA.fun(diam,weights),INDEX=list(id.plot,sp),FUN=sum,na.rm=T)
print(dim(BASP)) print(dim(BASP))
DATA.BASP <- data.frame(id.plot= rownames(BASP),BASP) DATA.BASP <- data.table(id.plot= rownames(BASP),BASP)
names( DATA.BASP) <- c("id.plot",colnames(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 #### MERGE with indivudal tree
data.indiv <- data.frame(id.tree,sp,id.plot,diam,weights) ## use library(data.table)
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)]))
if(!is.na(weight.full.plot)){ if(!is.na(weight.full.plot)){
t.arbre.BASP[t(outer(data.merge$sp,names(data.merge[,-(1:5)]),FUN="=="))] <- data.indiv <- data.table(id.tree=id.tree,sp=sp,
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))) id.plot=id.plot,diam=diam,
}else{ BA.indiv=BA.fun(diam,rep(weight.full.plot,length(diam))))
t.arbre.BASP[t(outer(data.merge$sp,names(data.merge[,-(1:5)]),FUN="=="))] <- setkeyv(data.indiv,"id.plot")
t.arbre.BASP[t(outer(data.merge$sp,names(data.merge[,-(1:5)]),FUN="=="))] -BA.fun(data.merge$diam,data.merge$weights) print("second table created")
} data.merge <- merge(data.indiv,DATA.BASP)
data.res <- data.frame(data.merge$id.tree,t(t.arbre.BASP)) print("merge done")
names(data.res) <- c("id.tree",colnames(BASP)) ## mat.arbre.BASP <- as.matrix(subset(data.merge,select=sp.name))
return( data.res) ## 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 #### function with X Y coordinates based on a neighborhood of radius R
......
merge.data.PARACOU.R
+ 8
− 12
View file @ b5b1ce84
...@@ -21,23 +21,19 @@ data.paracou$treeid <- apply(data.paracou[,1:4],1,paste,collapse="."); ## Create ...@@ -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))] 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!!! ## 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) 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 ##draw map of each plots
plot(data.paracou[["x"]],data.paracou[["y"]]) plot(data.paracou[["x"]],data.paracou[["y"]])
......
merge.data.US.R
+ 44
− 9
View file @ b5b1ce84
...@@ -28,12 +28,13 @@ data.us$year <- data.us$Interval ## number of year between measuremen ...@@ -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$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$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$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$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; data.us$sp.name <- NA;
v <- species.clean$SPCD; for(i in 1:length(unique(data.us$sp))) { ## 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])] } ## 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 ## ECOREGION
...@@ -67,18 +68,52 @@ data.us <- merge(data.us,data.frame(plot=names(perc.dead),perc.dead=perc.dead), ...@@ -67,18 +68,52 @@ data.us <- merge(data.us,data.frame(plot=names(perc.dead),perc.dead=perc.dead),
################### ###################
## Remove data with dead == 1 ## Remove data with dead == 1
table(data.us$dead) 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.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") 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 ## 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"]]), ecoregion.name <- unique(data.us[["eco_codemerged"]])
sp=as.vector(data.us[["sp"]]), id.plot=as.vector(data.us[["plot"]]), subplot.1 <- unique(data.us[eco_codemerged %in% ecoregion.name[1:5],subplot])
weights=1/(10000*data.us[["PlotSize"]]), weight.full.plot=NA) 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 ## 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 data.BA.SP[is.na(data.BA.SP)] <- 0; data.BA.SP[,-1][data.BA.SP[,-1]<0] <- 0
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment