Using BiomasaFP and BIOMASS together

Martin Sullivan

2021-06-16

The BIOMASS R package provides an alternative tool for calculating aboveground biomass. Unlike BiomasaFP, it hasn't been designed to run of ForestPlots.net outputs, which means it requires the user to do a bit more coding, but will be easier to use if your data are not in ForestPlots.net as it doesn't assume that there are particular columns in the data. While for biomass calculations you could use either package, not all functionality overlaps. Only the BIOMASS package provides error propagated estimates of aboveground biomass, while only BiomasaFP provides measures of forest dynamics such as woody productivity.

While you will probably use one package or the other for your analysis, there are a few useful functions that BIOMASS provides that can be integrated, with a bit of work, into an analysis with BiomasaFP.

Wood density lookup

The ForestPlots.net database has an inbuilt copy of the Chave/ Zanne et al. 2009 wood density database, which you can access using the query library. This wood density data is one of the required inputs for the mergefp function:

dat<-mergefp(trees,md,wd)

But the BIOMASS R package also provides a way of querying this database through the getWoodDensity function. You might want to do this if you are having problems running the ForestPlots query for wood density, or if you want to change some of the parameters for matching wood density (e.g. not matching at family level, not matching at continent level).

To use the getWoodDensity function, our first challange is to split the species column into seperate genus and species columns:

dat$Genus<-sapply(dat$Species,function(x)strsplit(x," ")[[1]][1])
dat$Sp<-sapply(dat$Species,function(x)strsplit(x," ")[[1]][2])
head(dat$Species)
#> [1] "Protium hebetatum" "Protium hebetatum" "Protium hebetatum"
#> [4] "Protium hebetatum" "Caryocar glabrum"  "Caryocar glabrum"
head(dat$Genus)
#> [1] "Protium"  "Protium"  "Protium"  "Protium"  "Caryocar" "Caryocar"
head(dat$Sp)
#> [1] "hebetatum" "hebetatum" "hebetatum" "hebetatum" "glabrum"   "glabrum"

