# sitree: An Individual Tree Open Source Simulator

#### Clara Antón Fernández

#### 03-08-2021

# Introduction

This tutorial gives examples on how you can use sitree to simulate individual tree growth, death, and harvest through time, and to compare the effects of changing functions (e.g. use a different mortality function) in the results of the simulation.

For a general introduction and description of the package take a look at the ‘Introduction to SiTree’ vignette.

# A little note on Reference classes

SiTree uses reference classes to keep memory usage at bay. Bare in mind that reference classes such as trList and trListDead are modified in place. For example,

    ## Always bare in mind that RF trList and trListDead are modified in place
    ## Do not copy

    Polygon <- setRefClass("Polygon", fields = c("sides"))
    square <- Polygon$new(sides = 4)

    square.to.triangle <- function(a){
      ## This makes the object triangle point to square, it does not create a copy
      ## since square is a RC object
      triangle <- square
      ## when we modifiy triangle we are modifying square.
      triangle$sides <- 3
      return(a+2)
    }

    ## square has 4 sides
    square$sides

    ## [1] 4

    ## when we call the square.to.triangle function the 'square' object
    ## gets modified, eventhough we don't pass it as argument. That is because
    ## it is referenced inside the function
    square.to.triangle(2)

    ## [1] 4

    square$sides

    ## [1] 3

    ## but if we do
    square <- Polygon$new(sides = 4)

    square.to.triangle <- function(a){
        triangle <- square$copy()
        triangle$sides <- 3
        return(a+2)
    }

    square$sides

    ## [1] 4

    square.to.triangle(2)

    ## [1] 4

    ## the object remains unchanged
    square$sides

    ## [1] 4

# Methods available for trList and trListDead (advanced)

Most regular users will never need to use these methods.

Let's see now some examples of the methods implemented for these RF.

    library(sitree)
    res <- sitree (tree.df   = tr,
                     stand.df  = fl,
                     functions = list(
                         fn.growth     = 'grow.dbhinc.hgtinc',
                         fn.mort       = 'mort.B2007',
                         fn.recr       = 'recr.BBG2008',
                         fn.management = 'management.prob',
                         fn.tree.removal = 'mng.tree.removal',
                         fn.modif      = NULL,
                         fn.prep.common.vars = 'prep.common.vars.fun'
                     ),
                     n.periods = 5,
                     period.length = 5,
                     mng.options = NA,
                     print.comments = FALSE,
                     fn.dbh.inc = "dbhi.BN2009",
                     fn.hgt.inc =  "height.korf",
                     fun.final.felling = "harv.prob",
                     fun.thinning      = "thin.prob",
                     per.vol.harv = 0.83
                     )

     ## getTrees(i, j)  -- obtains the information of the i trees, on the j periods,
     ## by default it selects all. It does not display it, it passes the value.
     ## It returns a list with elements plot.id, treeid, dbh.mm, height.dm, yrs.sim,
     ## tree.sp

     get.some.trees <- res$live$getTrees(1:3, 2:5)

     ## extractTrees(i)  -- extracts the i trees, it removed the trees from the
     ## original object and it passes the information. It returns a list.

     dead <- res$live$extractTrees(4:7)

     ## addTrees(x) -- x should be a list

     res$live$addTrees(dead)

     ## last.time.alive. It checks when was the last DBH measured.
     new.dead.trees <- trListDead$new(
          data = dead,
          last.measurement = cbind(
            do.call("dead.trees.growth"
                  , args = list(
                      dt     = dead,
                      growth = data.frame(dbh.inc.mm     = rep(3, 4),
                                          hgt.inc.dm  = rep(8, 4)),
                      mort   = rep(TRUE, 4),
                      this.period = "t2")
                    ),
            found.dead = "t3"
          ),
          nperiods = res$live$nperiods
          )

     lta <- new.dead.trees$last.time.alive()

     ## which in this case differs from the data stored under the last.measurement
     ## field because we have defined it artificially above as "t3"
     lta

    ##    5    6    8    9
    ## "t5" "t5" "t5" "t5"

     new.dead.trees$last.measurement$found.dead

    ## [1] "t3" "t3" "t3" "t3"

     ## But we can remove the data from the periods after it was found dead
     new.dead.trees$remove.next.period("t3")
     new.dead.trees$remove.next.period("t4")
     new.dead.trees$remove.next.period("t5")
     ## and now results do match
     lta <- new.dead.trees$last.time.alive()
     ## last time it was alive was in "t2"
     lta

    ##    5    6    8    9
    ## "t2" "t2" "t2" "t2"

     ## ant it was found dead in "t3"
     new.dead.trees$last.measurement$found.dead

    ## [1] "t3" "t3" "t3" "t3"

