inst/doc/example.R

In efdm: Simulate Forest Resources with the European Forestry Dynamics Model

## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
options(max.print=50)
data.table::setDTthreads(1)

## -----------------------------------------------------------------------------
library(efdm)

## ----setup, warning = FALSE, message = FALSE----------------------------------
library(dplyr)
library(ggplot2)
library(sf)

## -----------------------------------------------------------------------------
data(example)

## -----------------------------------------------------------------------------
actprob <- example$actprob
actprob

## -----------------------------------------------------------------------------
transprobs_ff <- expand.grid(vol0=1:15, age0=1:35, vol1=1, age1=1, prob=1)
transprobs_ff

## ----ff-----------------------------------------------------------------------
ff <- define_activity("ff", c("vol", "age"), transprobs_ff)

## ----pairdata-----------------------------------------------------------------
example$noman_pairs

## ----statespace---------------------------------------------------------------
statespace <- actprob %>% select(c(region, soil, sp, vol, age))
statespace

## -----------------------------------------------------------------------------
transprobs_noman <- estimatetransprobs(c("vol", "age"), # the dynamic variables of the activity
                     example$noman_pairs,  # pair data for no management
                     statespace=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"))

## ----noman--------------------------------------------------------------------
noman <- define_activity("noman", c("vol", "age"), transprobs_noman)

## -----------------------------------------------------------------------------
transprobs_noman %>%
  arrange(vol0, age0, soil, region, sp, vol1, age1)

## ----thin---------------------------------------------------------------------
transprobs_thin <- estimatetransprobs(c("vol", "age"),
                     example$thin_pairs, # pair data for thinning
                     statespace=statespace,
                     factors=c("soil"),
                     by=c("region","sp"),
                     prior=prior_grow("age"))
thin <- define_activity("thin", c("vol", "age"), transprobs_thin)

## -----------------------------------------------------------------------------
state0 <- example$initial_state
state0

## -----------------------------------------------------------------------------
activities <- list(noman, thin, ff)
states1 <- runEFDM(state0, actprob, activities, 20)

## -----------------------------------------------------------------------------
states1 %>%
  arrange(soil, region, sp, vol, age, time, activity)

## -----------------------------------------------------------------------------
example$vol_coef %>%
  arrange(vol, sp, species)

## -----------------------------------------------------------------------------
example$drain_coef %>%
  arrange(-vol, sp, assort, activity)

## -----------------------------------------------------------------------------
example$income_coef

## -----------------------------------------------------------------------------
states1 <- states1 %>% mutate(time = factor(2016 + time*5))

## -----------------------------------------------------------------------------
volume <- merge(states1, example$vol_coef) %>%
  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))

## -----------------------------------------------------------------------------
states1$ageclass <- cut(states1$age, breaks=c(0,10,20,30,35), include.lowest = TRUE, 
                        #labels=c("0-50","51-100","101-150","150+"))
                        labels=c("-50","-100","-150","150+"))
states1$region <- factor(states1$region, labels = c("South", "Middle", "North"))
ggplot(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))

## -----------------------------------------------------------------------------
euro <- merge(example$drain_coef, example$income_coef) %>% mutate(euro = euro*drain)
removal <- merge(states1, euro) %>% mutate(income=euro*area)

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))

## -----------------------------------------------------------------------------
transprobs_ff_age_species <- statespace %>%
  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(
                         sp0=='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
                                        # in northern boreal vegetation zone 20% of final felled area
                                        # change to spruce dominated forest
                         sp0=='spruce' && region=='South' ~ c(0.3, 0.7),
                         sp0=='spruce' && region=='Middle' ~ c(0.2, 0.8),
                         sp0=='spruce' && region=='North' ~ c(0, 1))))
ff_age_species <- define_activity("ff", c("vol", "age", "sp"), transprobs_ff_age_species)

## -----------------------------------------------------------------------------
activities2 <- list(noman, thin, ff_age_species)
states2 <- runEFDM(state0, actprob, activities2, 20)