Before we run the getWoodDensity function, an important step is to reduce the data to just unique tree IDs (if you don't do this, bad things will happen when you try to run mergefp).

dat2<-unique(dat[,c("TreeID","PlotID","PlotCode","PlotViewID","Genus","Sp","Family")])
nrow(dat)
#> [1] 8099
nrow(dat2)
#> [1] 1309

Finally, we are ready to run getWoodDensity.

wd.biomass<-getWoodDensity(dat2$Genus,dat2$Sp,stand=dat2$PlotViewID,family=dat2$Family)
#> The reference dataset contains 16467 wood density values
#> Your taxonomic table contains 361 taxa
str(wd.biomass)
#> 'data.frame':    1309 obs. of  7 variables:
#>  $ family : chr  "Burseraceae" "Caryocaraceae" "Lecythidaceae" "Lecythidaceae" ...
#>  $ genus  : chr  "Protium" "Caryocar" "Eschweilera" "Eschweilera" ...
#>  $ species: chr  "hebetatum" "glabrum" "tessmannii" "coriacea" ...
#>  $ meanWD : num  0.568 0.654 0.789 0.852 0.533 ...
#>  $ sdWD   : num  0.0941 0.0708 0.0708 0.0708 0.0708 ...
#>  $ levelWD: chr  "genus" "species" "species" "species" ...
#>  $ nInd   : int  31 7 1 12 3 4 5 13 5 32 ...

Now append the mean wood density for each stem to our data, and then restrict the columns so they match the ForestPlots.net wood density output.

dat2$WD<-wd.biomass$meanWD
wd2<-data.frame("PlotID"=dat2$PlotID,"PlotCode"=dat2$PlotCode,"PlotViewID"=dat2$PlotViewID,"TreeID"=dat2$TreeID,"WD"=dat2$WD)

We can now proceed with mergefp as before.

dat3<-mergefp(trees,md,wd2)

Height-diameter allometry

BIOMASS has nice functions for fitting height-diameter models and plotting fits of different models against the input data.

The first thing we have to do is to covert the diameter data from mm into cm (BiomasaFP height-diameter functions do this in the function, but BIOMASS expects the inputs in cm).

dat$D.cm<-dat$D4/10

You can fit height-diameter models across multiple plots at once in BIOMASS, but you only get the nice graphical output if you look at one plot at a time, so we'll look at each plot in turn.

with(dat[dat$PlotCode=="JEN-11",],modelHD(D.cm,Height))
#> To build a HD model you must use the parameter 'method' in this function
#>      method  color      RSE    RSElog  Average_bias
#> 1      log1   blue 4.169132 0.1634860  0.0015729103
#> 2      log2  green 4.067784 0.1619454  0.0008976036
#> 3   weibull orange 4.054683        NA -0.0001152028
#> 4 michaelis purple 4.043162        NA -0.0001060615
with(dat[dat$PlotCode=="CAP-09",],modelHD(D.cm,Height))
#> To build a HD model you must use the parameter 'method' in this function
#>      method  color      RSE    RSElog  Average_bias
#> 1      log1   blue 6.039085 0.2232796  0.0039239120
#> 2      log2  green 5.999312 0.2218205  0.0031984494
#> 3   weibull orange 5.981388        NA -0.0002227037
#> 4 michaelis purple 5.960755        NA -0.0007356204

We can use the RSE and bias values to select the best model, and also use the graphical outputs.

Let's go with the Michaelis-Menten model to use for both plots. Refit the models specifying which type of model to use, and extract the coefficients.

m1<-with(dat[dat$PlotCode=="JEN-11",],modelHD(D.cm,Height,method="michaelis"))$coefficients
m2<-with(dat[dat$PlotCode=="CAP-09",],modelHD(D.cm,Height,method="michaelis"))$coefficients

Now our challenge is to get this output into a format that can be used with BiomasaFP. The BiomasaFP R package takes a height parameters dataset with the tree ID, a, b and c parameters for each tree, and the type of model for each tree. The code below sets that up:

h.params2<-data.frame("TreeID"=unique(dat$TreeID),"a_par"=NA,"b_par"=NA,"Model"=NA)
h.params2$a_par[h.params2$TreeID%in%dat$TreeID[dat$PlotCode=="JEN-11"]]<-m1[1,1]
h.params2$b_par[h.params2$TreeID%in%dat$TreeID[dat$PlotCode=="JEN-11"]]<-m1[2,1]
h.params2$Model[h.params2$TreeID%in%dat$TreeID[dat$PlotCode=="JEN-11"]]<-"MicMent"
h.params2$a_par[h.params2$TreeID%in%dat$TreeID[dat$PlotCode=="CAP-09"]]<-m2[1,1]
h.params2$b_par[h.params2$TreeID%in%dat$TreeID[dat$PlotCode=="CAP-09"]]<-m2[2,1]
h.params2$Model[h.params2$TreeID%in%dat$TreeID[dat$PlotCode=="CAP-09"]]<-"MicMent"
# Need a c_par variable
h.params2$c_par<-NA

head(h.params2)
#>    TreeID    a_par    b_par   Model c_par
#> 1 1186624 46.80975 22.88332 MicMent    NA
#> 2 1186625 46.80975 22.88332 MicMent    NA
#> 3 1186626 46.80975 22.88332 MicMent    NA
#> 4 1186627 46.80975 22.88332 MicMent    NA
#> 5 1187406 46.80975 22.88332 MicMent    NA
#> 6 1187407 46.80975 22.88332 MicMent    NA

Now we can run the SummaryAGWP function with the new height data.

dynamics.mic<-SummaryAGWP(dat,AGBChv14,height.data=h.params2)

summary(dynamics.mic)
#>    PlotViewID      Census.No         AGB.ha         Stems.ha    
#>  Min.   :112.0   Min.   :1.000   Min.   :300.5   Min.   :486.0  
#>  1st Qu.:112.0   1st Qu.:2.250   1st Qu.:309.1   1st Qu.:494.2  
#>  Median :112.0   Median :4.000   Median :313.8   Median :597.0  
#>  Mean   :124.5   Mean   :4.286   Mean   :343.8   Mean   :563.2  
#>  3rd Qu.:147.0   3rd Qu.:5.750   3rd Qu.:399.2   3rd Qu.:601.8  
#>  Max.   :147.0   Max.   :9.000   Max.   :426.6   Max.   :620.0  
#>                                                                 
#>    AGWPrec.ha       Recruit.ha      AGBmort.ha      Mortality.ha  
#>  Min.   :0.1553   Min.   : 3.00   Min.   : 3.364   Min.   : 9.00  
#>  1st Qu.:0.4953   1st Qu.: 8.50   1st Qu.: 7.305   1st Qu.:11.00  
#>  Median :1.1568   Median :20.00   Median :13.409   Median :16.50  
#>  Mean   :0.9788   Mean   :16.58   Mean   :15.678   Mean   :17.83  
#>  3rd Qu.:1.4059   3rd Qu.:24.00   3rd Qu.:16.063   3rd Qu.:22.75  
#>  Max.   :1.7759   Max.   :28.00   Max.   :49.721   Max.   :32.00  
#>  NA's   :2        NA's   :2       NA's   :2        NA's   :2      
#>   AGWPsurv.ha     SurvivingStems.ha UnobsAGWPmort.ha     Mean.WD      
#>  Min.   : 5.087   Min.   :  0.0     Min.   :0.03923   Min.   :0.6619  
#>  1st Qu.: 7.444   1st Qu.:469.8     1st Qu.:0.07197   1st Qu.:0.6659  
#>  Median :10.708   Median :573.0     Median :0.11713   Median :0.6683  
#>  Mean   :15.401   Mean   :469.7     Mean   :0.26307   Mean   :0.6679  
#>  3rd Qu.:16.631   3rd Qu.:586.0     3rd Qu.:0.28086   3rd Qu.:0.6710  
#>  Max.   :47.034   Max.   :609.0     Max.   :0.93066   Max.   :0.6714  
#>  NA's   :2                          NA's   :2                         
#>  Recruitment.stem.year Mortality.stem.year UnobsAGWPrec.ha    UnobsRecruits.ha 
#>  Min.   :0.004414      Min.   :0.005830    Min.   :0.004012   Min.   :0.07946  
#>  1st Qu.:0.006499      1st Qu.:0.008194    1st Qu.:0.008933   1st Qu.:0.16853  
#>  Median :0.010767      Median :0.012093    Median :0.010545   Median :0.19751  
#>  Mean   :0.012352      Mean   :0.014804    Mean   :0.012099   Mean   :0.22257  
#>  3rd Qu.:0.019355      3rd Qu.:0.018850    3rd Qu.:0.014637   3rd Qu.:0.26515  
#>  Max.   :0.023914      Max.   :0.036578    Max.   :0.022154   Max.   :0.39922  
#>  NA's   :2             NA's   :2           NA's   :2          NA's   :2        
#>  CensusInterval     AGWP.ha        AGWP.ha.year     PlotCode        
#>  Min.   :0.929   Min.   : 0.000   Min.   :3.737   Length:14         
#>  1st Qu.:1.174   1st Qu.: 6.272   1st Qu.:5.492   Class :character  
#>  Median :2.033   Median :10.982   Median :6.264   Mode  :character  
#>  Mean   :3.007   Mean   :14.275   Mean   :5.966                     
#>  3rd Qu.:3.903   3rd Qu.:15.470   3rd Qu.:6.672                     
#>  Max.   :7.978   Max.   :49.751   Max.   :7.326                     
#>  NA's   :2                        NA's   :2