Forest Inventory: Height estimation and plot summarise

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

library(forestmangr)
library(dplyr)
data(exfm9)
data_ex <- exfm9
data_ex
#> # A tibble: 900 x 14
#>   MAP   PROJECT STRATA GENCODE STRATA_AREA PLANTING_DATE SPACING  PLOT
#>   <chr>   <int>  <int> <fct>         <int> <date>        <fct>   <int>
#> 1 FORE~       1      2 FLO014           45 2007-04-01    3x3         1
#> 2 FORE~       1      2 FLO014           45 2007-04-01    3x3         1
#> 3 FORE~       1      2 FLO014           45 2007-04-01    3x3         1
#> 4 FORE~       1      2 FLO014           45 2007-04-01    3x3         1
#> 5 FORE~       1      2 FLO014           45 2007-04-01    3x3         1
#> 6 FORE~       1      2 FLO014           45 2007-04-01    3x3         1
#> # ... with 894 more rows, and 6 more variables: MEASUREMENT_DATE <date>,
#> #   PLOT_AREA <int>, PIT <int>, DBH <dbl>, TH <dbl>, OBS <fct>

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".

dom_height(df=data_ex,th="TH",dbh="DBH",plot="PLOT",obs="OBS",dom="D")
#> # A tibble: 10 x 2
#>    PLOT    DH
#>   <int> <dbl>
#> 1     1  25.0
#> 2     2  23.8
#> 3     3  20.7
#> 4     4  20.0
#> 5     5  19.6
#> 6     7  24.4
#> # ... with 4 more rows

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:

data_ex <- dom_height(data_ex,"TH","DBH","PLOT","OBS","D",merge_data = TRUE)
head(as.data.frame(data_ex))
#>       MAP PROJECT STRATA GENCODE STRATA_AREA PLANTING_DATE SPACING PLOT
#> 1 FOREST1       1      2  FLO014          45    2007-04-01     3x3    1
#> 2 FOREST1       1      2  FLO014          45    2007-04-01     3x3    1
#> 3 FOREST1       1      2  FLO014          45    2007-04-01     3x3    1
#> 4 FOREST1       1      2  FLO014          45    2007-04-01     3x3    1
#> 5 FOREST1       1      2  FLO014          45    2007-04-01     3x3    1
#> 6 FOREST1       1      2  FLO014          45    2007-04-01     3x3    1
#>   MEASUREMENT_DATE PLOT_AREA PIT  DBH   TH OBS    DH
#> 1       2012-09-04       810   1 15.0 23.8   N 25.04
#> 2       2012-09-04       810   2 13.0 23.8   N 25.04
#> 3       2012-09-04       810   3 15.0 24.7   N 25.04
#> 4       2012-09-04       810   4 13.5 23.3   N 25.04
#> 5       2012-09-04       810   5 15.0 24.3   N 25.04
#> 6       2012-09-04       810   6 14.0 22.4   N 25.04

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:

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)
#>       MAP PROJECT STRATA GENCODE STRATA_AREA PLANTING_DATE SPACING PLOT
#> 1 FOREST1       1      2  FLO014          45    2007-04-01     3x3    1
#> 2 FOREST1       1      2  FLO014          45    2007-04-01     3x3    1
#> 3 FOREST1       1      2  FLO014          45    2007-04-01     3x3    1
#> 4 FOREST1       1      2  FLO014          45    2007-04-01     3x3    1
#> 5 FOREST1       1      2  FLO014          45    2007-04-01     3x3    1
#> 6 FOREST1       1      2  FLO014          45    2007-04-01     3x3    1
#>   MEASUREMENT_DATE PLOT_AREA PIT  DBH   TH OBS    DH       CL Henricksen
#> 1       2012-09-04       810   1 15.0 23.8   N 25.04 24.60879   22.98647
#> 2       2012-09-04       810   2 13.0 23.8   N 25.04 23.42624   20.75027
#> 3       2012-09-04       810   3 15.0 24.7   N 25.04 24.60879   22.98647
#> 4       2012-09-04       810   4 13.5 23.3   N 25.04 23.74891   21.31799
#> 5       2012-09-04       810   5 15.0 24.3   N 25.04 24.60879   22.98647
#> 6       2012-09-04       810   6 14.0 22.4   N 25.04 24.05250   21.87975

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:

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:

 data_ex <- data_ex %>% 
  mutate( TH_EST = ifelse(is.na(TH), CL, TH ), CL=NULL,Henricksen=NULL )
