Skip to content
GitLab
Select Git revision
  • master default
  • BRIEFCOM
  • ctfs
  • test.dplyr
  • package
Find file BlameHistoryPermalink
merge.data.SWEDEN.R 7.64 KiB
1 
### MERGE sweden DATA
2 
rm(list = ls()); source("./R/format.function.R"); library(reshape);
3
4 
######################### READ DATA read individuals tree data
5 
#data.swe1 <- read.csv("./data/raw/DataSweden/Swe_NFI_1.csv",header=T,stringsAsFactors=F)
6 
#data.swe2 <- read.csv("./data/raw/DataSweden/Swe_NFI_2a.csv",header=T,stringsAsFactors=F)
7 
#data.swe3 <- read.csv("./data/raw/DataSweden/Swe_NFI_3.csv",header=T,stringsAsFactors=F)
8 
data.swe <- read.table("./data/raw/DataSweden/Swe_NFI_all.txt",header=T,stringsAsFactors=F,sep="\t")
9 
head(data.swe)
10 
data.swe$tree.id <- apply(cbind(data.swe[["TreeNr"]],data.swe[["TractNr"]],data.swe[["PlotNr"]],data.swe[["Year"]]),1,paste,collapse="_")
11 
data.swe$plot.id <- apply(cbind(data.swe[["TractNr"]],data.swe[["PlotNr"]],data.swe[["Year"]]),1,paste,collapse="_")
12 
#data.swe <- data.swe[order(data.swe$tree.id),]
13
14 
## Species names are in the xlsx files if required (we already have sp codes)
15 
#dim(data.swe)
16 
#table(table(data.swe$TreeID))
17 
#table(table(paste(data.swe$TractNr,data.swe$PlotNr,data.swe$TreeNr,data.swe$Year)))
18 
#data.swe <- rbind(data.swe1, data.swe2, data.swe3);
19 
#rm(data.swe1, data.swe2, data.swe3)
20 
#data.swe <- data.swe[order(data.swe$TreeID,data.swe$PlotInvent),] ## Shows the TreeID = "" first
21 
#data.swe <- data.swe[order(paste(data.swe$TractNr,data.swe$PlotNr,data.swe$TreeNr)),] ## Shows the TreeID = "" first
22
23 
## Format to desired form
24 
data.swe$PlotInventID <- NULL; data.swe$TreeID <- NULL;
25 
data.swe2 <- data.swe[rep(1:nrow(data.swe),each=2),]
26 
data.swe2$dbh1 <- as.vector(rbind(data.swe$dia_t1,data.swe$dia_t2))
27 
data.swe2$dbh2 <- as.vector(rbind(data.swe$dia_t2,data.swe$dia_t3))
28 
data.swe2$dia_t1 <- data.swe2$dia_t2 <- data.swe2$dia_t3 <- NULL
29 
data.swe2$vol1 <- as.vector(rbind(data.swe$vol_t1,data.swe$vol_t2))
30 
data.swe2$vol2 <- as.vector(rbind(data.swe$vol_t2,data.swe$vol_t3))
31 
data.swe2$vol_t1 <- data.swe2$vol_t2 <- data.swe2$vol_t3 <- NULL
32 
data.swe2$dryw1 <- as.vector(rbind(data.swe$DryW_t1,data.swe$DryW_t2))
33 
data.swe2$dryw2 <- as.vector(rbind(data.swe$DryW_t2,data.swe$DryW_t3))
34 
data.swe2$DryW_t1 <- data.swe2$DryW_t2 <- data.swe2$DryW_t3 <- NULL
35
36 
data.swe2$Diameter <- data.swe2$Volume <- data.swe2$BrhAge <- NULL
37 
data.swe <- data.swe2
38 
rm(data.swe2)
39
40 
###################################### MASSAGE TRAIT DATA Compute maximum height per species plus sd from observed
41 
###################################### height to add variables to the traits data base Because we have two heights,
42 
###################################### then take the max of the two heights and then bootstrap
43 
## Obtain heights in dataset - I assume that they use the volume of a cylinder V = pi*r^2*h
44 
data.swe$ht1 <- data.swe$vol1/(pi*(0.5*data.swe$dbh1/1000)^2)
45 
data.swe$ht2 <- data.swe$vol2/(pi*(0.5*data.swe$dbh2/1000)^2)
46
47 
res.quant.boot <- t(sapply(levels(factor(data.swe[["TreeSpecies"]])), FUN = f.quantile.boot,
48 
    R = 1000, x = log10(apply(data.swe[, c("ht1", "ht2")], 1, max, na.rm = T)),
49 
    fac = factor(data.swe[["TreeSpecies"]])))