# Taking a look at the results of the simulation

We first use the example functions available in the sitree package to run a 15 periods simulation.

    set.seed(2017)

    n.periods <- 50
    res <- sitree (tree.df   = tr,
                    stand.df  = fl,
                    functions = list(
                        fn.growth     = 'grow.dbhinc.hgtinc',
                        fn.mort       = 'mort.B2007',
                        fn.recr       = 'recr.BBG2008',
                        fn.management = 'management.prob',
                        fn.tree.removal = 'mng.tree.removal',
                        fn.modif      = NULL,
                        fn.prep.common.vars = 'prep.common.vars.fun'
                    ),
                    n.periods = n.periods,
                    period.length = 5,
                    mng.options = NA,
                    print.comments = FALSE,
                    fn.dbh.inc = "dbhi.BN2009",
                    fn.hgt.inc =  "height.korf",
                    fun.final.felling = "harv.prob",
                    fun.thinning      = "thin.prob",
                    per.vol.harv = 0.83
                    )

Let’s start with some basic graphics of the data.

## dbh distribution through time

    dbh.mm <- res$live$data$dbh.mm
    dbh.mm.short <- reshape(dbh.mm,
                            varying = paste0("t", 0:n.periods),
                            timevar = "period",
                            direction = "long", sep = "")
    head(dbh.mm.short)

    ##     period   t id
    ## 1.0      0 234  1
    ## 2.0      0 173  2
    ## 3.0      0  86  3
    ## 4.0      0 230  4
    ## 5.0      0 236  5
    ## 6.0      0 234  6

    dbh.mm.short$t[dbh.mm.short$t == 0] <- NA
    library(ggplot2)
    ggplot(dbh.mm.short, aes(x = t)) + geom_histogram() + ylab('dbh.mm') +
       facet_wrap(~ period) + theme_minimal()

    ## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

    ## Warning: Removed 255723 rows containing non-finite values (stat_bin).

## Volume distribution through time

And since we already have the code to calculate volume we can easily plot the evolution of standing volume or harvested volume

    vol <- data.frame(matrix(NA, ncol = n.periods+1, nrow = length(res$plot.data$plot.id)))
    names(vol) <- paste0("t", 0:n.periods)
    for (i.period in 0:n.periods){
        sa <- prep.common.vars.fun (
            tr = res$live, fl= res$plot.data,
            i.period, this.period = paste0("t", i.period),
            common.vars = NULL,
            period.length = 5 )

        vol[, (i.period +1)] <- sa$res$vuprha.m3.ha
        ## This is volume per ha, if we prefer just m3 we multiply by ha2total
        ## ha2total is the number of ha represented by the plot
        vol[, (i.period +1)] <- vol[, (i.period +1)] * sa$fl$ha2total
    }
    vol.m3.short <- reshape(vol,
                            varying = paste0("t", 0:n.periods),
                            timevar = "period",
                            direction = "long",
                        sep = "")
    vol.m3.short$t[vol.m3.short$t == 0] <- NA
    total.standing.volume <- aggregate(t ~ period, data = vol.m3.short, FUN = sum)
    total.standing.volume$vol.mill.m3  <- total.standing.volume$t/1e6

    ggplot(total.standing.volume,
           aes(period, vol.mill.m3))   + geom_line()+
      ylab("standing volume (mill m3)") + theme_minimal()