head(data_ex)
#>       MAP PROJECT STRATA GENCODE STRATA_AREA PLANTING_DATE SPACING PLOT
#> 1 FOREST1       1      2  FLO014          45    2007-04-01     3x3    1
#> 2 FOREST1       1      2  FLO014          45    2007-04-01     3x3    1
#> 3 FOREST1       1      2  FLO014          45    2007-04-01     3x3    1
#> 4 FOREST1       1      2  FLO014          45    2007-04-01     3x3    1
#> 5 FOREST1       1      2  FLO014          45    2007-04-01     3x3    1
#> 6 FOREST1       1      2  FLO014          45    2007-04-01     3x3    1
#>   MEASUREMENT_DATE PLOT_AREA PIT  DBH   TH OBS    DH TH_EST
#> 1       2012-09-04       810   1 15.0 23.8   N 25.04   23.8
#> 2       2012-09-04       810   2 13.0 23.8   N 25.04   23.8
#> 3       2012-09-04       810   3 15.0 24.7   N 25.04   24.7
#> 4       2012-09-04       810   4 13.5 23.3   N 25.04   23.3
#> 5       2012-09-04       810   5 15.0 24.3   N 25.04   24.3
#> 6       2012-09-04       810   6 14.0 22.4   N 25.04   22.4

To estimate the volume with bark, we’ll take a previously fitted equation, and save it in a data frame:

tabcoef_vwb <- data.frame(b0=-9.595863,b1=1.889372,b2=0.9071631)
tabcoef_vwb
#>          b0       b1        b2
#> 1 -9.595863 1.889372 0.9071631

And do the same for volume without bark:

tabcoef_vwob <- data.frame(b0=-9.808975,b1=1.918007,b2=0.908154)
tabcoef_vwob
#>          b0       b1       b2
#> 1 -9.808975 1.918007 0.908154

Now we’ll estimate volume, basal area and age:

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)
#>       MAP PROJECT STRATA GENCODE STRATA_AREA PLANTING_DATE SPACING PLOT
#> 1 FOREST1       1      2  FLO014          45    2007-04-01     3x3    1
#> 2 FOREST1       1      2  FLO014          45    2007-04-01     3x3    1
#> 3 FOREST1       1      2  FLO014          45    2007-04-01     3x3    1
#> 4 FOREST1       1      2  FLO014          45    2007-04-01     3x3    1
#> 5 FOREST1       1      2  FLO014          45    2007-04-01     3x3    1
#> 6 FOREST1       1      2  FLO014          45    2007-04-01     3x3    1
#>   MEASUREMENT_DATE PLOT_AREA PIT  DBH   TH OBS    DH TH_EST        CSA
#> 1       2012-09-04       810   1 15.0 23.8   N 25.04   23.8 0.01767146
#> 2       2012-09-04       810   2 13.0 23.8   N 25.04   23.8 0.01327323
#> 3       2012-09-04       810   3 15.0 24.7   N 25.04   24.7 0.01767146
#> 4       2012-09-04       810   4 13.5 23.3   N 25.04   23.3 0.01431388
#> 5       2012-09-04       810   5 15.0 24.3   N 25.04   24.3 0.01767146
#> 6       2012-09-04       810   6 14.0 22.4   N 25.04   22.4 0.01539380
#>    AGE       VWB      VWOB
#> 1 66.1 0.2011052 0.1761617
#> 2 66.1 0.1534627 0.1338786
#> 3 66.1 0.2079921 0.1822010
#> 4 66.1 0.1616611 0.1411803
#> 5 66.1 0.2049342 0.1795194
#> 6 66.1 0.1670810 0.1460599

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.

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))   
#>   PLOT AGE STRATA_AREA PLOT_AREA     DBH       q  TH_EST    DH Indv
#> 1    1  66          45       810 14.1722 14.2196 24.0004 25.04   90
#> 2    2  66          45       810 14.4494 14.5235 23.3723 23.76   90
#> 3    3  63          45       810 12.8539 12.8902 19.9767 20.70   90
#> 4    4  63          51       810 12.2833 12.3285 18.9652 19.98   90
#> 5    5  63          51       810 12.6389 12.6939 18.8658 19.60   90
#> 6    7  66          45       810 14.7584 14.8135 23.9202 24.40   90
#>     Indvha      G    G_ha     VWB   VWB_ha    VWOB  VWOB_ha
#> 1 1111.111 1.4292 17.6450 16.5503 204.3245 14.4789 178.7514
#> 2 1111.111 1.4744 18.2028 16.6473 205.5216 14.5737 179.9222
#> 3 1111.111 1.1614 14.3388 11.5018 141.9973 10.0318 123.8493
#> 4 1111.111 1.0744 13.2638 10.2121 126.0759  8.8958 109.8248
#> 5 1111.111 1.1390 14.0618 10.7609 132.8508  9.3818 115.8243
#> 6 1111.111 1.5339 18.9371 17.6205 217.5366 15.4333 190.5350