50
51 
## create data base
52 
data.max.height <- data.frame(code = rownames(res.quant.boot), Max.height.mean = res.quant.boot[,
53 
    1], Max.height.sd = res.quant.boot[, 2], Max.height.nobs = res.quant.boot[, 3], stringsAsFactors =FALSE)
54 
rm(res.quant.boot)
55 
write.csv(data.max.height,file='./data/process/data.max.height.swe.csv')
56
57 
########################################## FORMAT INDIVIDUAL TREE DATA
58 
## STOP HERE!!!
59 
data.swe$obs.id <- 1:nrow(data.swe)
60 
data.swe$census.id <- rep(1:2,nrow(data.swe)/2) ## Using first census in the row
61 
data.swe$G <- 10*(data.swe$dbh2-data.swe$dbh1)/data.swe$Interval ## diameter growth in mm per year
62 
data.swe$G[which(data.swe$InitDBH == 0 | data.swe$FinalDBH == -999)] <- NA
63 
data.swe$year <- data.swe$Interval ## number of year between measuremen
64 
data.swe$D <- data.swe[["InitDBH"]]; data.swe$D[data.swe$D == 0] <- NA ;## diameter in cm
65 
data.swe$dead <- as.numeric(data.swe$FinalDBH > 0) ## dummy variable for dead tree 0 alive 1 dead
66 
data.swe$sp <- as.character(data.swe[["Species"]]) ## species code
67 
data.swe$plot <- (data.swe[["PLOT_ID"]]) ## plot code
68 
data.swe$htot <- rep(NA,length(data.swe[["Species"]])) ## height of tree in m - MISSING
69 
data.swe$tree.id <- gsub("_",".",data.swe$PLOTTREE); ## tree unique id
70 
data.swe$sp.name <- NA;
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
for(i in 1:length(unique(data.swe$sp))) { v <- species.clean$SPCD data.swe$sp.name[which(data.swe$sp == unique(data.swe$sp)[i])] <- species.clean$COMMON_NAME[which(v == unique(data.swe$sp)[i])] } ###################### ## ECOREGION ################### ## merge greco to have no ecoregion with low number of observation greco <- read.csv(file = "./data/raw/DataCanada/EcoregionCodes.csv", header = T, sep = "\t") table(data.swe$Ecocode) ## Some ecoregions still have small # of individuals, so either drop off for analysis later on or wait for Quebec data to come in # # library(RColorBrewer); mycols <- brewer.pal(10,"Set3"); # ecoreg <- unclass(data.swe$eco_code); # plot(data.swe[["CX"]][order(ecoreg)],data.swe[["CY"]][order(ecoreg)],pty=".",cex=.2, col = rep(mycols,as.vector(table(ecoreg)))); # legend("bottomright", col = mycols, legend = levels(data.swe$eco_code), pch = rep(19,length(levels(ecoreg))),cex=2) # points(data.swe[["CX"]][ecoreg == 9],data.swe[["CY"]][ecoreg == 9],pty=".",cex=.2, col = "black"); ## Highlight the region with 55 sites # ## PA1219 looks to be similar to PA1209; merge them together # data.swe$eco_codemerged <- combine_factor(data.swe$eco_code, c(1:8,6,9)) ###################### ## PERCENT DEAD ################### ## variable percent dead/cannot do with since dead variable is missing ## compute numer of dead per plot to remove plot with disturbance perc.dead <- tapply(data.swe[["dead"]],INDEX=data.swe[["plot"]],FUN=function.perc.dead) # ## VARIABLE TO SELECT PLOT WITH NOT BIG DISTURBANCE KEEP OFTHER VARIABLES IF AVAILABLE (disturbance record) data.swe <- merge(data.swe,data.frame(plot=names(perc.dead),perc.dead=perc.dead), by = "plot", sort=FALSE) ########################################################### ### PLOT SELECTION FOR THE ANALYSIS ################### ## Remove data with dead == 1 table(data.swe$dead) ## Nothing to remove colnames(data.swe)[c(3,1,11,13)] <- c("sp","plot","w","ecocode") vec.abio.var.names <- c("MAT","MAP") vec.basic.var <- c("tree.id","sp","sp.name","plot","ecocode","D","G","dead","year","htot","Lon","Lat","perc.dead") data.tree <- subset(data.swe,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.swe[["tree.id"]]), diam=as.vector(data.swe[["D"]]), sp=as.vector(data.swe[["sp"]]), id.plot=as.vector(data.swe[["plot"]]), weights=1/(10000*data.swe[["SubPlot_Size"]]), 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.swe[["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],1,sum,na.rm=TRUE) data.BA.SP$BATOT.COMPET <- BATOT.COMPET; rm(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=data.swe[["tree.id"]],ecocode=data.swe[["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.swe <- list(data.tree=data.tree,data.BA.SP=data.BA.sp,data.traits=data.traits) save(list.spain,file="./data/process/list.swe.Rdata")