## Harvested volume

    vol <- data.frame(matrix(NA, ncol = n.periods + 1,
                             nrow = length(res$plot.data$plot.id)))
    ## res$removed$data only contains the "history of the tree", but we need
    ## the dbh and height of the tree at harvest time
    names(vol) <- paste0("t", 0:n.periods)
    removed <- recover.last.measurement(res$removed)

    for (i.period in 0:n.periods){
        harv.vol <- prep.common.vars.fun (
            tr = res$removed,
            fl = res$plot.data,
            i.period,
            this.period = paste0("t", i.period),
            common.vars = NULL,
            period.length = 5 )

        vol[, (i.period +1)] <- harv.vol$res$vuprha.m3.ha
        ## This is volume per ha, if we prefer just m3.
        vol[, (i.period +1)] <- vol[, (i.period +1)] * sa$fl$ha2total
    }
    names(vol) <- paste0(substr(names(vol), 1, 1), ".", substr(names(vol), 2, 3))

    vol.m3.short <- reshape(vol,
                            varying = paste0("t.", 0:n.periods),
                            timevar = "period",
                            idvar = "id",
                            direction = "long"
                            )
    vol.m3.short$t[vol.m3.short$t == 0] <- NA
    harv.total <- aggregate(t ~ period, data = vol.m3.short, FUN = sum)

    ggplot(harv.total, aes(period, t))   + geom_line()+
      ylab("harvested volume ( m3)") + theme_minimal()

## Volume of dead trees

    vol <- data.frame(matrix(NA, ncol = n.periods + 1,
                             nrow = length(res$plot.data$plot.id)))
    names(vol) <- paste0("t", 0:n.periods)
    dead <- recover.last.measurement(res$dead)

    for (i.period in 0:n.periods){
      vol[, (i.period +1)] <-
        prep.common.vars.fun (
          tr = dead,
          fl= res$plot.data,
          i.period,
          this.period = paste0("t", i.period),
          common.vars = NULL,
          period.length = 5 )$res$vuprha.m3.ha

        ## This is volume per ha, if we prefer just m3.
        vol[, (i.period +1)] <- vol[, (i.period +1)] * sa$fl$ha2total
    }

    names(vol) <- paste0(substr(names(vol), 1, 1), ".", substr(names(vol), 2, 3))
    vol.m3.short <- reshape(vol,
                            varying = paste0("t.", 0:n.periods),
                            timevar = "period",
                            idvar = "id",
                            direction = "long"
                            )

    head(vol.m3.short)

    ##     period t id
    ## 1.0      0 0  1
    ## 2.0      0 0  2
    ## 3.0      0 0  3
    ## 4.0      0 0  4
    ## 5.0      0 0  5
    ## 6.0      0 0  6

    ## let's plot dead trees by plot, for the first 10 plots
    ## to look at the variation
    ggplot(vol.m3.short[vol.m3.short$id %in% 1:10,],
           aes(period, t, group = id, col = as.factor(id)))   +
      geom_line()+
      ylab("Dead trees volume (mill m3)") + theme_minimal()

## Stand age

It is also interesting to look at the evolution of stand age

    age <- res$plot.data$stand.age.years
    age.short<- reshape(age,
                        varying = paste0("t", 0:(n.periods-1)),
                        timevar = "period",
                        idvar = "id",
                        direction = "long",
                        sep = ""
                        )
    head(age.short)

    ##     period  t id
    ## 1.0      0 33  1
    ## 2.0      0 62  2
    ## 3.0      0 67  3
    ## 4.0      0 69  4
    ## 5.0      0  5  5
    ## 6.0      0 67  6

    ggplot(age.short, aes(x = t)) + geom_histogram() + ylab('dbh.mm') +
       facet_wrap(~ period) + theme_minimal()

    ## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# Switching functions

But the real power of sitree is how easily one can test new functions and its effects on important forestry variables. Let's see what happens when we switch mortality functions.

So far we have been using ‘mort.B2007’, Bollandås (2007) mortality equation. we will now try older Norwegian mortality functions. Those developed by Eid and Tuhus (2001).

    ET2001 <- function (tr, fl, common.vars, this.period, ...) {
        dbh.mm <- tr$data[["dbh.mm"]][, this.period]
        p.functions <-
            data.frame(a0 = c( 8.06,   8.49, 4.89,  5.16),
                       b1 = c(-6.7, -14.27, -2.528, -7.35),
                       b2 = c(-0.03,  -0.05, 0,     -0.02),
                       b3 = c(-0.03, -0.08, 0,     0 ),
                       b4 = c(-0.01, -0.08, 0,     0 )
                       )
        rownames(p.functions) <- c("spruce", "pine", "birch", "other")

        logit <- p.functions[common.vars$spp, "a0"] +
            p.functions[common.vars$spp, "b1"] * (dbh.mm/10)^(-1) +
            p.functions[common.vars$spp, "b2"] * common.vars$PBAL.m2.ha +
            p.functions[common.vars$spp, "b3"] * fl$SI.m[common.vars$i.stand] +
            p.functions[common.vars$spp, "b4"] * common.vars$pr.spp.ba$spru
        mort.B <- 1- (1/(1 + exp(-logit)))
        mort <- ifelse(mort.B >= runif(nrow(tr$data$dbh.mm),  0, 1), TRUE, FALSE)
        return(mort)
    }

