# Vignette BIOMASS

# 1 Load BIOMASS and datasets

Install BIOMASS (to be done once)

    install.packages("BIOMASS")

Load the package

    library(BIOMASS)
    require(knitr) # To build tables in this document

    ## Loading required package: knitr

Load the two datasets stored in the package

    data(KarnatakaForest)
    str(KarnatakaForest)

    ## 'data.frame':    61965 obs. of  8 variables:
    ##  $ plotId : Factor w/ 96 levels "BSP1","BSP10",..: 1 1 1 1 1 1 1 1 1 1 ...
    ##  $ treeId : Factor w/ 65283 levels "BSP1_001","BSP1_002",..: 1 2 4 5 6 7 8 9 10 11 ...
    ##  $ family : Factor w/ 81 levels "Acanthaceae",..: 18 61 29 64 6 64 30 29 64 26 ...
    ##  $ genus  : Factor w/ 256 levels "Acacia","Acrocarpus",..: 233 256 21 134 252 134 9 21 134 82 ...
    ##  $ species: Factor w/ 320 levels "alata","albiflavescens",..: 28 204 163 32 16 32 161 163 32 189 ...
    ##  $ D      : num  3.5 3.82 16.87 4.14 5.41 ...
    ##  $ lat    : num  14.4 14.4 14.4 14.4 14.4 ...
    ##  $ long   : num  74.9 74.9 74.9 74.9 74.9 ...

    #
    data(NouraguesHD)
    str(NouraguesHD)

    ## 'data.frame':    1051 obs. of  7 variables:
    ##  $ plotId : chr  "Plot1" "Plot1" "Plot1" "Plot1" ...
    ##  $ genus  : chr  "indet" "Qualea" "Dicorynia" "Protium" ...
    ##  $ species: chr  "indet" "rosea" "guianensis" "cf_guianense" ...
    ##  $ D      : num  11.5 11.6 83.9 15 36.8 13.5 17.8 17.8 15.9 17.8 ...
    ##  $ H      : num  12 16 40 18 27 20 24 21 22 24 ...
    ##  $ lat    : num  4.07 4.07 4.07 4.07 4.07 ...
    ##  $ long   : num  -52.7 -52.7 -52.7 -52.7 -52.7 ...

Select 10 plots for illustrative purpose

    selecPlot<-KarnatakaForest$plotId%in%c("BSP2","BSP12","BSP14","BSP26","BSP28","BSP30","BSP34","BSP44","BSP63","BSP65")
    KarnatakaForestsub<-droplevels(KarnatakaForest[selecPlot,])

# 2 Retrieve wood density

## 2.1 Check and retrieve taxonomy

First, check for any typo in the taxonomy

    Taxo<-correctTaxo(genus=KarnatakaForestsub$genus,species=KarnatakaForestsub$species)

    ## [1] "Calling http://taxosaurus.org/retrieve/a55f49e4a33763f1c378558ff50cf02f"
    ## [1] "Calling http://taxosaurus.org/retrieve/c9515144813c820a0ab24e209267ca6b"
    ## [1] "Calling http://taxosaurus.org/retrieve/87a56ea8949a2511fc8d271e2cd863d2"
    ## [1] "Calling http://taxosaurus.org/retrieve/16f5614d03129a9345352dd934c24e9d"
    ## [1] "Calling http://taxosaurus.org/retrieve/173f238d8a6fc827604610272a33ab79"
    ## [1] "Calling http://taxosaurus.org/retrieve/2870f350c757980c530c021d219bd12a"
    ## [1] "Calling http://taxosaurus.org/retrieve/1c38fec8de19cdb4abf94410a490764a"

    KarnatakaForestsub$genusCorr<-Taxo$genusCorrected
    KarnatakaForestsub$speciesCorr<-Taxo$speciesCorrected

If needed, retrieve APG III families and orders from genus names

    APG<-getTaxonomy(KarnatakaForestsub$genusCorr, findOrder =T)
    KarnatakaForestsub$familyAPG<-APG$family
    KarnatakaForestsub$orderAPG<-APG$order

