European Forestry Dynamics Model (EFDM) is a model for large scale scenario analysis. The model is an area-based matrix model where forest area is stratified with respect to, for example, use, management strategy and development (growth). The model can consider also other land uses, see Scenario 3. In a model run forest area moves in ‘a matrix’ span with the stratifying factors.
The efdm
package includes functions to
and an example dataset for illustration purposes.
In this vignette, we demonstrate ‘efdm’ in running a scenario with the example dataset. Typical workflow in a EFDM scenario analysis starts with preparing the datasets required as an input. Preparing the datasets is usually the most time consuming part of the scenario analysis, and the work flow depends on the data available. Therefore, the preparation is not explained in this vignette. However, the attached datasets offer an example for preparing own datasets.
The efdm
is loaded with
library(efdm)
The following packages are used in this vignette to manipulate data and plot results. However, the data manipulation and reporting the results could be done using other packages or even outside the R environment.
library(dplyr)
library(ggplot2)
library(sf)
In this vignette we are using example
, a collection of example datasets to describe working with efdm
. The dataset covers entire Finland and the initial state represents year 2016. However, the dataset is a toy data which is meant for learning the model concept.
data(example)
In this example we have chosen
In EFDM forest management is represented by activities and activity probabilities. In the first example we have three activities
The activity probabilities (actprob
below) are the probabilities for applying a specific management activity in a specific forest state described by the stratification variables. The activity probabilities are used to allocate forest area to different activities. A valid activity probability specifies activity probabilities for all possible states and activities. All probabilities should be positive and the sum of activity probabilities for each state should be 1.
<- example$actprob
actprob
actprob#> region soil sp vol age noman thin ff
#> 1 South Mineral other 1 1 0.976 0.023 0.001
#> 2 South Mineral other 1 2 0.981 0.019 0.000
#> 3 South Mineral other 1 3 0.994 0.006 0.000
#> 4 South Mineral other 1 4 0.991 0.008 0.001
#> 5 South Mineral other 1 5 0.764 0.024 0.212
#> 6 South Mineral other 1 6 0.974 0.025 0.001
#> [ reached 'max' / getOption("max.print") -- omitted 6294 rows ]
Each activity has a transition probability matrix. The transition probability describes movements of forest areas by the simulation time steps. It is used to move a forest state to the next states. A valid transition probability should specify transition probability for each state where the activity can be used. For example if final felling only applies to forests with some minimum volume, then there is no need to specify transitions from states with less than the minimum volume.
The final felling activity is represented by a transition matrix that has probability 1 to move any forest to a state where volume and age are in the smallest class.
The target state is independent of the starting state. It will always be vol=1
and age=1
, the smallest class in the stratification.
<- expand.grid(vol0=1:15, age0=1:35, vol1=1, age1=1, prob=1)
transprobs_ff
transprobs_ff#> vol0 age0 vol1 age1 prob
#> 1 1 1 1 1 1
#> 2 2 1 1 1 1
#> 3 3 1 1 1 1
#> 4 4 1 1 1 1
#> 5 5 1 1 1 1
#> 6 6 1 1 1 1
#> 7 7 1 1 1 1
#> 8 8 1 1 1 1
#> 9 9 1 1 1 1
#> 10 10 1 1 1 1
#> [ reached 'max' / getOption("max.print") -- omitted 515 rows ]
In efdm
the final felling activity is defined by giving
<- define_activity("ff", c("vol", "age"), transprobs_ff) ff
If there are no management activities or treatments in the forest, “no management” activity is used to make the forest grow. In this example, we assume that growth changes the values of volume and age and is dependent on region, soil type and tree species. Moreover, we assume that the growth is so different in between regions and tree species, that it should be estimates completely separately for each class.
To estimate the transition matrix, we use “pair data” obtained from two consecutive measurements of permanent sample plots in the national forest inventory. For estimatetransprobs
, the estimation function implemented in efdm
, the pair data is a data.frame
of pairs of dynamic variables and possible stratification variables, which affect the growth but are not affected by growth.
$noman_pairs
example#> region soil sp vol0 vol1 age0 age1
#> 1 South Mineral spruce 10 7 13 14
#> 2 South Mineral spruce 8 8 20 21
#> 3 South Mineral other 9 10 5 6
#> 4 South Mineral other 1 1 1 2
#> 5 South Mineral other 7 8 2 3
#> 6 South Mineral other 7 10 5 6
#> 7 South Mineral other 3 5 31 32
#> [ reached 'max' / getOption("max.print") -- omitted 15535 rows ]
In order to use estimatetransprobs
we still need two parts, a state space and a prior. The state space and prior are required because the pair data might not have observations for all possible states, but the transition matrix should be available for all states.
The state space is the collection of all possible states of stratifying variables. Since we have an activity probability for each state we can use the actprob
to obtain the statespace
<- actprob %>% select(c(region, soil, sp, vol, age))
statespace
statespace#> region soil sp vol age
#> 1 South Mineral other 1 1
#> 2 South Mineral other 1 2
#> 3 South Mineral other 1 3
#> 4 South Mineral other 1 4
#> 5 South Mineral other 1 5
#> 6 South Mineral other 1 6
#> 7 South Mineral other 1 7
#> 8 South Mineral other 1 8
#> 9 South Mineral other 1 9
#> 10 South Mineral other 1 10
#> [ reached 'max' / getOption("max.print") -- omitted 6290 rows ]
In this example dataset the number of states in the state space is the product of the number of classes in each stratifying variable, that is \(3*2*2*15*35 = 6300\).
The prior in estimatetransprobs
works by adding one observation to the pair data for each starting state. The efdm
has a few standard prior choices.
prior_grow(variable)
grows variable
by one class."nochange"
, forest state is not changed"uninformative"
, an observation with tiny weight is added to every possible target state.For the no management activity we are using prior_grow("age")
stating that there is always at least one observation where the age grows by one class but the volume is not changing.
Now we are ready to estimate the transition probabilities for no management
<- estimatetransprobs(c("vol", "age"), # the dynamic variables of the activity
transprobs_noman $noman_pairs, # pair data for no management
examplestatespace=statespace,
factors=c("soil"), # Information for different soil types
# is used with smaller weight
by=c("region","sp"), # Separate estimation for each
# region and species
prior=prior_grow("age"))
and to define the activity noman
<- define_activity("noman", c("vol", "age"), transprobs_noman) noman
The resulting transition probabilities look like
%>%
transprobs_noman arrange(vol0, age0, soil, region, sp, vol1, age1)
#> vol0 age0 soil region sp vol1 age1 prob
#> 1 1 1 Mineral Middle other 1 2 0.629107981
#> 2 1 1 Mineral Middle other 2 2 0.145539906
#> 3 1 1 Mineral Middle other 3 2 0.122065728
#> 4 1 1 Mineral Middle other 4 2 0.089201878
#> 5 1 1 Mineral Middle other 5 2 0.009389671
#> 6 1 1 Mineral Middle other 6 2 0.004694836
#> [ reached 'max' / getOption("max.print") -- omitted 13180 rows ]
We use exactly the same procedure for thinning that we used for no management. The only differences are the activity name and the pair data used.
<- estimatetransprobs(c("vol", "age"),
transprobs_thin $thin_pairs, # pair data for thinning
examplestatespace=statespace,
factors=c("soil"),
by=c("region","sp"),
prior=prior_grow("age"))
<- define_activity("thin", c("vol", "age"), transprobs_thin) thin
EFDM is an area based model. The initial state of the scenario is a data.frame
where for each state there is an area
(zeros may be omitted)
<- example$initial_state
state0
state0#> region soil sp vol age area
#> 1 South Mineral other 1 1 73518.749
#> 2 Middle Mineral other 1 1 105532.353
#> 3 North Mineral other 1 1 49583.432
#> 4 South Peat other 1 1 15035.772
#> 5 Middle Peat other 1 1 33996.388
#> 6 North Peat other 1 1 2039.203
#> 7 South Mineral spruce 1 1 186533.582
#> 8 Middle Mineral spruce 1 1 63984.635
#> [ reached 'max' / getOption("max.print") -- omitted 3271 rows ]
Next, we run the EFDM scenario for 20 time steps (100 years) starting with initial state state0
, activity probabilities actprob
and a list of activities.
<- list(noman, thin, ff)
activities <- runEFDM(state0, actprob, activities, 20) states1
runEFDM
produces a data.frame
of areas allocated to each activity at each time step.
%>%
states1 arrange(soil, region, sp, vol, age, time, activity)
#> soil region sp vol age area activity time
#> 1 Mineral Middle other 1 1 105.53235 ff 0
#> 2 Mineral Middle other 1 1 101627.65637 noman 0
#> 3 Mineral Middle other 1 1 3799.16472 thin 0
#> 4 Mineral Middle other 1 1 99.11577 ff 1
#> 5 Mineral Middle other 1 1 95448.48821 noman 1
#> 6 Mineral Middle other 1 1 3568.16778 thin 1
#> [ reached 'max' / getOption("max.print") -- omitted 209007 rows ]
Using so called result coefficients is a way to obtain results with respect to other forest properties than area. In this scenario we have growing stock volume, drain and harvest income as result variables.
The volume coefficients convert volume class vol
and dominant species sp
into growing stock volume volume
(m3/ha) in four tree species groups (species
): pine, spruce, broadleaves and all. The volume coefficients were estimated for the classes based on the species composition in the forest inventory data as class averages.
$vol_coef %>%
examplearrange(vol, sp, species)
#> vol sp species volume
#> 1 1 other all 0.50632300
#> 2 1 other broadleaf 0.12501893
#> 3 1 other pine 0.30152225
#> 4 1 other spruce 0.07978181
#> 5 1 spruce all 1.20000000
#> 6 1 spruce broadleaf 0.37598659
#> 7 1 spruce pine 0.09533751
#> 8 1 spruce spruce 0.72867590
#> 9 2 other all 3.39429400
#> 10 2 other broadleaf 0.74911648
#> 11 2 other pine 2.16464343
#> 12 2 other spruce 0.48053409
#> [ reached 'max' / getOption("max.print") -- omitted 108 rows ]
Drain is the harvest accumulation which is linked to management activities (activity), and separated by the volume class (vol) and dominant species (sp). It was also estimated from the observed changes in the forest inventory data as class averages. The drain is given as m3 per ha (drain) and timber assortment (assort): pulp/saw wood.
$drain_coef %>%
examplearrange(-vol, sp, assort, activity)
#> vol sp drain assort activity
#> 1 15 other 203.57043 pulp ff
#> 2 15 other 61.07113 pulp thin
#> 3 15 other 197.54969 saw ff
#> 4 15 other 59.26491 saw thin
#> 5 15 spruce 165.51780 pulp ff
#> 6 15 spruce 49.65534 pulp thin
#> 7 15 spruce 318.13899 saw ff
#> 8 15 spruce 95.44170 saw thin
#> 9 14 other 172.94581 pulp ff
#> 10 14 other 51.88374 pulp thin
#> [ reached 'max' / getOption("max.print") -- omitted 110 rows ]
Income is loosely based on the actual statistics of Finnish timber assortment prices over last years in unit eur/m3.
$income_coef
example#> assort activity euro
#> 1 saw ff 56
#> 2 saw thin 44
#> 3 pulp ff 20
#> 4 pulp thin 15
Simulation timesteps are mutated to mid-years of simulation steps (the simulation begins from year 2016):
<- states1 %>% mutate(time = factor(2016 + time*5)) states1
Growing stock volume is estimated by merging the area distribution by simulation steps (states1) and volume coefficients. The total growing stock volume (m3) is a result of multiplication of area (ha) and coefficient (m3/ha). Finally the result is visualized in a figure.
<- merge(states1, example$vol_coef) %>%
volume mutate (volume=area*volume) %>% #area is multiplied with m3/ha volume the get total volume in m3
filter(species != "all")
ggplot(volume) +
scale_fill_viridis_d(end = 0.9) +
geom_bar(aes(x=time, weight=volume/1000000000, fill=species)) +
labs(y=NULL,title=expression(paste("Growing stock, bil.",m^3)), x="Year", fill="") +
theme(legend.position = "bottom", axis.text.x = element_text(angle = -90))
Estimating age distribution does not require result coefficients. However, the age classes are defined for 50-year classes instead of original 5-year classes, and then plotted for only three time steps: years 2016, 2066 and 2116.
$ageclass <- cut(states1$age, breaks=c(0,10,20,30,35), include.lowest = TRUE,
states1#labels=c("0-50","51-100","101-150","150+"))
labels=c("-50","-100","-150","150+"))
$region <- factor(states1$region, labels = c("South", "Middle", "North"))
states1ggplot(subset(states1, time %in% c(2016,2066,2116))) +
geom_bar(aes(x=ageclass, weight=area/1000000, fill=region)) +
scale_fill_viridis_d(end = 0.9) +
facet_grid(cols=vars(time)) +
labs(y=NULL, title="Area, mill.ha", x="Ageclass", fill=NULL) +
theme(legend.position = "bottom", axis.text.x = element_text(angle = -90))
Income is estimated by first converting the drain (m3/ha) into income (eur/ha). Then multiplication with area (ha) gives the euros.
<- merge(example$drain_coef, example$income_coef) %>% mutate(euro = euro*drain)
euro <- merge(states1, euro) %>% mutate(income=euro*area)
removal
ggplot(subset(removal,!time %in% c(2116))) +
geom_bar(aes(x=time, weight=income/5000000000, fill=assort)) +
scale_fill_viridis_d(end = 0.9) +
labs(y=NULL,title = "Income, bil.€/year", x="5-year intervals", fill=NULL ) +
theme(legend.position = "bottom", axis.text.x = element_text(angle = -90))
In this example the tree species changes after final felling. Therefore we redefine final felling activity taking into account in addition to volume and age also the dominant species. The change probability depends on region and dominant species. Volume and age act as before. They move to the smallest classes (vol1=1 and age1=1).
<- statespace %>%
transprobs_ff_age_species select(vol0=vol, age0=age, sp0=sp, region) %>%
%>%
unique group_by(region, sp0, vol0, age0) %>%
summarize(data.frame(vol1=1, age1=1, sp1=c('other','spruce'),
prob=case_when(
=='other' && region=='South' ~ c(1, 0),
sp0=='other' && region=='Middle' ~ c(1, 0),
sp0=='other' && region=='North' ~ c(0.8, 0.2), #in other dominated forest
sp0# in northern boreal vegetation zone 20% of final felled area
# change to spruce dominated forest
=='spruce' && region=='South' ~ c(0.3, 0.7),
sp0=='spruce' && region=='Middle' ~ c(0.2, 0.8),
sp0=='spruce' && region=='North' ~ c(0, 1))))
sp0#> `summarise()` has grouped output by 'region', 'sp0', 'vol0', 'age0'. You can override using the `.groups` argument.
<- define_activity("ff", c("vol", "age", "sp"), transprobs_ff_age_species) ff_age_species
The list of activities includes the same noman and thin as in Scenario 1, and the above defined final felling. The same initial state as in Scenario 1 is run with the new list of activities and stored in variable states2.
<- list(noman, thin, ff_age_species)
activities2 <- runEFDM(state0, actprob, activities2, 20) states2
The proportion of spruce dominated forest area is estimated for time steps 0, 10 and 20 (years 2016, 2066 and 2116). And the result is presented by the vegetation zones as a map.
# Compute proportion of spruces in each time and region
<- states2 %>%
prop_spruces group_by(region, time) %>%
summarise(proportion = 100*sum((sp=="spruce")*area)/sum(area)) %>%
filter(time %in% c(0, 10, 20)) %>%
mutate(time = 2016 + 5*time)
#> `summarise()` has grouped output by 'region'. You can override using the `.groups` argument.
# Add bio-geographical regions (MetsaKasvVyoh provided in efdm package) to the data
<- merge(MetsaKasvVyoh, prop_spruces)
prop_spruces %>% ggplot() +
prop_spruces geom_sf(aes(fill=proportion)) +
facet_wrap(vars(time)) +
scale_fill_gradient(low = "white", high = "forestgreen") +
labs(fill="Proportion" ,title="Proportion of spruce dominated forests, %") +
guides(fill = guide_colourbar(frame.colour = "grey20",ticks.colour = "grey20")) +
theme(legend.position = "bottom", axis.text.x = element_text(angle = -90))
The map provided by: Finnish Environment Institute
First we add two land use classes to the statespace. Agriculture is stratified according to region and soil type, while other land use is only stratified according to region.
<- statespace
statespace3 $landuse <- "forest"
statespace3<- unique(statespace3 %>% mutate(sp=0, vol=0, age=0, landuse="agriculture"))
agriculture <- unique(statespace3 %>% mutate(soil=0, sp=0, vol=0, age=0, landuse="other"))
other <- rbind(statespace3, agriculture, other) statespace3
We use separate activities for deforestation to each land use class. Variables (vol, age, sp) not used by the agriculture are set to 0. Soil type and region are not changing as a result of deforestation to agriculture.
<- statespace3 %>%
transprobs_defor_to_agri filter(landuse=="forest") %>%
select(vol0=vol, age0=age, sp0=sp,landuse0=landuse) %>%
%>%
unique mutate(vol1=0, age1=0, sp1=0, landuse1="agriculture", prob=1)
<- define_activity("defor_to_agri",
defor_to_agri c("vol", "age", "sp", "landuse"),
transprobs_defor_to_agri)
Deforestation to other land use also changes soil type to 0.
<- statespace3 %>%
transprobs_defor_to_other filter(landuse=="forest") %>%
select(soil0=soil, vol0=vol, age0=age, sp0=sp,landuse0=landuse) %>%
%>%
unique mutate(soil1=0, vol1=0, age1=0, sp1=0, landuse1="other", prob=1)
<- define_activity("defor_to_other",
defor_to_other c("soil", "vol", "age", "sp", "landuse"),
transprobs_defor_to_other)
Afforestation only applies to agriculture. Volume and age classes start from 1 and the area is split evenly to spruce and other species.
<- statespace3 %>%
transprobs_aff filter(landuse=="agriculture") %>%
select(soil, vol0=vol, age0=age, region, sp0=sp, landuse0=landuse) %>%
mutate(vol1=1, age1=1, sp1="other", landuse1="forest", prob=1)
<- rbind(transprobs_aff %>% mutate(sp1="spruce", prob=0.5),
transprobs_aff %>% mutate(sp1="other", prob=0.5))
transprobs_aff <- define_activity("affor",
affor c("vol", "age", "sp", "landuse"),
transprobs_aff)
A donothing activity is used for non-forest land uses, when there is nothing forestry related going on.
<- define_activity("donothing", character())
donothing <- list(noman, thin, ff, defor_to_other, defor_to_agri, affor, donothing) activities3
adding land use information to state
<- state0
state03 $landuse <- "forest"
state03<- rbind(state03, agriculture %>% mutate(area=rep(c(1552000/2,997000/2,76000/2),each=2)),
state03 %>% mutate(area=c(721000,507000,112000))) other
Activity probabilities for new land uses
<- actprob
actprob3 $landuse <- "forest"
actprob3# Define probabilities for new activities
$defor_to_other <- 0.0002
actprob3$defor_to_agri <- 0.00025
actprob3$affor <- 0
actprob3$donothing <- 0
actprob3# Normalise probabilities after adding new activities
<- c("noman", "thin", "ff", "defor_to_other", "defor_to_agri", "affor", "donothing")
actnames <- actprob3[actnames]/rowSums(actprob3[actnames])
actprob3[actnames] # Define activity probabilities also for agriculture and other land uses
<- rbind(actprob3,
actprob3 %>% mutate(thin=0, ff=0, noman=0, defor_to_other=0, defor_to_agri=0, affor=0.0005, donothing=0.9995),
agriculture %>% mutate(thin=0, ff=0, noman=0, defor_to_other=0, defor_to_agri=0, affor=0, donothing=1)) other
<- runEFDM(state03, actprob3, activities3, 20) states3
#pie charts of land use areas
<- states3 %>%
LUareas group_by(time, landuse) %>%
summarise(area = sum(area)) %>%
ungroup()
#> `summarise()` has grouped output by 'time'. You can override using the `.groups` argument.
<- LUareas %>%
LUareas mutate(time = factor(time, labels = seq(2016,2116,5)))
<- sum(state03$area)
totalarea %>% filter(time %in% c("2016","2066","2116")) %>%
LUareas mutate(area = area/totalarea) %>%
ggplot() +
geom_bar(aes(x="",y=area, fill=landuse), stat="identity", width=1, color = "white") +
scale_fill_viridis_d(end = 0.9) +
coord_polar("y", start=0) +
facet_wrap(vars(time)) +
labs(y=NULL,fill="Land use", x=NULL)+
theme(axis.text.x = element_blank(),
axis.ticks = element_blank()) +
geom_text(aes(x="",y = rep(c(0.03, 0.33, 0.95), 3), label = round(area*100,1)),
nudge_x=0.25, size=2.25,
color=rep(c("grey20","white","white"), 3), fontface="bold") +
theme(legend.position = "bottom")