The Package ForestElementsR
Peter Biber and Astor Toraño Caicoya
library(ForestElementsR) # Attach ForestElementsR
# Other packages required in the code examples below
library(dplyr)
library(ggplot2)
library(purrr)1 Introduction
Many methods from the fields of forest growth and yield, as well as forest inventory, that belong to the core contents of forest related education were not widely available in a contemporary form so far. To cover this deficit and to continue doing so in the long run is the driving idea behind the creation and development of the package ForestElementsR. The concept of the package is mainly based on experiences and procedures made and developed at the former Chair of Forest Growth and Yield at the Technical University of Munich. Therefore, the method and object collection covered by the current version of ForestElementsR is slightly biased towards German and Bavarian conditions. This is, however, not intended to be a limitation for the the package’s future development. For this reason - too much pre-definition leads to restriction - ForestElementsR does not provide readily defined workflows. It does, in contrast, offer a collection of elementary (hence the name of the package) building blocks that can be combined into any workflow required by users for their specific requirements. If you are a developer, feel encouraged to use the methods and objects provided by ForestElementsR in your own forest-related R packages, e.g. for evaluating specific types of forest inventories, research plots, or the output of forest growth models. This is exactly what we, the authors of ForestElementsR, are currently using the package for.
Most of the functionality we provide with the package is focused on producing information that is typically connected to the fields of forest growth and yield, forest mensuration, and forest inventory. Therefore, there is a certain bias towards goal variables which are related to the production of wood in one way or the other. However, even data from classic forest surveys can be evaluated in many ways that allow conclusions about forest structure and diversity. A few such methods are implemented in the current version of ForestElementsR, and we intend to add more in the future. We deem this important, because such approaches make a better use of forest data with regard to multifunctional forest management and planning.
It is an important design requirement of ForestElementsR not to be just a conglomerate of functions and objects that exist independently of each other. We hope users will see overarching concepts when they find similarities in the naming of our functions and objects, and the composition of function calls. Other, even more important bundling features are the package’s open Species Coding System and the hierarchically organized family of objects that represent (samples from) forest stands
With ForestElementsR we would like to address a broad audience. Practitioners might want to apply it for many purposes from using single functions in a “forest pocket calculator” style to evaluate forest stand surveys by arranging functions in a workchain way. Scientists will hopefully find the package useful for carrying out routine evaluations of forest data (like volume estimates, site indexing) that are regular prerequisites for in-dephth scientific analyzes. Students could profit from applying ForestElementR’s functionality to the extensive set of example data that comes with the package. We are confident that this would help them to develop an understanding and intuition about quantitative key methods in forestry as well as about the quantities themselves.
We should add that, besides our own developments, we also have quite shamelessly implemented methods developed and published by other authors that have proven themselves useful in our daily work. We have tried our best to find and cite the original publications correctly in order to give credit to whom it belongs.
Before we proceed, let us put out an Important Warning: While the functions and objects we provide with ForestElemensR behave exactly as described in the documentation, professional training in the fields of forestry, forest management, forest science is required to apply them correctly and to understand what their output actually means. If you are a hobby forester, you are more than welcome to make use of this package, but be sure to consult a professional before you draw important conclusions from what you get out of your evaluations. We also would like to point out another trivial wisdom that, however, seems to be ignored on a regular basis: If your data are bad, don’t expect to get good results.
In the following sections we give an overview about the most important
features of ForestElementsR. We do not delve deeply into all
details as these can be taken from the specific documentations provided
for each function. For two major features, namely the species coding system, and the
representation of yield tables, we
provide their own vignettes.
2 Object Families
ForestElementsR comprises three main families of S3 objects that are intended to facilitate the users’ as well as the developers’ lives. These relate to the representations of forest plots and stands, the species coding system, and ForestElementR’s yield table implementation.
2.1 Representations of Plots and Stands
The basic and most versatile class for representing plots or stands
is fe_stand. For how to make such an object from data call
the documentation ?fe_stand. As each tree in an fe_stand
object can have its own representation number (i.e. an information about
how many trees it represents per ha), it can represent (a cutout from) a
forest stand as well as any kind of sample with varying representation
like, e.g., angle count samples (where the representation number of a
tree depends on its dbh) or n-tree samples (where typically, the
outermost tree is only half represented). Also a concentric circle based
design can be represented; however, we provide specialized objects for
these. Let’s have a closer look at such an object:
# spruce_beech_1_fe_stand is one of the included example objects,
# for an overview, type "?example_data"
spruce_beech_1_fe_stand
#> $stand_id
#> [1] "spruce_beech_1"
#>
#> $area_ha
#> [1] 0.49
#>
#> $trees
#> # A tibble: 225 × 12
#> tree_id species_id layer_key time_yr age_yr dbh_cm height_m
#> <chr> <tm_wwk_shrt> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 1 1 2022 75 26.3 27.5
#> 2 2 1 1 2022 75 34.8 30.9
#> 3 3 1 1 2022 75 34.1 32.1
#> 4 4 1 1 2022 75 33.6 32.7
#> 5 5 1 1 2022 75 27.5 29.7
#> 6 6 1 1 2022 75 30 30
#> 7 7 1 1 2022 75 32.8 31.2
#> 8 8 1 1 2022 75 34.8 30.9
#> 9 9 1 1 2022 75 33.6 29.9
#> 10 10 1 1 2022 75 38.5 31.5
#> # ℹ 215 more rows
#> # ℹ 5 more variables: crown_base_height_m <dbl>, crown_radius_m <dbl>,
#> # removal <lgl>, ingrowth <lgl>, n_rep_ha <dbl>
#>
#> $small_trees
#> Dataframe mit 0 Spalten und 0 Zeilen
#>
#> attr(,"class")
#> [1] "fe_stand"An fe_stand object is evidently a list-like object; its most important element is the data frame (tibble) “trees”, which contains all trees that were measured individually and have a dbh. Note, that the element small_trees is still experimental and not yet used by any function we provide with ForestElementsR. It is intended to, in the future, provide a technical home for trees so small that they do not have a dbh (i.e. trees with heights < 1.3 m) or are under the size threshold for being measured individually, but are counted in size classes. Due to the multitude existing approaches, we haven’t decided on a representation in fe_stand objects, yet. The only current requirement for the small_trees element is that it has to be a data frame. In our example above, this data frame is empty. We have written fe_stand versions for the S3 generic functions summary and plot:
oo <- options(fe_spec_lang = "eng") # display species names in English
spruce_beech_1_fe_stand |> summary() # Give a summary of stand-level values
#> # A tibble: 2 × 9
#> # Groups: time_yr [1]
#> time_yr species_id stem_number_ha basal_area_m2_ha d_q_cm d_dom_cm h_q_m
#> <dbl> <tm_wwk_shrt> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 2022 Norway spruce 306. 25.9 32.8 38.3 30.8
#> 2 2022 European beech 153. 9.85 28.6 36.0 28.1
#> # ℹ 2 more variables: h_dom_m <dbl>, v_m3_ha <dbl>
spruce_beech_1_fe_stand |> plot() # Make an appropriate plot
options(oo) # Reset species name display optionIf spatial information about the location of the stand (sample)
and/or the single trees is available, the class
fe_stand_spatial can be used. It is a child of fe_stand.
All spatial functionality is taken from the package sf, therefore,
fe_stand_spatial accepts a broad range of coordinate definitions, from
simple cartesian coordinates to more or less all common geocoordinate
definitions. The output of summary and plot
for such an object is virtually the same as for its parent class
fe_stand:
oo <- options(fe_spec_lang = "sci") # display scientific species names
# Example: A mixed mountain forest plot in the Bavarian Alps;
# for an overview, type "?example_data"
mm_forest_1_fe_stand_spatial |> summary() # Give a summary of stand-level values
#> # A tibble: 22 × 9
#> # Groups: time_yr [5]
#> time_yr species_id stem_number_ha basal_area_m2_ha d_q_cm d_dom_cm
#> <int> <tm_wwk_lng> <dbl> <dbl> <dbl> <dbl>
#> 1 1975 Picea abies 143. 21.4 43.6 60.3
#> 2 1975 Abies alba 101. 6.16 27.8 42.5
#> 3 1975 Fagus sylvatica 161. 3.83 17.4 29.5
#> 4 1975 Acer pseudoplatanus 101. 3.55 21.1 27.2
#> 5 1984 Picea abies 71.6 14.1 50.2 64.0
#> 6 1984 Abies alba 53.7 3.47 28.7 46.6
#> 7 1984 Fagus sylvatica 89.5 3.22 21.4 36.7
#> 8 1984 Acer pseudoplatanus 59.7 2.54 23.3 31.4
#> 9 1995 Picea abies 71.6 15.8 52.9 67.6
#> 10 1995 Abies alba 41.8 3.87 34.4 50.2
#> # ℹ 12 more rows
#> # ℹ 3 more variables: h_q_m <dbl>, h_dom_m <dbl>, v_m3_ha <dbl>
mm_forest_1_fe_stand_spatial |> plot() # Make an appropriate plot
options(oo) # Reset species name display optionHowever, the spatial information about the outline of the stand (plot) and the positions of the trees can also be accessed graphically. In our example, the object represents a rectangular research plot where the coordinates are not geolocated:
mm_forest_1_fe_stand_spatial$outline |>
ggplot() +
geom_sf(fill = NA) +
geom_sf(data = mm_forest_1_fe_stand_spatial$tree_positions) +
xlab("x (m)") + ylab("y (m)")
An important special child of fe_stand_spatial is ´the
fe_ccircle_spatial object. It is a generic container for
sample plots with a concentric circle design. A multitude of such
objects can be used for representing entire sample forest inventories.
The latter is actually the intended use of
fe_ccircle_spatial which can be done e.g. in R packages
that build up on ForestElementsR. The summary and
plot functions are also available for these objects:
oo <- options(fe_spec_lang = "ger") # display German species names
spruce_pine_ccircle_spatial |> summary()
#> # A tibble: 2 × 9
#> # Groups: time_yr [1]
#> time_yr species_id stem_number_ha basal_area_m2_ha d_q_cm d_dom_cm h_q_m
#> <dbl> <tm_wwk_lng> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 2020 Fichte 207. 13.8 29.1 31.4 23.8
#> 2 2020 Kiefer 642. 32.0 25.2 43.8 33.3
#> # ℹ 2 more variables: h_dom_m <dbl>, v_m3_ha <dbl>
spruce_pine_ccircle_spatial |> plot(dbh_scale = 3) # oversize dbh x3
options(oo) # Reset species name display optionWe provide user friendly functions for generating all the objects
described above. These functions have the same name as the class of the
object to be created. E.g. with fe_stand() you can generate
an object of class fe_stand from your own data. This will
only be successful if the candidate object passes a validation that
makes sure it is not ill-defined. For details, check the documentations
of the functions fe_stand(),
fe_stand_spatial(), fe_ccircle_spatial(), and
fe_ccircle_spatial_notrees() (a class designed for special
cases e.g, permanent forest inventory samples where one plot does not
have any trees at a given inventory period).
2.2 Species Coding System
Unfortunately, the different forest research institutions in Central Europa are typically using different tree species codings, which are typically not or only partly compatible with each other. This creates considerable friction e.g. in cross sectional analyses of data that come from different institutions. Therefore, ForestElementsR comes with a generic open species coding system that allows to cast between different codings. As we provide its own vignette, which covers this species coding system in detail, we give only an outline here.
The system is designed in a way that coding casts work as long as there is no logical problem in the specific situation to solve, even when the two codings of interest are not entirely compatible. That means, e.g., when incompatibilities between two codings exist for very uncommon species or species groups only, they will not hamper the majority of casts. When a species coding cast loses information, a warning will be issued, but the cast will be executed. This is the case when a species cast leads from an individually coded species in the original coding into a species group code in the goal coding. Currently, six species codings are implemented in ForestElementsR. Five of them are codings that are in use in practice, probably most prominently the coding of the German National Forest Inventory (Riedel et al. 2017). All specific codings are linked to an internal master coding, which guarantees consistency. While the authors have a faint hope that this master coding could eventually become a common standard, they are also realists and understand the longevity of traditions.
Many functions in ForestElementsR that utilize species
codings, e.g. the volume function v_gri(), will try to
convert any given species coding into the most similar supported one for
which parameters are available. In this way, useful functions that were
originally developed at one institution and are therefore compatible
with a specific species coding, become more widely applicable without
forcing users to do tedious technical pull-ups beforehand. When
examining the code examples in the previous
section, you may have noticed that R’s options()
mechanism is used there to select the language which is used to display
species names. Let us demonstrate this here again:
# Define a vector of species ids compatible with the coding used by the
# Bavarian Forest Service
species_ids <- as_fe_species_bavrn_state(c(10, 20, 21, 24, 68, 62, 60))
# How they are displayed depends on the current setting of the option
# fe_spec_lang
species_ids
#> <fe_species_bavrn_state[7]>
#> [1] 10 20 21 24 68 62 60
# Check the option
options("fe_spec_lang") # if NULL or "code" the number codes are displayed
#> $fe_spec_lang
#> NULL
# Set to English (and keep the prevous option settings in oo)
oo <- options(fe_spec_lang = "eng")
species_ids
#> <fe_species_bavrn_state[7]>
#> [1] Norway spruce Scots pine eastern white pine Arolla pine
#> [5] sweet cherry small-leaved lime European beech
# Set to scientific names
options(fe_spec_lang = "sci")
species_ids
#> <fe_species_bavrn_state[7]>
#> [1] Picea abies Pinus sylvestris Pinus strobus Pinus cembra
#> [5] Prunus avium Tilia cordata Fagus sylvatica
# Set to German
options(fe_spec_lang = "ger")
species_ids
#> <fe_species_bavrn_state[7]>
#> [1] Fichte Kiefer Strobe Zirbe Vogelkirsche
#> [6] Winterlinde Buche
# Set to codes
options(fe_spec_lang = "code")
species_ids
#> <fe_species_bavrn_state[7]>
#> [1] 10 20 21 24 68 62 60
# Reset to the previous setting before enforcing English
options(oo)See the vignette Tree Species Codings in ForestElementsR for an in-detail description of the species coding system of ForestElementsR
2.3 Yield Tables
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. For about three decades,
however, they have been increasingly frowned upon for several reason,
including 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 since changed,
and cannot be considered as stable as they were 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. An initial 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. Currently, it covers the yield tables contained
in BayMinELF (1990). All readily available
tables are listed in the ForestElementsR’s documentation. You
can look them up with ?yield_tables. It is also possible
and intended that users can make their own preferred yield tables usable
with ForestElementsR even though they are not part of the
collection that comes with the package. The key object here is
fe_yield_table, and the function
fe_yield_table() was designed for generating such objects
in a user friendly way. See the vignette Yield Tables in ForestElementsR for a
detailed description of ForestElementR’s yield table
system.
The yield table system allows for an easy visualization of yield
tables. We demonstrate that by example of the widely known yield table
for Norway spruce by Wiedemann (1942).
Applying plot() to an object of class
fe_yield_table displays the ‘site index fan’, i.e. the mean
stand height over age curves for the site indexes (si) covered by the
table.
fe_ytable_spruce_wiedemann_moderate_1936_42 |> plot()
However, the visualization is not limited to the site index fan of a yield table. Actually, any variable that is covered by a given yield table can be visualized. See below how to obtain the names of all these variables:
# Get all variables available in the yield table of interest
vars <- names(fe_ytable_spruce_wiedemann_moderate_1936_42$values)
vars
#> [1] "h_q_m_si_plus_025" "n_ha"
#> [3] "ba_m2_ha" "d_q_cm"
#> [5] "v_m3_ha" "pai_m3_ha_yr"
#> [7] "pai_perc_yr" "mai_m3_ha_yr"
#> [9] "red_v_m3_ha" "red_pai_m3_ha_yr"
#> [11] "red_mai_m3_ha_yr" "red_pre_yield_m3_ha_10yr"
#> [13] "h_q_m" "tvp_m3_ha"The following code plots the age-curves of the standing volume
(variable name v_m3_ha); while the resulting figure is not
embedded in this vignette due to memory constraints, the reader is
encouraged to run this example.
# select standing volume ...
fe_ytable_spruce_wiedemann_moderate_1936_42 |> plot(variable = "v_m3_ha")With the code example below, we plot the periodic annual increment
(variable name pai_m3_ha_yr).
# ... and the periodic annual incement for plotting
fe_ytable_spruce_wiedemann_moderate_1936_42 |> plot(variable = "pai_m3_ha_yr")
Historical fact: As the diagram obove clearly shows, Wiedemann’s yield tables were not constructed by fitting functions to electronic data, but in a rather visual and manual way to point clouds on millimeter paper. At the time, viable alternatives were not available. Nevertheless, ForestElementsR keeps yield table values as provided (as long as they are consistent) and does not attempt to smooth or harmonize them. We deem that important to ensure consistency with previous applications of widespread yield tables. For more information about the yield table system of ForestElementsR see the vignette Yield Tables in ForestElementsR.
3 Low Level Methods
The functions we call low level methods are arguably the
most important constituents of ForestElementsR. These functions
do not directly work on objects of the fe_stand family, but
they are typically the elementary ingredients of such. The purpose of
these low level methods is to be as generic as possible. This makes them
the main tools for users who want to use ForestElementsR in a
“forestry pocket calculator” style. As their input, they typically
expect vectors of tree variables. Their output typically is a single
value that is an aggregate of the input version or a vector that
corresponds to each element of the input vector(s). Randomly chosen
examples for such functions are d_q (quadratic mean
diameter), d_100 (dominant diameter), h_q,
h_100 (quadratic mean height, dominant height), and
site_index (site index according to a yield table). If low
level methods require species information, which is usually the case if
they rely on species specific parameters, a vector of a supported
species coding (i.e. an object from the fe_species family) has
to be their first input variable. Examples for this are
v_gri (merchantable wood volume),
shannon_index (species diversity index after Shannon and Weaver (1948)),
h_standard_bv (Bavarian standard height curve system after
Kennel (1972)). We demonstrate a little
workchain in “forestry pocket calculator” style:
Assuming we have sampled eight trees in a forest stand and recorded their diameters at breast height, dbh, in cm.
dbh <- c(30.5, 30.9, 33.3, 35.1, 19.4, 19.5, 28.1, 25.5)
For making things a bit more complex, we also assume that these trees have been sampled with the angle count method and a factor of 4. This means that each tree represents a basal area of 4 m²/ha. Therefore, the number of trees that each tree represents per ha depends on its diameter at breast height:
# Calculate each tree's basal area in m²
ba <- (dbh / 100)^2 * pi/4
# Angle count factor 4/ba gives each tree's individual representation number
# per ha
nrep_ha <- 4 / ba
nrep_ha
#> [1] 54.74827 53.34002 45.92843 41.33861 135.32145 133.93710 64.49967
#> [8] 78.32308Given this information, we have functions for calculating classic stand diameters:
# Quadratic mean diameter
dq <- d_q(dbh, nrep_ha)
dq
#> [1] 25.8988
# Note that the quadratic mean diameter is always greater than the arithmetic
# mean as calculated here
stats::weighted.mean(dbh, w = nrep_ha)
#> [1] 25.26209
# Assmann's dominant diameter (quadratic mean diameter of the 100 thickest stems
# per ha)
d_100(dbh, nrep_ha)
#> [1] 33.76636
# Weise's dominant diameter (quadratic mean diameter of 20% thickest stems)
d_dom_weise(dbh, nrep_ha)
#> [1] 33.27737As a next step we would like to calculate the trees’ wood volumes. However, since their heights were not measured, we used a few orienting height measurements of trees with diameters close to the quadratic mean diameter dq of 25.9 cm. From this measurements we estimated the corresponding height, hq, to be 23.6 m. Let’s further assume that all trees in our sample are European beech (Fagus sylvatica) trees from a 72-year-old monospecific stand. This information allows us to estimate heights for all trees with the standard height curve system used by the Bavarian State Forest Service (Kennel 1972):
# Make a vector of 8 species ids = 5 and transform it into a tum_wwk_short
# species coding vector (see vignette about Species Codings in ForestElementsR),
# where code 5 stands for European beech.
# For our monospecific sample, we could make things even simpler (i.e. providing
# only one species id value for all trees together to the functions for height
# and volume estimates below), but we prefer to show the more generic approach
# here, where each tree has its own corresponding species id.
species_id <- rep(5, 8) |> as_fe_species_tum_wwk_short()
species_id
#> <fe_species_tum_wwk_short[8]>
#> [1] 5 5 5 5 5 5 5 5
oo <- options(fe_spec_lang = "eng") # display species names in English
species_id
#> <fe_species_tum_wwk_short[8]>
#> [1] European beech European beech European beech European beech European beech
#> [6] European beech European beech European beech
hq <- 23.6 # Height value corresponding to the quadratic mean diameter dq
age <- 72 # Stand age
# Estimate a height for each tree in the stand based on the information above
# with the Bavarian standard height curve system
h <- h_standard_bv(species_id, dbh, age, dq, hq)
h # Vector of height estimates corresponding to the dbh vector
#> [1] 24.64269 24.72119 25.15889 25.45388 21.49704 21.53731 24.13349 23.49549
options(oo) # Switch back options to previous valuesHaving now dbh and height values for all trees, we can calculate
their volumes. To that end, ForestElementsR provides the
function v_gri() that estimates a trees’ standing
merchantable wood volume over bark in m³:
v <- v_gri(species_id, dbh, h)
v
#> [1] 0.9195415 0.9479353 1.1274727 1.2724106 0.3071555 0.3111566 0.7583854
#> [8] 0.6016652Summing up the single tree volumes multiplied by the corresponding tree’s representation number per ha, we obtain the ha-related stand volume (as estimated from our sample of 11):
v_ha <- sum(v * nrep_ha)
v_ha
#> [1] 384.569For estimating the volume increment, we could use a yield table
(being aware that yield tables substantially underestimate the actual
growth observed in Central Europe over the last decades). For our stand
we will use the yield table for European beech by Wiedemann (1931) as
published in BayMinELF (2018). This table
is available as an fe_yield_table object in
ForestElementsR under the name
fe_ytable_beech_wiedemann_moderate_1931. In the following
code examples, we will not provide too many details about how to
actually use yield tables in ForestElementsR. See the vignette
Yield Tables in ForestElementsR for
detailed explanations. The first thing we have to do is to find out the
site index (ger. “Bonität”) of our stand according to the selected yield
table. For this, we require the stand’s quadratic mean height, hq = 23.6
m, and the stand’s age, which is 72 years:
si <- site_index(
age = age,
size = hq,
ytable = fe_ytable_beech_wiedemann_moderate_1931,
si_variable = "h_q_m"
)
si
#> [1] 1.518717This site index of 1.5 refers to the traditional way of site indexing
where 1.0 expresses the best site quality covered by a yield table.
Typically, 4.0, 5.0, or 6.0, depending on the table, would indicate the
weakest site conditions. In the vignette Yield Tables in ForestElementsR, we
describe alternative supported ways of expressing site indexes. With the
site index now known, the function ytable_lookup() can be
used to obtain the stand’s current increment in m³/ha/yr according to
the table.
inc <- ytable_lookup(age = age, si = si, variable = "pai_m3_ha_yr",
ytable = fe_ytable_beech_wiedemann_moderate_1931)
inc
#> [1] 10.84781This increment estimate, 10.8 m³/ha/yr, refers to a fully stocked stand. To find out whether it should be corrected due to the actual stand density, the stocking level (ger. “Bestockungsgrad”) of the stand of interest should be determined. We require the stand’s basal area, its age, and the site index for doing so. As we have assumed that our trees are an angle count sample with factor 4, each tree represents a basal area of 4 m²/ha. Therefore, the stand’s basal area amounts to 32 m²/ha. The stocking level results in:
# Basal area per hectare from angle count sample
ba_ha <- length(dbh) * 4
ba_ha
#> [1] 32
# Stocking level
stl <- stocking_level(ba = ba_ha, age = age, si = si,
ytable = fe_ytable_beech_wiedemann_moderate_1931)
stl
#> [1] 1.045007With a stocking level of about 1, our stand is fully stocked compared to the yield table. While some yield tables come with their own increment correction tables for deviating stocking levels, we have decided not to implement these in ForestElementsR so far. Instead, we advocate for the simple, classic way of multiplying the original increment estimate by the stocking level if the latter is less than one. An easy generic way for applying this rule in R would be as follows:
inc * min(1, stl) # correct the increment only if stocking level < 1
#> [1] 10.84781A simple alternative way for estimating a stand increment in our example situation could be applying the height and diameter growth functions of the German National Forest Inventory (Riedel et al. 2017) to the quadratic mean diameter tree (ger. “Grundflächenmittelstamm”):
# Current volume of the quadratic mean diameter tree (dbh = dq, height = hq)
species_id <- as_fe_species_tum_wwk_short(5) # Beech in tum_wwk_short coding
vol_current <- v_gri(species_id, dq, hq)
# future (+5 years) dq and hq estimated with German NFI3 functions
dq_future <- d_age_gnfi3(species_id, age + 5, dq, age)
hq_future <- h_age_gnfi3(species_id, age + 5, hq, age)
# future volume of the quadratic mean diameter tree (dq-tree)
vol_future <- v_gri(species_id, dq_future, hq_future)
# Estimated annual volume increment of the dq-tree
# Division by 5 years, because we want an annual value
iv_mean_tree <- (vol_future - vol_current) / 5
# The dq-tree currently represents so many trees per ha ...
n_ha <- sum(nrep_ha)
# ... therefore, we estimate the stand's increment as
iv_stand <- iv_mean_tree * n_ha
iv_stand
#> [1] 14.78204It does not come as a surprise that the latter increment estimate is higher than the one obtained from the yield table. The growth functions from the German National Forest Inventory (NFI) were calibrated with data that reflect the recent reality of forest growth, in contrast with Wiedemann’s yield table which was developed using data that reflect the growth of over nine decades or more ago. See Pretzsch et al. (2014), and Pretzsch et al. (2023) for more information about forest growth trends in Central Europe. In general, increment estimates that are not directly derived from repeated measurements must always be interpreted with care and professional knowledge. Both methods used in the code examples above, have their specific strengths and weaknesses. Do not use them if you are not aware of these.
We would like to draw your attention to another family of low level
functions, which annable the evaluating of forest structure and
diversity. In this field, we see a huge potential for assessing
ecosystem services beyond wood production and biodiversity based on
standard forest mensuration data (see Biber et
al. (2021), and Biber et al.
(2020)). Let us have a look at the function
shannon_index() that allows you to calculate the classic
species diversity index after Shannon and Weaver
(1948). In the simplest case, the function
shannon_index() requires just a vector of tree species
codes as an input. See the function’s documentation for more options. We
apply it to example data that come with ForestElementsR (type
?example_data for an overview), namely a monospecific
stand, an even-aged mixed stand of two species, and a selection forest
(ger. “Plenterwald”):
# Monospecific stand
# 1. Extract the tree data frame frame from an fe_stand_object
trees <- norway_spruce_1_fe_stand$trees
# 2. Run shannon_index() on its species_id column (technically a vector)
shannon_index(trees$species_id)
#> [1] 0
# Two-species mixed stand
trees <- spruce_beech_1_fe_stand$trees
shannon_index(trees$species_id)
#> [1] 0.6365142
# Selection forest
trees <- selection_forest_1_fe_stand$trees
shannon_index(trees$species_id)
#> [1] 1.149218Clearly, there is a plausible gradient with the monospecific stand
having zero species diversity, the selection forest having the highest
diversity, and the even-aged mixed stand being in between. The species
profile index by Pretzsch (2009) uses the
same basic concept, but adds a vertical component. Therefore, the
species profile index is higher, the more vertically structured a stand
is, and the more species it contains. In a mono-layered, monospecific
stand, this index is zero. In ForestElementsR, it is
implemented as the function species_profile, whose minimum
input requirements are a vector of species codes, and a corresponding
vector of tree heights. We apply it to the same example stands as
above:
# Monospecific stand
trees <- norway_spruce_1_fe_stand$trees
species_profile(trees$species_id, trees$height_m)
#> [1] 0.06521436
# Two-species mixed stand
trees <- spruce_beech_1_fe_stand$trees
species_profile(trees$species_id, trees$height_m)
#> [1] 0.8928083
# Selection forest
trees <- selection_forest_1_fe_stand$trees
species_profile(trees$species_id, trees$height_m)
#> [1] 1.988993Predictably, the gradient we obtain is similar to what we find for the shannon index above. However, the species profile index of the monospecific stand is somewhat greater than zero, indicating that the stand is at least, to some degree, vertically structured. The values we obtain for the other two forest stands are quite typical for the stand types they represent. A lot more could be said about indicators for stand structure and diversity, but this is beyond the scope of this vignette. Interested readers are referred to the publications we mentioned above. We see an increasing demand for such indicators in forest practice and science, and plan to extend the corresponding family of functions in ForestElementsR in the future.
4 High Level Methods
With the term high level methods we mean functions that
expect objects of the fe_stand class or its children as one
of their inputs. Internally and usually, these functions are tailored
combinations of low level methods,
offering users more aggregated workflows that integrate seamlessly with
fe_standand child objects. Still, we consider the currently
existing suite of high level methods far from complete, but useful
already. As a design requirement, all high level functions that operate
on the single trees of a fe_stand (or child) object have in
common the argument tree_filter. By default, this argument
is TRUE which means that the function will use all trees
stored in the trees data frame inside the
fe_stand object (see Section
2.1). However, any valid R expression that evaluates to a logical
value for each row of the trees data frame will work. This
provides a high level of adaptability to varying data evaluation
requirements. Let us demonstrate this using the function
stand_sums_static and the example object
mm_forest_1_fe_stand_spatial which is a cutout of a mixed
mountain forest stand that has been surveyed several times since 1975.
The function stand_sums_static provides an overview of the
most important static stand sum and mean values (from mean diameters to
the wood volume per ha) for each available survey of an
fe_stand object:
oo <- options(fe_spec_lang = "eng")
# Include all trees in evaluation with stand_sums_static
mm_forest_1_fe_stand_spatial |> stand_sums_static()
#> # A tibble: 22 × 9
#> # Groups: time_yr [5]
#> time_yr species_id stem_number_ha basal_area_m2_ha d_q_cm d_dom_cm h_q_m
#> <int> <tm_wwk_lng> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1975 Norway spruce 143. 21.4 43.6 60.3 31.2
#> 2 1975 silver fir 101. 6.16 27.8 42.5 23.9
#> 3 1975 European beech 161. 3.83 17.4 29.5 19.8
#> 4 1975 sycamore maple 101. 3.55 21.1 27.2 19.6
#> 5 1984 Norway spruce 71.6 14.1 50.2 64.0 34.3
#> 6 1984 silver fir 53.7 3.47 28.7 46.6 26.4
#> 7 1984 European beech 89.5 3.22 21.4 36.7 24.4
#> 8 1984 sycamore maple 59.7 2.54 23.3 31.4 22.8
#> 9 1995 Norway spruce 71.6 15.8 52.9 67.6 33.6
#> 10 1995 silver fir 41.8 3.87 34.4 50.2 26.7
#> # ℹ 12 more rows
#> # ℹ 2 more variables: h_dom_m <dbl>, v_m3_ha <dbl>
# For applying a filter, we must refer to the column names of the "trees" data
# frame inside the fe_stand object. These are
mm_forest_1_fe_stand_spatial$trees |> names()
#> [1] "tree_id" "species_id" "layer_key"
#> [4] "time_yr" "age_yr" "dbh_cm"
#> [7] "height_m" "crown_base_height_m" "crown_radius_m"
#> [10] "removal" "ingrowth" "n_rep_ha"
# Assume, we want to include only the trees that were not removed ...
mm_forest_1_fe_stand_spatial |> stand_sums_static(tree_filter = !removal)
#> # A tibble: 22 × 9
#> # Groups: time_yr [5]
#> time_yr species_id stem_number_ha basal_area_m2_ha d_q_cm d_dom_cm h_q_m
#> <int> <tm_wwk_lng> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1975 Norway spruce 71.6 12.6 47.3 60.9 31.8
#> 2 1975 silver fir 53.7 2.93 26.4 44.0 24.4
#> 3 1975 European beech 89.5 2.31 18.1 32.1 21.0
#> 4 1975 sycamore maple 59.7 2.22 21.8 28.5 20.0
#> 5 1984 Norway spruce 71.6 14.1 50.2 64.0 34.3
#> 6 1984 silver fir 41.8 3.29 31.7 47.5 27.2
#> 7 1984 European beech 89.5 3.22 21.4 36.7 24.4
#> 8 1984 sycamore maple 53.7 2.40 23.9 31.6 23.1
#> 9 1995 Norway spruce 71.6 15.8 52.9 67.6 33.6
#> 10 1995 silver fir 41.8 3.87 34.4 50.2 26.7
#> # ℹ 12 more rows
#> # ℹ 2 more variables: h_dom_m <dbl>, v_m3_ha <dbl>
# ... and now only the removal trees
mm_forest_1_fe_stand_spatial |> stand_sums_static(tree_filter = removal)
#> # A tibble: 10 × 9
#> # Groups: time_yr [4]
#> time_yr species_id stem_number_ha basal_area_m2_ha d_q_cm d_dom_cm h_q_m
#> <int> <tm_wwk_lng> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1975 Norway spruce 71.6 8.86 39.7 58.4 30.5
#> 2 1975 silver fir 47.7 3.23 29.4 40.2 23.4
#> 3 1975 European beech 71.6 1.52 16.4 25.7 17.9
#> 4 1975 sycamore maple 41.8 1.33 20.2 24.9 19.0
#> 5 1984 silver fir 11.9 0.179 13.8 17.3 12.5
#> 6 1984 sycamore maple 5.97 0.132 16.8 16.8 18.5
#> 7 2004 Norway spruce 41.8 10.1 55.6 66.9 34.6
#> 8 2004 European beech 35.8 1.72 24.7 42.5 24.7
#> 9 2004 sycamore maple 11.9 0.440 21.7 24.2 19.7
#> 10 2015 Norway spruce 5.97 0.866 43 43 30.7
#> # ℹ 2 more variables: h_dom_m <dbl>, v_m3_ha <dbl>
# Focus on Norway spruce and the surveys after 2000 only. Note, for filtering
# on species names you must use the format() function
# 1. Remaining trees only
mm_forest_1_fe_stand_spatial |>
stand_sums_static(
tree_filter =
format(species_id, "eng") == "Norway spruce" & time_yr > 2000 & !removal
)
#> # A tibble: 2 × 9
#> # Groups: time_yr [2]
#> time_yr species_id stem_number_ha basal_area_m2_ha d_q_cm d_dom_cm h_q_m
#> <int> <tm_wwk_lng> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 2004 Norway spruce 35.8 6.28 47.3 70.1 34.0
#> 2 2015 Norway spruce 35.8 5.78 45.4 72.5 33.9
#> # ℹ 2 more variables: h_dom_m <dbl>, v_m3_ha <dbl>
# 2. Removal trees only
mm_forest_1_fe_stand_spatial |>
stand_sums_static(
tree_filter =
format(species_id, "eng") == "Norway spruce" & time_yr > 2000 & removal
)
#> # A tibble: 2 × 9
#> # Groups: time_yr [2]
#> time_yr species_id stem_number_ha basal_area_m2_ha d_q_cm d_dom_cm h_q_m
#> <int> <tm_wwk_lng> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 2004 Norway spruce 41.8 10.1 55.6 66.9 34.6
#> 2 2015 Norway spruce 5.97 0.866 43 43 30.7
#> # ℹ 2 more variables: h_dom_m <dbl>, v_m3_ha <dbl>
options(oo)When an fe_stand or child object comprises more than one
survey, stand level periodic annual increments can be calculated. The
high level function for this task is stand_sums_dynamic.
Internally, it uses the high level function
stand_sums_static we introduced above, and the low level
function stand_level_increment. As all trees represented in
an fe_stand object must have a dbh entry, the basal area
increment will always be calculated. If enough information is available,
i.e. each tree also has a height entry (be it measured or estimated),
the volume increment will be calculated, too. The method of increment
calculation follows the classic approach \(iv
= (vs_2 + vr_2 - vs_1)/p\), where \(iv\) is the mean annual stand volume
increment during a period of \(p\)
years between two subsequent surveys. With \(vs_2\) and \(vs_1\) we indicate the standing volume at
the second and the first of two subsequent surveys, respectively. \(vr_2\) is the volume that has been removed
(or died) after the first survey up to, and including, the second survey
#2. For the basal area increment, the approach is exactly the same as
shown for the volume increment. We demonstrate
stand_sums_dynamic by example of the mixed mountain forest
stand introduced above:
oo <- options(fe_spec_lang = "eng")
# Calculate basal area and volume increments for an fe_stand object that covers
# more than one survey
ssd <- stand_sums_dynamic(mm_forest_1_fe_stand_spatial)
#> Joining with `by = join_by(species_id, time_yr)`
#> Joining with `by = join_by(time_yr, species_id)`
#> Joining with `by = join_by(time_yr, species_id)`
ssd |> print(n = Inf)
#> # A tibble: 22 × 4
#> time_yr species_id iba_m2_ha_yr iv_m3_ha_yr
#> <int> <tm_wwk_lng> <dbl> <dbl>
#> 1 1975 Norway spruce NA NA
#> 2 1984 Norway spruce 0.176 3.94
#> 3 1995 Norway spruce 0.146 1.52
#> 4 2004 Norway spruce 0.0745 1.58
#> 5 2015 Norway spruce 0.0337 0.311
#> 6 1975 silver fir NA NA
#> 7 1984 silver fir 0.0598 1.10
#> 8 1995 silver fir 0.0532 0.574
#> 9 2004 silver fir 0.0825 1.49
#> 10 2015 silver fir 0.0981 1.59
#> 11 1975 European beech NA NA
#> 12 1984 European beech 0.101 1.72
#> 13 1995 European beech 0.107 1.33
#> 14 2004 European beech 0.0964 1.32
#> 15 2015 European beech 0.170 2.13
#> 16 1975 sycamore maple NA NA
#> 17 1984 sycamore maple 0.0354 0.794
#> 18 1995 sycamore maple 0.0289 0.0709
#> 19 2004 sycamore maple 0.0369 0.671
#> 20 2015 sycamore maple 0.0577 0.462
#> 21 2004 European ash NA NA
#> 22 2015 European ash 0.0399 0.175
options(oo)As in this example height values are available for all trees in all
surveys, not only the basal area increment, but also the volume
increment is calculated. The increment entry at the first survey is
NA, because there is no preceding survey, and an increment
entry always relates to the period between the survey of interest and
the previous one. Species overarching increments can be easily
calculated from the output of stand_sums_dynamic just by
standard R methods, or even more conveniently with group_by
and summarise as provided by the dplyr
package:
# Species-overarching increments
stand_sums_dynamic(mm_forest_1_fe_stand_spatial) |>
group_by(time_yr) |>
summarise(
iba_m2_ha_yr = sum(iba_m2_ha_yr, na.rm = TRUE),
iv_m3_ha_yr = sum(iv_m3_ha_yr, na.rm = TRUE)
)
#> Joining with `by = join_by(species_id, time_yr)`
#> Joining with `by = join_by(time_yr, species_id)`
#> Joining with `by = join_by(time_yr, species_id)`
#> # A tibble: 5 × 3
#> time_yr iba_m2_ha_yr iv_m3_ha_yr
#> <int> <dbl> <dbl>
#> 1 1975 0 0
#> 2 1984 0.372 7.56
#> 3 1995 0.335 3.49
#> 4 2004 0.290 5.06
#> 5 2015 0.399 4.67
# Zero increments for 1975 are obtained, because no previous survey is
# available, so no meaningful increment can be calculated there
Arguably, stand_sums_static and
stand_sums_dynamic are the most important high level
functions currently implemented in ForestElementsR. Other
prominent ones are survey_overview, and
species_shares. We also provide specific implementations of
the generic functions plot and summary for
fe_stand and child classes; the summary method
being basically a wrapper around stand_sums_dynamic. All
the high level functions mentioned so far have in common the
tree_filter argument (see above) that allows users to
isolate specific sub cohorts of trees for evaluation. A high level
function that does not require such a filter is
get_area_ha. For any fe_stand or child object
where this is meaningful, it returns the total area (in ha) covered by
the stand or plot in the terrain:
# fe_stand
get_area_ha(norway_spruce_1_fe_stand)
#> [1] 0.49
# fe_stand_spatial
get_area_ha(mm_forest_1_fe_stand_spatial)
#> [1] 0.16761
# fe_ccircle_spatial (the area of the outermost circle is returned)
get_area_ha(spruce_pine_ccircle_spatial)
#> [1] 0.05So far, with this vignette we have provided an overview of the features of the package ForestElementsR we deem most important for standard users. Such a vignette, however, is not the place to describe each function in detail. You find this information in the specific documentation of the functions.
5 Expert Methods
What we call expert methods are functions that are not intended to be
used by standard users who want to directly work with forest data, but
by developers. An important group of such functions are the constructors
for the S3 classes (i.e. the object
families) we implemented, e.g. fe_stand,
fe_yield_table, etc. Following standard practice, the names
of these constructors begin with “new_” followed by the name of the
class, e.g. new_fe_stand, or
new_fe_yield_table. Only persons who are well aware of what
they are doing should use these constructors, as they do not prevent
users from creating ill-formed objects. Checking objects for adherence
to the classes’ requirements is the task of the validators, another
important group of expert functions. Their names begin with “validate_”
followed by the name of the class, e.g. validate_fe_stand,
or validate_fe_yield_table. If you are using the
constructors for creating objects, it will most probably be an excellent
idea to use a validator afterwards. As a standard user you should always
use the helper functions that simply have the same name as the class of
interest for creating an object (e.g. the functions
fe_stand, fe_yield_table). Internally, these
apply the appropriate constructor and the validator in the correct way
and try to be as helpful to users as possible.
6 Technical Specifics
ForestElementsR uses the S3 object system which is arguably the most widely used of R’s approaches to object oriented programming (Wickham 2019). We deemed this approach ideal for using the advantages of object oriented programming in a lean pragmatic way.
7 Acknowledgments
The authors would like to thank the Bavarian Ministry for Nutrition, Agriculture, Forestry, and Tourism for funding the projects FeNEU: Ein innovatives Instrument für die forstliche Planung in Bayern (E062) and Ertragskundliche Betreuung der langfristigen Versuche (W007).