--- title: "Yield Tables in ForestElementsR" author: "Peter Biber" output: rmarkdown::html_document: toc: true toc_depth: 3 toc_float: true bibliography: REFERENCES.bib vignette: > %\VignetteIndexEntry{Yield Tables in ForestElementsR} %\VignetteEncoding{UTF-8} %\VignetteEngine{knitr::rmarkdown} editor_options: markdown: wrap: 72 --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` ## 1 Introduction Yield tables might be the most important tools forest science has ever provided to forest practice. They represent an invaluable wealth of empirical data that has been condensed into tabular descriptions of forest stand growth that are easy to work with. Since about three decades, however, they have been increasingly frowned upon for several reasons, among them their limited applicability to mixed species stands, their inherent restriction to even aged stands only, and their representation of forest growth under environmental conditions that have meanwhile changed, and cannot be considered as stable as they were at the time when the data behind the tables were collected. Nevertheless, when handled wisely, yield tables can still be highly valuable assets in our tool kit. Therefore, we developed a generic implementation of yield tables in *ForestElementsR* that is designed to work equally well with yield tables of different origin with different specifications and descriptive variables. A collection of yield tables for important tree species in Central Europe has already been implemented and comes readily with *ForestElements*; this collection will certainly keep growing. All readily available tables are listed in the documentation, you can look them up with ```?yield_tables```. The purpose of this vignette is to explain how to work with our yield table implementation. ```{r} # Before we start, let's attach the package library(ForestElementsR) ``` ## 2 Yield Table Representation in ForestElementsR Let's have a look to a classic yield table, Wiedemann's 1943 table for Scots pine (*Pinus sylvestris*). The bulk information used for our implementation was the version published by the Bavarian State Forest Administration [@hilfstafeln_bavaria_2018], some information that was missing there, was taken from the version in Schober's yield table collection [@ertragstafeln_schober_1975] in addition. As any yield table in ForestElementsR, the Wiedemann table for Scots pine is an S3 object of class ```fe_yield_table```. its representation in the package is called ```fe_ytable_pine_wiedemann_moderate_1943```. ```{r} class(fe_ytable_pine_wiedemann_moderate_1943) ``` Such an object is essentially a list containing some metadata followed by the actual data. A basic overview can be obtained with: ```{r} fe_ytable_pine_wiedemann_moderate_1943 |> summary() ``` The first six list elements represent the metadata, the seventh element, ```values```, comprises the actual data. ### 2.1 Metadata Let's have a look at the metadata first: ```{r} fe_ytable_pine_wiedemann_moderate_1943[1:6] # Elements 1-6 contain metadata ``` The list element ```name_orig``` contains the established original name of the yield table in the language it was originally published in. The second item, ```name_international``` is an English translation of the original name. The item ```site_indexes``` is a numeric vector of all site index values that are actually represented in the table with values. The Wiedemann pine table is using the traditional European approach of *relative site indexing*, i.e. site index 1.0 stands for the best site conditions, while 4.0 or greater represent considerably weak sites. A site index of 1.5 would be in the middle between the best and the second best site represented in the table. In its original form, this traditional site index numbering actually uses roman numerals before and arabic ones after the decimal point (e.g. II.3 = 2.3), but implementing this simply was not worth the effort. More important, our implementation works equally with yield tables that follow the approach of absolute site indexing (which is preferable from our point of view, but not frequently found in German yield tables). Have a look at the ```site_indexes``` entry of the yield table by Ernst Assmann and Friedrich Franz for Norway spruce (as taken from @hilfstafeln_bavaria_2018): ```{r} fe_ytable_spruce_assmann_franz_mean_yield_level_1963$site_indexes ``` These site indexes represent a stand's dominant height in meters at an age of 100 years. So, a stand of site index 40 after Assmann and Franz would be expected to be 40 m high when it is 100 years old. As you will see below, we have implemented the option to express site quality in terms of such an absolute site index also for yield tables with a relative site index system. Back to the Wiedemann pine table, the next element of the metadata is called ```site_index_variable```. This is a character vector with - typically - one element that provides the name of the variable which is used to determine the site index. As for most German yield tables, for the Wiedemann table it is the quadratic mean height in meters, abbreviated as "h_q_m". This means when the quadratic mean height (and the age) of a Scots pine stand is known, its site index can be determined with this yield table. Again, the Assmann-Franz table can serve as an example for other possibilities: ```{r} fe_ytable_spruce_assmann_franz_mean_yield_level_1963$site_index_variable ``` This table was designed for site indexing based on the stands dominant height, here defined as the quadratic mean height of the 100 tallest trees per hectare, called "h_100_m" in our implementation. But in addition, this table also contains the usual quadratic mean stand height, "h_q_m", so we allow both variables to be used for finding a stand's site index. The next metadata element we called ```mai_variable```. This relates to the concept of expressing site indexes in terms of the mean annual increment (the German term for this is "DGZ-Bonitierung"). The idea of that approach is to use a yield table with the classic input (stand height, age), but to express the site quality in terms of the expected mean annual increment ("m.a.i.", German abbreviation "dGz") at a given age or at its maximum (see more details below). The last metadata item, ```age_coverage``` is a numeric vector of all ages for which the table actually provides data. Note, that this is an "outer" age coverage, i.e. not all site indexes given in the table must necessarily cover the whole range of values. Often, poor site indexes start in the table at a higher age than better ones. Sometimes, the age interval might be different for different site indexes. E.g. the Wiedemann pine table usually provides values in five year age intervals, but in ten-year intervals only for site indexes 5.5 and 6.0. This is also not a problem with the age coverage to be given in the metadata; it always relates to the highest temporal resolution available in the table. ### 2.2 Actual Yield Table Data The yield table data themselves are stored in the list element called ```values```. This itself is a named list of matrices, one matrix for each variable given in the yield table, the list element names being the names of the variables. Note, that ```fe_yield_table``` objects are not required to have all an identical set of variables. In order to find out which variables are given in a yield table, you could obtain their names as follows: ```{r} fe_ytable_pine_wiedemann_moderate_1943$values |> names() ``` The order of the variables is not related to their importance, e.g. the first variable listed here, *h_q_m_si_pluis_025*, is the lower threshold height for a given site index as by default provided in the Bavarian state forest yield table editions (@hilfstafeln_bavaria_2018, @hilfstafeln_bavaria_1990). For electronic site indexing this variable is not required, we kept it, however, test and documentation purposes. Of central importance, in contrast, are the variables which are used for the actual site indexing. In the Wiedemann pine table, this is ```h_q_m"```, the quadratic mean stand height, as indicated in the metadata. As for any other variable in the yield table, it can be simply displayed as follows: ```{r} fe_ytable_pine_wiedemann_moderate_1943$values$h_q_m ``` The columns of the matrix correspond to the site indexes given in the metadata. They are represented in the column names with the prefix "si_". Each matrix row stands for a specific stand age, whereby the stand ages are the row names. The NA values are places where no value is given in the table, typically with low site indexes at younger ages or with broader age intervals. All other variables in the table are given as a matrices of exactly the same structure, take e.g. the periodic annual increment ```pai_m3_ha_yr```: ```{r} fe_ytable_pine_wiedemann_moderate_1943$values$pai_m3_ha_yr ``` Before we come to visualizing yield table contents, let us shortly explain the variable names used in the Wiedemann pine table, as we are sticking as close as possible to this naming concept in any yield table we make a part of ForestElementsR. Note, however, that not all of these variables are necessarily contained in all tables we have included so far; some tables might also contain additional variables and (unfortunately) use other than metric units. For creating a working ```fe_yield_table``` object it is not necessary to stick to a certain naming concept although we highly recommend to do so. Let's explain the variable names given above, beginning whith those that are common in most yield tables: * **h_q_m**: Quadratic mean stand height (m) * **d_q_cm**: Quadratic mean diameter at breast height (cm) * **n_ha**: Number of trees per ha * **ba_m2_ha**: Stand basal area (m²/ha) * **v_m3_ha** : Stand wood volume (m³/ha) * **tvp_m³_ha**: Total volume production (m³/ha) * **mai_m³_ha_yr**: Mean annual volume increment (m³/ha/yr), i.e. the total volume production divided by the stand age * **pai_m3_ha_yr**: Periodic mean annual increment (m³/ha/yr), i.e. the difference of the total volume production for two subsequent points in time divided by the corresponding time interval * **pai_perc_yr**: Annual relative periodic increment (%/yr) As the Wiedemann pine table was imported from the edition provided by the Bavarian State Forest Service [@hilfstafeln_bavaria_2018], there are some peculiarities: Some of the variable names carry the prefix "red_". This means, while volume-related values are usually, to be understood as m³ standing wood over bark, the prefix "red_" ("reduced") indicates, that such wood volumes are to be understood under bark and after subtracting losses due to the harvest. This explains the variable names **red_v_m3_ha**, **red_mai_m3_ha_yr**, and **red_pai_m3_ha_yr**. Another peculiarity of the Bavarian table edition is the scarce information about the removal stand. For a given stand age, it provides only the harvest volume for the subsequent decade, also under bark, and after harvest, hence the name **red_pre_yield_m3_ha_10yr**. More information about the removal stand in yield tables will be provided (if available from other table editions such as @ertragstafeln_schober_1975) in future versions of *ForestElementsR*. ## 3 Visualizing Yield Table Contents For visualization purposes whe have written a plot method for ```fe_yield_table``` objects. ```{r, fig.width=4.9, fig.height=4.9, fig.align="center"} # Plot the quadratic mean height (site index fan) fe_ytable_pine_wiedemann_moderate_1943 |> plot() ``` When the plot function is called without further arguments, it displays the first variable listed as ```site_index_variable``` in the object's metadata. Typically this results in a display of a "site index fan". All other variables contained in the ```fe_yield_table``` object can be plotted by specifiying their name in the call: ```{r, fig.width=4.9, fig.height=4.9, fig.align="center"} # Plot the periodic annual increment fe_ytable_pine_wiedemann_moderate_1943 |> plot(variable = "pai_m3_ha_yr") ``` The plots are *ggplot* objects whose appearance can be partly modified in a post-hoc manner if required, e.g. ```{r, fig.width=4.9, fig.height=4.9, fig.align="center"} # Plot the tree number per ha and modify the plot's appearance # 1 catch the plot v_plot <- fe_ytable_pine_wiedemann_moderate_1943 |> plot(variable = "n_ha") # 2 make adjustments as allowed by ggplot2, e.g. v_plot + ggplot2::theme_classic() + # Appearance similar to R base graphics ggplot2::scale_y_log10() + # Logarithmic scale for vertical axis ggplot2::ylab("Tree Number per hectare") # More explaining axis label ``` ## 4 Practical Work With Yield Tables in ForestElementsR ### 4.1 Site Indexing The fundamental task in working with yield tables is finding the site index of a given stand. This requires two input variables, the stand's age, and the stand's height, defined as (one of) the yield table's site index variable. For the Wiedemann pine table this is the quadratic mean height in meters. The function to be used is called ```site_index()```. Let's have a look at some examples: ```{r} ytab <- fe_ytable_pine_wiedemann_moderate_1943 # store the yield table in a # variable with a shorter name # for convenience # Assume a stand age of 73 years and a quadratic mean height of 19.7 m si <- site_index(age = 72, size = 19.7, ytable = ytab, si_variable = "h_q_m") si # The stand's relative site index according to the table si |> round(digits = 1) # usually relative site indexes are rounded to one digit # Same height, but twenty years younger site_index(age = 52, size = 19.7, ytable = ytab, si_variable = "h_q_m") |> round(digits = 1) # Same height, but thirty years older site_index(age = 102, size = 19.7, ytable = ytab, si_variable = "h_q_m") |> round(digits = 1) ``` Input values for stand ages that are not exactly given in the table, as in the examples above, are made use of by linear interpolation between the values given in the corresponding matrix. Site index values that are outside the site index fan, are obtained by linear extrapolation. This is the case in the second example, where we obtain a site index of 0.9 which reflects better site conditions than assumed by the best site index given in the table (i.e. 1.0). If a stand age provided by the user is beyond the age coverage of the yield table (as provided in the ```fe_yield_table``` object's slot ```age_coverage```), a warning is issued, and the age that is internally used is set to the nearest value covered in the table: ```{r} site_index(age = 171, size = 30.4, ytable = ytab, si_variable = "h_q_m") site_index(age = 12, size = 9.2, ytable = ytab, si_variable = "h_q_m") ``` Obviously, site indexing can be done on input vectors, in the following example we are doing this by way of the function ```map2_dbl``` from the *purrr* package: ```{r} age <- c(35, 25, 120, 75, 42, 53) h_q <- c(14.3, 9.1, 22.2, 13.6, 11.0, 21.7) si_vec <- purrr::map2_dbl( .x = age, .y = h_q, .f = function(x, y, yt, si_v) site_index(x, y, yt, si_v), yt = ytab, si_v = "h_q_m" ) |> round(digits = 1) si_vec ``` In order to convert this site index into a mean annual increment (mai) site index (German: "DGZ-Bonität"), there are two options. First, the conversion can be made to the expected mai at a given age, typically 100 years. This is achieved with the function ```si_to_mai_age```: ```{r} # Obtain the expected mean annual increment at age 100 for a stand with the # relative site index 1.2 si_to_mai_age( si = 1.2, mai_variable = "mai_m3_ha_yr", age = 100, ytable = ytab ) # Do the same for the site index vector we generated above si_vec purrr::map_dbl( .x = si_vec, .f = function(x, ma_v, age, yt) si_to_mai_age(x, ma_v, age, yt), ma_v = "mai_m3_ha_yr", age = 100, yt = ytab ) ``` The second option is not to determine the mai site index at a given age, but to take the maximum mai, that does typically occur not at the same age among different site qualities (earlier on good sites, later on poor sites). The function ```si_to_mai_max``` was written for that purpose: ```{r} # Obtain the expected maximum mean annual increment for a stand with the # relative site index 1.2 si_to_mai_max( si = 1.2, mai_variable = "mai_m3_ha_yr", ytable = ytab ) # Do the same for the site index vector we generated above si_vec purrr::map_dbl( .x = si_vec, .f = function(x, ma_v, yt) si_to_mai_max(x, ma_v, yt), ma_v = "mai_m3_ha_yr", yt = ytab ) ``` Note that in some cases, especially when extrapolating with site indexes that are above the best site covered in a yield table, the obtained mai max site index might be smaller than the mai site index obtained for an age that is near the maximum mai. Unfortunately, such subtle inconsistencies are in the very nature of yield tables in general, and cannot be avoided. In order to convert relative site indexes into absolute ones, e.g. the expected stand height at an age of 100 years, you can use the function ```ytable_lookup``` which is the general means to draw values out of yield tables, once the site index is known: ```{r} # Obtain the expected stand height at age 100 for a stand with the relative # site index 1.2 (this actually converts a relative into an absolute site index) si_abs <- ytable_lookup(age = 100, si = 1.2, variable = "h_q_m", ytable = ytab) si_abs si_abs |> round(digits = 1) # Let's to this again for the whole site index vector from above. # Reasonably, we obtain warnings, if we use site indexes beyond range of the # yield table si_abs_vec <- purrr::map_dbl( .x = si_vec, .f = function(x, age, var, yt) ytable_lookup(age, x, var, yt), age = 100, var = "h_q_m", yt = ytab ) si_abs_vec |> round(digits = 1) ``` ### 4.2 Extracting Values From Yield Tables The key function for extracting desired values from yield tables, once the site index is known, is ```ytable_lookup``` as already shown above when the task was to obtain absolute site indexes from relative ones. Let us give some more examples here: ```{r} # Get the periodic annual increment: ytable_lookup(age = 100, si = 1.2, variable = "pai_m3_ha_yr", ytable = ytab) # Try to go beyond the table's age coverage (raises warning) ytable_lookup(age = 170, si = 1.2, variable = "pai_m3_ha_yr", ytable = ytab) # ... standing volume ytable_lookup(age = 73, si = 2.4, variable = "v_m3_ha", ytable = ytab) # Use a site index above the table's coverage (raises warning) ytable_lookup(age = 73, si = 0.4, variable = "v_m3_ha", ytable = ytab) # ... basal area ytable_lookup(age = 41, si = 3.4, variable = "ba_m2_ha", ytable = ytab) # Use a site index below the table's coverage ytable_lookup(age = 41, si = 6.2, variable = "ba_m2_ha", ytable = ytab) ``` This way, any variable that is contained in the ```fe_yield_table``` object can be easily accessed. Vectorized application is available with standard R functionality (we prefer the functions from the *purrr* package, but the *apply* family of standard R works as well). ## 5 Importing Yield Table Data and Make them fe_yield_table Objects The key function for transoforming your own data into valid and fully functional ```fe_yield_table``` objects is ```fe_yield_table```. We will add some more explanations here; for the time being, be referred to the function's documentation. You can access it with ```?fe_yield_table```. ## References