README.Rmd

---
output: github_document 
editor_options: 
  chunk_output_type: inline
---
<!-- README.md is generated from README.Rmd. Please edit that file -->

```{r, include = FALSE}
knitr::opts_chunk$set(
  cache = TRUE,
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-"
)

set.seed(1)
```

# <img src="https://i.imgur.com/vTLlhbp.png" align="right" height=88 /> Calculate biomass

[![lifecycle](https://img.shields.io/badge/lifecycle-experimental-orange.svg)](https://www.tidyverse.org/lifecycle/#experimental)
[![Travis build status](https://travis-ci.org/forestgeo/fgeo.biomass.svg?branch=master)](https://travis-ci.org/forestgeo/fgeo.biomass)
[![Coverage status](https://coveralls.io/repos/github/forestgeo/fgeo.biomass/badge.svg)](https://coveralls.io/r/forestgeo/fgeo.biomass?branch=master)
[![CRAN status](https://www.r-pkg.org/badges/version/fgeo.biomass)](https://cran.r-project.org/package=fgeo.biomass)

The goal of fgeo.biomass is to calculate biomass using [ForestGEO](https://forestgeo.si.edu/) data and equations from either the [BIOMASS package](https://CRAN.R-project.org/package=BIOMASS) or the [allodb package](https://forestgeo.github.io/allodb/). 

* The BIOMASS package is applicable to tropical forests. It was first [published on CRAN in 2016](https://cran.r-project.org/) and on [Methods on Ecology and Evolution in 2017](https://besjournals.onlinelibrary.wiley.com/doi/abs/10.1111/2041-210X.12753). fgeo.biomass provides the main features of BIOMASS with a simpler interface, consistent with all [fgeo packages](https://forestgeo.github.io/fgeo/).

* The allodb package is work in progress, and aims to provide expert-selected allometric equations, both for tropical and temperate forests. fgeo.biomass provides a simple interface to automate the process of finding the right equation(s) for each stem and computing biomass. 

## Installation

Install the development version of **fgeo.biomass** with:

```
# install.packages("devtools")
devtools::install_github("forestgeo/fgeo.biomass")
```

## Setup

In addition to the fgeo.biomass package we will use dplyr and ggplot2 for data wrangling and plotting.

```{r, message=FALSE}
library(ggplot2)
library(dplyr)
library(fgeo.biomass)
```

## fgeo.biomass wrapping BIOMASS

We'll use data from the [Barro Colorado Island, Panama](https://forestgeo.si.edu/sites/neotropics/barro-colorado-island) (BCI). We first pick alive trees and drop missing `dbh` values as we can't calculate biomass for them.

```{r}
if (!requireNamespace("bciex", quietly = TRUE)) {
  stop(
    "For this example, you must first install the bciex package with:\n",
    "devtools::install_github('forestgeo/bciex')"
  )
}

bci_tree <- as_tibble(bciex::bci12t7mini) %>% 
  filter(status == "A", !is.na(dbh))
bci_tree
```

We also need species data.

```{r}
bci_species <- as_tibble(bciex::bci_species)
bci_species
```

`add_tropical_biomass()` adds biomass to your census data.

```{r}
biomass <- add_tropical_biomass(bci_tree, bci_species)
biomass
```

You may also provide a specific `region` or `latitude` and `longitude`.

```{r}
biomass <- add_tropical_biomass(
  bci_tree, 
  bci_species,
  latitude = 9.154965, 
  longitude = -79.845884
)

biomass %>% 
  select(biomass, everything())
```

`propagate_errors()` allows you to propagate errors.

```{r}
str(
  propagate_errors(biomass)
)
```

`model_height()` allows you to create a height model, which you can use to propagate height errors. This is what the entire pipeline looks like:

```{r}
model <- model_height(bci_tree)

errors <- bci_tree %>% 
  add_tropical_biomass(bci_species) %>% 
  propagate_errors(height_model = model)

str(errors)
```

If you pass `latitude` and `longitude` to `add_tropical_biomass(), and then you pass a `height_model` to `propagate_errors()`, then you will need to ignore the coordinates. On an interactive session, you should see something like this:

![](https://i.imgur.com/dhHCYJN.png)

```{r}
if (interactive()) {
  errors <- bci_tree %>% 
    add_tropical_biomass(
      bci_species, 
      latitude = 9.154965, 
      longitude = -79.845884
    ) %>% 
    propagate_errors(height_model = model)
  
  str(errors)
}
```

`add_wood_density()` adds wood density to your census data. It is not limited to tropical forests, and has support for all of these regions: `r glue::glue_collapse(fgeo.biomass:::wd_regions(), sep = ", ", last = ", and ")`.

```{r}
wood_density <- add_wood_density(bci_tree, bci_species)

wood_density %>% 
  select(starts_with("wd_"), everything())
```

The BIOMASS package provides a tool to correct taxonomic names. fgeo.biomass does not include that feature. You may use BIOMASS directly or the more focused [taxize package](https://cran.r-project.org/web/packages/taxize/taxize.pdf).

##  fgeo.biomass wrapping allodb

## Warning

These features are not ready for research. We are now building a [Minimum Viable Product](https://en.wikipedia.org/wiki/Minimum_viable_product), with just enough features to collect feedback from alpha users and redirect our effort. The resulting biomass is still meaningless.

We'll use the `add_biomass()` with these inputs:

1. A ForestGEO-like _stem_ or _tree_ table.
2. A _species_ table (internally used to look up the Latin species names from the species codes in the `sp` column of the census table).

We'll use data from the [Smithsonian Conservation Biology Institute, USA](https://forestgeo.si.edu/sites/north-america/smithsonian-conservation-biology-institute) (SCBI). We first pick alive trees and drop missing `dbh` values as we can't calculate biomass for them.

```{r}
census <- fgeo.biomass::scbi_tree1 %>% 
  filter(status == "A", !is.na(dbh))

census
```

We now use `add_biomass()` to add biomass to our census dataset. 

```{r}
species <- fgeo.biomass::scbi_species

with_biomass <- census %>% 
  add_biomass(species, site = "SCBI")
```

We are warned that we are using a tree-table (as opposed to a stem-table), and informed about how to interpret the resulting `biomass` values for trees and shrubs.

Some equations couldn't be found. There may be two reasons:

* Some stems in the data belong to species with no matching species in allodb.
* Some stems in the data belong to species that do match species in allodb but the available equations were designed for a dbh range that doesn't include actual dbh values in the data.

Here are the most interesting columns of the result:

```{r}
with_biomass %>% 
  select(treeID, species, biomass)
```

Let's now visualize the relationship between `dbh` and b`biomass` by `species` (black points), in comparison with `agb` (above ground biomass) values calculated with allometric equations for tropical trees (grey points). 

```{r, fig.height=14}
with_biomass %>% 
  # Convert agb from [Mg] to [kg]
  mutate(agb_kg = agb * 1e3) %>% 
  ggplot(aes(x = dbh)) +
  geom_point(aes(y = agb_kg), size = 1.5, color = "grey") +
  geom_point(aes(y = biomass), size = 1, color = "black") +
  facet_wrap("species", ncol = 4) +
  ylab("Reference `agb` (grey) and calculated `biomass` (black) in [kg]") +
  xlab("dbh [mm]") +
  theme_bw()
```

Above, the species for which `biomass` couldn't be calculated show no black points, although they do show grey reference-points.

To better understand the distribution of `biomass` values for each species we can use a box-plot.

```{r}
with_biomass %>% 
  ggplot(aes(species, biomass)) +
  geom_boxplot() +
  ylab("biomass [kg]") +
  coord_flip()
```

For some species the maximum `dbh` for which `biomass` was calculated is much lower than the maximum `dbh` value for which the reference `agb` was calculated. This is because most equations in __allodb__ are defined for a specific range of `dbh` values. Eventually __allodb__ might provide equations beyond the `dbh` limits currently available.

To explore this issue, here we use `add_component_biomass()` which allows us to see intermediary results that `add_biomass()` doesn't show.

```{r}
detailed_biomass <- suppressWarnings(suppressMessages(
  add_component_biomass(census, species, site = "SCBI")
))

# Maximum `dbh` values by species
max_by_species <- detailed_biomass %>% 
  select(species, dbh_max_mm) %>% 
  group_by(species) %>% 
  arrange(desc(dbh_max_mm)) %>% 
  filter(row_number() == 1L) %>% 
  ungroup()

# `dbh` is above the maximum limit, so `biomass` is missing (agb has a value)
detailed_biomass %>% 
  filter(dbh > 1000) %>% 
  select(-dbh_max_mm) %>% 
  left_join(max_by_species) %>% 
  mutate(agb_kg = agb * 1e3) %>%
  select(species, biomass, agb, dbh, dbh_max_mm) %>% 
  arrange(species) %>%
  print(n = Inf)
```

## Biomass via BIOMASS versus allodb

```{r}
temperate_biomass <- add_biomass(census, species, site = "scbi")

# Warning: Aplying tropical equations to a temperate forest for comparison
tropical_biomass <- add_tropical_biomass(census, species)

dbh_biomsss <- tibble(
  dbh = temperate_biomass$dbh,
  species = temperate_biomass$species,
  temperate_biomass = temperate_biomass$biomass, 
  tropical_biomass = tropical_biomass$biomass
)
```

```{r, fig.height=14}
dbh_biomsss %>% 
  ggplot(aes(x = dbh)) +
  geom_point(aes(y = tropical_biomass), size = 1.5, color = "grey") +
  geom_point(aes(y = temperate_biomass), size = 1) +
  facet_wrap("species", ncol = 4) +
  ylab("Biomass [kg] (via the BIOMASS (grey) and allodb (black) packages)") +
  xlab("dbh [mm]") +
  theme_bw()
```

## General information

* [Getting help](SUPPORT.md).
* [Contributing](CONTRIBUTING.md).
* [Contributor Code of Conduct](CODE_OF_CONDUCT.md).