## 2.2 Wood density

Retrieve wood density using the plot level average if no genus level information is available

    dataWD<-getWoodDensity(genus=KarnatakaForestsub$genusCorr,
                 species=KarnatakaForestsub$speciesCorr,
                 stand=KarnatakaForestsub$plotId)

    ## The reference dataset contains 16467 wood density values
    ## Your taxonomic table contains 196 taxa

The same but using the family average and adding other wood density values as references (here invented for the example)

    LocalWoodDensity<-data.frame(genus=c("Ziziphus","Terminalia","Garcinia"),
                                 species=c("oenopolia","bellirica","indica"),
                                 wd=c(0.65,0.72,0.65))

    dataWD<-getWoodDensity(genus=KarnatakaForestsub$genusCorr,
                 species=KarnatakaForestsub$speciesCorr,
                 family=KarnatakaForestsub$familyAPG,
                 stand=KarnatakaForestsub$plotID,
                 addWoodDensityData=LocalWoodDensity)

    ## The reference dataset contains 16470 wood density values
    ## Your taxonomic table contains 196 taxa

Below the number of wood density value estimated at the species, genus and plot level:

    # At species level
    sum(dataWD$levelWD=="species")

    ## [1] 2304

    # At genus level
    sum(dataWD$levelWD=="genus")

    ## [1] 2768

    # At plot level
    sum(!dataWD$levelWD%in%c("genus","species"))

    ## [1] 370

# 3 Build height-diameter models

You may compare different models at once

    HDmodel <- modelHD(D=NouraguesHD$D,
                       H =NouraguesHD$H,
                       drawGraph=TRUE,
                       useWeight=TRUE)

    Which model would you like to select to model your data ?

    ##1 : Log 1 (blue)
    ##----- RSE = 4.1893 (RSElog = 0.2211)
    ##----- Average bias =  0.0042

    ##2 : Log 2 (green)
    ##----- RSE = 4.1017 (RSElog = 0.2194)
    ----- Average bias =  0.003

    3 : Log 3 (red)
    ----- RSE = 4.1038 (RSElog = 0.2195)
    ----- Average bias =  0.003

    4 : Weibull (orange)
    ----- RSE =  4.1716
    ----- Average bias =  0.0052

    5 : Michaelis - Menten (purple)
    ----- RSE =  4.1662
    ----- Average bias =  0.0148

    1: 3

Compute the local H-D model with the lowest RSE

    HDmodel<-modelHD(D=NouraguesHD$D,
                     H=NouraguesHD$H,
                     method="log2",
                     useWeight =TRUE)

Compute models specific to given stands

    HDmodelPerPlot <- by(NouraguesHD,NouraguesHD$plotId,
                         function(x) modelHD(D=x$D,H=x$H, method="weibull",useWeight =T),
                         simplify=FALSE)
    RSEmodels<-sapply(HDmodelPerPlot,function(x) x$RSE)
    Coeffmodels<-lapply(HDmodelPerPlot,function(x) x$coefficients)
    ResHD<-data.frame(Plot=names(unlist(RSEmodels)),
                      a=round(unlist(sapply(Coeffmodels,"[",1)),3),
                      b=round(unlist(sapply(Coeffmodels,"[",2)),3),
                      c=round(unlist(sapply(Coeffmodels,"[",3)),3),
                      RSE=round(unlist(RSEmodels),3))
    kable(ResHD, row.names = F)

Plot | a | b | c | RSE
--- | --- | --- | --- | ---
Plot1 | 43.673 | 29.629 | 0.987 | 4.503
Plot2 | 369.917 | 46732.555 | 0.367 | 3.865

# 4 Retrieve height data

Retrieve height data from a local Height-diameter model (Note that using a HD model built on French guianan trees for Indian trees is only for illustrative purpose here)

    dataHlocal<-retrieveH(D=KarnatakaForestsub$D,
                          model =HDmodel)