We will test the performance of this new function by simply comparing the development of the stand basal area for 80 years. Since both `ET2001` and `mort.B2007` and implemented stochastically, we will average the stand basal area of 20 runs. This could also be implemented easily in parallel. We first run the simulation with the `ET2001` function:

    SBA.m2.ha <- list()
    set.seed(2017)
    for ( i.to in 1:20){

      result.sitree.ET2001 <-
        sitree (tree.df   = tr,
                stand.df  = fl,
                functions = list(
                  fn.growth     = 'grow.dbhinc.hgtinc',
                  fn.mort       = 'ET2001',
                  fn.recr       = 'recr.BBG2008',
                  fn.management = NULL,
                  fn.tree.removal = NULL,
                  fn.modif      = NULL,
                  fn.prep.common.vars = 'prep.common.vars.fun'
                ),
                n.periods = 16,
                period.length = 5,
                mng.options = NA,
                print.comments = FALSE,
                fn.dbh.inc    = 'dbhi.BN2009',
                fn.hgt.inc    = 'height.korf'
                )

      ## We extract basal area from sitree.summary
      SBA.m2.ha[[i.to]] <-
        attr(sitree.summary(
          result.sitree.ET2001,
          plots = c(1),
          by.stand = TRUE,
          plot.all.together = TRUE,
          plot = FALSE), "data")$SBA.m2.ha
    }

    library(data.table)
    SBA.m2.ha.ET2001 <- rbindlist(SBA.m2.ha)
    SBA.m2.ET2001 <- aggregate(t ~ period + plot.id, data = SBA.m2.ha.ET2001,
                                  FUN = mean)

We run the simulation now with mortality function `mort.B2007` implemented in sitree:

    SBA.m2.ha <- list()
    set.seed(2017)
    for ( i.to in 1:20){

      result.sitree <-
        sitree (tree.df   = tr,
                stand.df  = fl,
                functions = list(
                  fn.growth     = 'grow.dbhinc.hgtinc',
                  fn.mort       = 'mort.B2007',
                  fn.recr       = 'recr.BBG2008',
                  fn.management = NULL,
                  fn.tree.removal = NULL,
                  fn.modif      = NULL,
                  fn.prep.common.vars = 'prep.common.vars.fun'
                ),
                n.periods = 16,
                period.length = 5,
                mng.options = NA,
                print.comments = FALSE,
                fn.dbh.inc    = 'dbhi.BN2009',
                fn.hgt.inc    = 'height.korf'
                )
      ## Basal area
      SBA.m2.ha[[i.to]] <- attr(sitree.summary(
        result.sitree,
        plots = c(1),
        by.stand = TRUE,
        plot.all.together = TRUE,
        plot = FALSE), "data")$SBA.m2.ha

    }

    SBA.m2.ha.B2007 <- rbindlist(SBA.m2.ha)
    SBA.m2.B2007 <- aggregate(t ~ period + plot.id, data = SBA.m2.ha.B2007,
                                  FUN = mean)

    SBA.m2.B2007$mort <- "mort.B2007"
    SBA.m2.ET2001$mort <- "ET2001"

We can now visualize the stand basal area development for the first 10 plots using the different mortality functions:

    mort.data <- rbind(SBA.m2.B2007, SBA.m2.ET2001)
    mort.data$uid  <- as.factor(paste0(mort.data$plot.id, mort.data$mort))
    ggplot(mort.data[mort.data$plot.id %in% 1:10,],
           aes(  period, t, col = mort, group = uid)) +
      geom_line() + ylab ( "Stand basal area (m2)")+
      theme(legend.title = element_blank())

The resulting figure shows the differences in stand basal area between the two mortality sub-models for the first 10 plots.

# References

Bollandsås, O. 2007 Uneven-aged Forestry in Norway: Inventory and Management Models. [Ås, Norway]: Norwegian University of Life Sciences, Department of Ecology and Natural Resource Management.

Eid, T., and Tuhus E. 2001 Models for Individual Tree Mortality in Norway. Forest Ecology and Management 154 (1/2): 69–84.