---
title: "Forest Inventory: Height estimation and plot summarise"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Forest Inventory: Height estimation and plot summarise}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---
  
```{r, echo = FALSE, message = FALSE}
knitr::opts_chunk$set(collapse = T, comment = "#>")
knitr::opts_chunk$set(fig.width=7, fig.height=5)
options(tibble.print_min = 6L, tibble.print_max = 6L)
library(forestmangr)
library(dplyr)
```

For this example we'll use forest inventory data, calculate volume and other variables for each plot.

```{r}
library(forestmangr)
library(dplyr)
data(exfm9)
data_ex <- exfm9
data_ex
```

The first step is to estimate the height of non measured trees. We'll evaluate two hypsometric models. Henricksen's:
$$ Ln(H) = \beta_0 + \beta_1*Ln(H) $$ 

And Campos & Leite's model, which uses dominant height:
$$ Ln(H) = \beta_0 + \beta_1*\frac{1}{dbh} + \beta_2*Ln(DH) $$ 
In order to use this model, first we'll need to calculate the dominant height for each plot. To do this we'll use the `dom_height` function. In it we'll input the dataset, and height, dbh, plot and observation variables. The observation variable refers to the quality or classification of the tree, i.e., if it's normal, forked, dead or dominant. In addition to that, we'll also input the code used to distinguish dominant trees. In this dataset, the code is `"dom"`.
```{r}
dom_height(df=data_ex,th="TH",dbh="DBH",plot="PLOT",obs="OBS",dom="D")
```
Now that we've seen the dominant height of each plot, we can run the function again, but this time set the `merge_data` argument as `TRUE`, to bind the variable to the data directly:
```{r}
data_ex <- dom_height(data_ex,"TH","DBH","PLOT","OBS","D",merge_data = TRUE)
head(as.data.frame(data_ex))
```

Now we can fit the hypsometric models. We'll fit them using the `lm_table` function. the function `forestmangr::inv` will allow us to fit the Campos & Leite model without creating a new variable for 1/dbh:
```{r}
data_ex <- data_ex %>% 
  lm_table(log(TH) ~ inv(DBH) + log(DH),output="merge_est",est.name="CL") %>% 
  lm_table(log(TH) ~ log(DBH), output="merge_est",est.name="Henricksen") 
head(data_ex)
```
Ps: The `lm_table` function checks if the model has log in the y variable, and if it does, it removes it automatically when estimating variables. Because of that, there's no need to calculate the exponential for the estimated variables.

We'll evaluate the quality of the fitted values using `resid_plot`. Non measured trees will be removed automatically:
```{r, warning=FALSE}
resid_plot(data_ex, "TH", "CL","Henricksen")
```
Campos & Leite's model was superior, thus, it will be used.

Now we can estimate the height of non measured trees, using `dplyr::mutate` and `ifelse`:
```{r}
 data_ex <- data_ex %>% 
  mutate( TH_EST = ifelse(is.na(TH), CL, TH ), CL=NULL,Henricksen=NULL )
head(data_ex)
```

To estimate the volume with bark, we'll take a previously fitted equation, and save it in a data frame:
```{r}
tabcoef_vwb <- data.frame(b0=-9.595863,b1=1.889372,b2=0.9071631)
tabcoef_vwb
```
And do the same for volume without bark:
```{r}
tabcoef_vwob <- data.frame(b0=-9.808975,b1=1.918007,b2=0.908154)
tabcoef_vwob
```

Now we'll estimate volume, basal area and age:
```{r}
data_ex <- data_ex %>% 
    mutate(CSA = pi*DBH^2/40000,
         AGE = as.numeric(MEASUREMENT_DATE - PLANTING_DATE) / 30,
         VWB = exp(tabcoef_vwb$b0 + tabcoef_vwb$b1*log(DBH) + tabcoef_vwb$b2*log(TH_EST) ),
         VWOB = exp(tabcoef_vwob$b0 + tabcoef_vwob$b1*log(DBH) + tabcoef_vwob$b2*log(TH_EST) ) )
head(data_ex)
```


To summarise the plots, we'll use `plot_summarise`. We'll input the dataset, and variables such as plot, plot area, dbh, height, total area, volume with bark, volume without bark, dominant height and age.
With this, the function will calculate for each plot the mean diameter, quadratic diameter, mean height and mean dominant height. It will also calculate to total number of individuals, basal area and volume with and without bark for each plot, and extrapolate it to hectares.
```{r}
tab_invt <- plot_summarise(df=data_ex,plot="PLOT",plot_area="PLOT_AREA",dbh="DBH",                           th="TH_EST",total_area="STRATA_AREA",vwb="VWB",vwob="VWOB",dh="DH",age="AGE")
head(as.data.frame(tab_invt))   
```