## -----------------------------------------------------------------------------
# Compute proportion of spruces in each time and region
prop_spruces <- states2 %>%
  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)

# Add bio-geographical regions (MetsaKasvVyoh provided in efdm package) to the data
prop_spruces <- merge(MetsaKasvVyoh, prop_spruces)
prop_spruces %>% ggplot() +
  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))

## -----------------------------------------------------------------------------
statespace3 <- statespace
statespace3$landuse <- "forest"
agriculture <- unique(statespace3 %>% mutate(sp=0, vol=0, age=0, landuse="agriculture"))
other <- unique(statespace3 %>% mutate(soil=0, sp=0, vol=0, age=0, landuse="other"))
statespace3 <- rbind(statespace3, agriculture, other)

## -----------------------------------------------------------------------------
transprobs_defor_to_agri <- statespace3 %>%
  filter(landuse=="forest") %>%
  select(vol0=vol, age0=age, sp0=sp,landuse0=landuse) %>%
  unique %>%
  mutate(vol1=0, age1=0, sp1=0, landuse1="agriculture", prob=1)
defor_to_agri <- define_activity("defor_to_agri",
                                 c("vol", "age", "sp", "landuse"),
                                 transprobs_defor_to_agri)

## -----------------------------------------------------------------------------

transprobs_defor_to_other <- statespace3 %>%
  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)
defor_to_other <- define_activity("defor_to_other",
                                  c("soil", "vol", "age", "sp", "landuse"),
                                  transprobs_defor_to_other)

## -----------------------------------------------------------------------------
transprobs_aff <- statespace3 %>%
  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)
transprobs_aff <- rbind(transprobs_aff %>% mutate(sp1="spruce", prob=0.5),
                        transprobs_aff %>% mutate(sp1="other", prob=0.5))
affor <- define_activity("affor",
                         c("vol", "age", "sp", "landuse"),
                         transprobs_aff)

## -----------------------------------------------------------------------------
donothing <- define_activity("donothing", character())
activities3 <- list(noman, thin, ff, defor_to_other, defor_to_agri, affor, donothing)

## -----------------------------------------------------------------------------
state03 <- state0
state03$landuse <- "forest"
state03 <- rbind(state03,  agriculture %>% mutate(area=rep(c(1552000/2,997000/2,76000/2),each=2)),
                 other %>% mutate(area=c(721000,507000,112000)))

## -----------------------------------------------------------------------------
actprob3 <- actprob
actprob3$landuse <- "forest"
# Define probabilities for new activities
actprob3$defor_to_other <- 0.0002
actprob3$defor_to_agri <- 0.00025
actprob3$affor <- 0
actprob3$donothing <- 0
# Normalise probabilities after adding new activities
actnames <- c("noman", "thin", "ff", "defor_to_other", "defor_to_agri", "affor", "donothing")
actprob3[actnames] <- actprob3[actnames]/rowSums(actprob3[actnames])
# Define activity probabilities also for agriculture and other land uses
actprob3 <- rbind(actprob3, 
        agriculture %>% mutate(thin=0, ff=0, noman=0, defor_to_other=0, defor_to_agri=0, affor=0.0005, donothing=0.9995),
        other %>% mutate(thin=0, ff=0, noman=0, defor_to_other=0, defor_to_agri=0, affor=0, donothing=1))

## -----------------------------------------------------------------------------
states3 <- runEFDM(state03, actprob3, activities3, 20)

## -----------------------------------------------------------------------------
#pie charts of land use areas
LUareas <- states3 %>%
  group_by(time, landuse) %>%
  summarise(area = sum(area)) %>%
  ungroup()
LUareas <- LUareas %>%
  mutate(time = factor(time, labels = seq(2016,2116,5)))
totalarea <- sum(state03$area)
LUareas %>% filter(time %in% c("2016","2066","2116")) %>%
  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")