Retrieve height data from a Feldpaush et al. (2012) averaged model

    dataHfeld<-retrieveH(D=KarnatakaForestsub$D,
                         region ="SEAsia")

Retrieve height data from Chave et al. (2012) equation 6

    dataHchave<-retrieveH(D=KarnatakaForestsub$D,
                          coord=cbind(KarnatakaForestsub$long,KarnatakaForestsub$lat))

# 5 Estimate AGB

Organize data

    KarnatakaForestsub$WD=dataWD$meanWD
    KarnatakaForestsub$H=dataHlocal$H
    KarnatakaForestsub$Hfeld=dataHfeld$H

Compute AGB(Mg) per tree

    AGBtree<-computeAGB(D=KarnatakaForestsub$D,
                        WD=KarnatakaForestsub$WD,
                        H =KarnatakaForestsub$H)

Compute AGB(Mg) per plot

    AGBPlotList<-by(KarnatakaForestsub, KarnatakaForestsub$plotId,
                    function(x) computeAGB(D=x$D,WD=x$WD,H=x$H),
                    simplify=F)
    AGBplot<-sapply(AGBPlotList,sum)

Compute AGB(Mg) per tree without height information (Eq. 7 from Chave et al. (2014))

    AGBPlotListChave<-by(KarnatakaForestsub, KarnatakaForestsub$plotId,
                    function(x) computeAGB(D=x$D,WD=x$WD,coord =cbind(x$long, x$lat)),
                    simplify=F)
    AGBplotChave<-sapply(AGBPlotListChave,sum)

Compute AGB(Mg) per tree with Feldpausch et al. (2012) regional H-D model

    AGBPlotListFeld<-by(KarnatakaForestsub, KarnatakaForestsub$plotId,
                    function(x) computeAGB(D=x$D,WD=x$WD,H=x$Hfeld),
                    simplify=F)
    AGBplotFeld<-sapply(AGBPlotListFeld,sum)

# 6 Propagate AGB errors

Organize data

    KarnatakaForestsub$sdWD=dataWD$sdWD
    KarnatakaForestsub$HfeldRSE=dataHfeld$RSE

Propagate error for all tree at once using the local HD model constructed above (modelHD), i.e. non-independent allometric errors will be assigned to all trees at each iteration, independently of plots.

    resultMC<-AGBmonteCarlo(D=KarnatakaForestsub$D,WD=KarnatakaForestsub$WD,errWD = KarnatakaForestsub$sdWD,HDmodel=HDmodel,Dpropag ="chave2004")
    meanAGBperplot<-by(resultMC$AGB_simu,KarnatakaForestsub$plotId,function(x) mean(apply(x, 2, sum)))
    credperplot<-by(resultMC$AGB_simu,KarnatakaForestsub$plotId,function(x) quantile(apply(x,2,sum, na.rm = T), probs = c(0.025, 0.975)))
    credinf<-sapply(credperplot,"[",1)
    credsup<-sapply(credperplot,"[",2)
    ord<-order(meanAGBperplot)
    plot(meanAGBperplot[ord],pch=20,xlab="Plots",ylab="AGB (Mg/ha)",ylim=c(0,max(credsup)),las=1,cex.lab=1.3)
    segments(1:length(ord),credinf[ord],1:length(ord),credsup[ord],col="red")

Propagate error per plot using the local HD model constructed above (modelHD), i.e. independent allometric errors will be assigned to all trees at each iteration, between plots.

    resultMC<-by(KarnatakaForestsub, KarnatakaForestsub$plotId,
                 function(x) AGBmonteCarlo(D=x$D,WD=x$WD,H=x$H,errWD = x$sdWD,
                                           HDmodel=HDmodel,Dpropag ="chave2004"),
                 simplify=F)
    meanAGBperplot<-unlist(sapply(resultMC,"[",1))
    credperplot<-sapply(resultMC,"[",4)
    credinf<-sapply(credperplot,"[",1)
    credsup<-sapply(credperplot,"[",2)
    ord<-order(meanAGBperplot)
    plot(meanAGBperplot[ord],pch=20,xlab="Plots",ylab="AGB (Mg/ha)",ylim=c(0,max(credsup)),las=1,cex.lab=1.3)
    segments(1:length(ord),credinf[ord],1:length(ord),credsup[ord],col="red")

Per plot using the Feldpaush regional HD averaged model (code only given)

    resultMC<-by(KarnatakaForestsub, KarnatakaForestsub$plotId,
                 function(x) AGBmonteCarlo(D=x$D,WD=x$WD,errWD=x$sdWD, H=x$Hfeld,
                                           errH=x$HfeldRSE, Dpropag="chave2004"),
                 simplify=F)
    meanAGBperplot<-unlist(sapply(resultMC,"[",1))
    credperplot<-sapply(resultMC,"[",4)
    credinf<-sapply(credperplot,"[",1)
    credsup<-sapply(credperplot,"[",2)
    ord<-order(meanAGBperplot)
    plot(meanAGBperplot[ord],pch=20,xlab="Plots",ylab="AGB (Mg/ha)",ylim=c(0,max(credsup)),las=1,cex.lab=1.3)
    segments(1:length(ord),credinf[ord],1:length(ord),credsup[ord],col="red")

Per plot using the Chave et al. (2014) Equation 7 (code only given)

    resultMC<-by(KarnatakaForestsub, KarnatakaForestsub$plotId,
                 function(x)AGBmonteCarlo(D=x$D,WD=x$WD,errWD=x$sdWD,
                                          coord=cbind(x$long,x$lat),
                                          Dpropag="chave2004"),
                 simplify=F)
    meanAGBperplot<-unlist(sapply(resultMC,"[",1))
    credperplot<-sapply(resultMC,"[",4)
    credinf<-sapply(credperplot,"[",1)
    credsup<-sapply(credperplot,"[",2)
    ord<-order(meanAGBperplot)
    plot(meanAGBperplot[ord],pch=20,xlab="Plots",ylab="AGB (Mg/ha)",ylim=c(0,max(credsup)),las=1,cex.lab=1.3)
    segments(1:length(ord),credinf[ord],1:length(ord),credsup[ord],col="red")

# 7 Some tricks

## 7.1 Mixing measured and estimated height values

If you want to use a mix of directly-measured height and of estimated ones, you may do the following steps.

1 Build a vector of H and RSE where we assume an error of 0.5 m on directly measured trees

    NouraguesHD$Hmix<-NouraguesHD$H
    NouraguesHD$RSEmix<-0.5
    filt<-is.na(NouraguesHD$Hmix)
    NouraguesHD$Hmix[filt]<- retrieveH(NouraguesHD$D,model = HDmodel)$H[filt]
    NouraguesHD$RSEmix[filt]<-HDmodel$RSE

2 Apply the AGBmonteCarlo by setting the height values and their errors (which depend on wether the tree was directly measured or estimated)

    resultMC<-by(NouraguesHD, NouraguesHD$plotId,
                 function(x)AGBmonteCarlo(D=x$D,WD=x$WD,errWD=x$sdWD,
                                          H=NouraguesHD$Hmix,errH=NouraguesHD$RSEmix,
                                          Dpropag="chave2004"),
                 simplify=F)
    meanAGBperplot<-unlist(sapply(resultMC,"[",1))
    credperplot<-sapply(resultMC,"[",4)
    credinf<-sapply(credperplot,"[",1)
    credsup<-sapply(credperplot,"[",2)
    ord<-order(meanAGBperplot)
    plot(meanAGBperplot[ord],pch=20,xlab="Plots",ylab="AGB (Mg/ha)",ylim=c(0,max(credsup)),las=1,cex.lab=1.3)
    segments(1:length(ord),credinf[ord],1:length(ord),credsup[ord],col="red")

## 7.2 Add your tricks

Please contact Maxime (maxime.rejou@gmail.com) if you would like to add here a code that may be useful for users (code authorship will be respected)