# inst/doc/TestEquations.R In sitree: Single Tree Simulator

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

la <- function(a){
triangle <- square
triangle\$sides <- 3
return(a+2)
}

square\$sides
la(2)

square\$sides

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

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

square\$sides
la(2)
## the object remains unchanged
square\$sides

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

require(sitree)
res <- sitree (tree.df   = tr,
stand.df  = fl,
functions = list(
fn.growth     = 'grow.dbhinc.hgtinc',
fn.dbh.inc = "dbhi.BN2009",
fn.hgt.inc =  "height.korf",
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,
species.spruce = c(1, 2, 3),
species.pine = c(10, 11, 20, 21, 29),
species.harw = c(30, 31),
fun.final.felling = "harv.prob",
fun.thinning      = "thin.prob",
'BN2009',
'BBG2008', 'SBA.m2.ha', 'spp','pr.spru.ba', 'QMD.cm',
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.

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

## last.time.alive. It checks when was the last DBH measured.
last.measurement = cbind(
, args = list(
growth = data.frame(dbh.inc.mm     = rep(3, 4),
hgt.inc.dm  = rep(8, 4)),
mort   = rep(TRUE, 4),
this.period = "t2")
),
),
nperiods = res\$live\$nperiods
)

## which in this case differs from the data stored under the last.measurement
## field because we have defined it artificially above as "t3"
lta
## But we can remove the data from the periods after it was found dead
## and now results do match
## last time it was alive was in "t2"
lta
## ant it was found dead in "t3"

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

set.seed(2017)

res <- sitree (tree.df   = tr,
stand.df  = fl,
functions = list(
fn.growth     = 'grow.dbhinc.hgtinc',
fn.dbh.inc = "dbhi.BN2009",
fn.hgt.inc =  "height.korf",
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,
species.spruce = c(1, 2, 3),
species.pine = c(10, 11, 20, 21, 29),
species.harw = c(30, 31),
fun.final.felling = "harv.prob",
fun.thinning      = "thin.prob",
'BN2009',
'BBG2008', 'SBA.m2.ha', 'spp','pr.spru.ba', 'QMD.cm',
per.vol.harv = 0.83
)

## ------------------------------------------------------------------------
require(lattice)

dbh.mm <- res\$live\$data\$dbh.mm
dbh.mm.short <- reshape(dbh.mm,
varying = paste0("t", 0:5),
timevar = "period",
direction = "long", sep = "")
dbh.mm.short\$t[dbh.mm.short\$t == 0] <- NA
histogram( ~ t | period, data = dbh.mm.short, plot.points = FALSE,
ref = TRUE, auto.key = list(space = "right"), xmin = 50,
xlab = "dbh.mm")

## ------------------------------------------------------------------------
vol <- data.frame(matrix(NA, ncol = 6, nrow = length(res\$plot.data\$plot.id)))
names(vol) <- paste0("t", 0:5)
for (i.period in 0:5){
sa <- prep.common.vars.fun (
tr = res\$live, fl= res\$plot.data,
i.period, this.period = paste0("t", i.period),
common.vars = NULL, vars.required = "vuprha.m3.ha",
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:5),
timevar = "period",
direction = "long",
sep = "")
vol.m3.short\$t[vol.m3.short\$t == 0] <- NA
harv.total <- aggregate(t ~ period, data = vol.m3.short, FUN = sum)
xyplot( t/1e6 ~ period, data = harv.total, type = 'l',
ylab = "standing volume in mill. m3")

## ------------------------------------------------------------------------
vol <- data.frame(matrix(NA, ncol = 6, 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:5)
removed <- recover.last.measurement(res\$removed)

for (i.period in 0:5){
sa <- prep.common.vars.fun (
tr = res\$removed,
fl = res\$plot.data,
i.period,
this.period = paste0("t", i.period),
common.vars = NULL,
vars.required = "vuprha.m3.ha",
period.length = 5 )

vol[, (i.period +1)] <- sa\$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:5),
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)
xyplot( t/1e6 ~ period, data = harv.total, type = 'l',
ylim = c(0, max(harv.total\$t/1e6)+0.1),
ylab = "harvested volume in mill. m3")

## ------------------------------------------------------------------------
vol <- data.frame(matrix(NA, ncol = 6, nrow = length(res\$plot.data\$plot.id)))
names(vol) <- paste0("t", 0:5)

for (i.period in 0:5){
sa <- prep.common.vars.fun (
i.period, this.period = paste0("t", i.period),
common.vars = NULL, vars.required = "vuprha.m3.ha",
period.length = 5 )

vol[, (i.period +1)] <- sa\$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:5),
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)
xyplot( t/1e6 ~ period, data = harv.total, type = 'l',
ylim = c(0, max(harv.total\$t/1e6)+10),
ylab = "volume of dead trees in mill. m3")

## ------------------------------------------------------------------------
age <- res\$plot.data\$stand.age.years
age.short<- reshape(age,
varying = paste0("t", 0:40),
timevar = "period",
idvar = "id",
direction = "long",
sep = ""
)
age.short <- age.short[!age.short\$period > 5,]
histogram( ~ t | period, data = age.short, plot.points = FALSE,
ref = TRUE, auto.key = list(space = "right"), xmin = 50,
xlab = "stand age (in years)")

## ------------------------------------------------------------------------
ET2001 <- function (tr, fl, common.vars, this.period, ...)
{
if (!all(unique(common.vars\$spp) %in% c("spruce", "pine",
"birch", "other"))) {
message("spp should only contain values spruce, pine, birch, other")
}
dbh.mm <- tr\$data[["dbh.mm"]][, this.period]
p.functions <-
data.frame(a0 = c( 8.059,   8.4904, 4.892,  5.157),
b1 = c(-6.702, -14.266, -2.528, -7.3544),
b2 = c(-0.028,  -0.0462, 0,     -0.0199),
b3 = c(-0.0264, -0.0761, 0,     0 ),
b4 = c(-0.0132, -0.0761, 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(length(tr\$data[["plot.id"]]),
0, 1), TRUE, FALSE)
return(mort)
}

## ----eval = FALSE--------------------------------------------------------
#  set.seed(2017)
#  i.max <- 50
#  vol.1 <- data.frame(matrix(NA, ncol = 6, nrow = i.max))
#  for (i in (1:i.max)){
#      print(i)
#      res1 <- sitree (tree.df   = tr,
#                       stand.df  = fl,
#                       functions = list(
#                           fn.growth     = 'grow.dbhinc.hgtinc',
#                           fn.dbh.inc = "dbhi.BN2009",
#                           fn.hgt.inc =  "height.korf",
#                           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,
#                       species.spruce = c(1, 2, 3),
#                       species.pine = c(10, 11, 20, 21, 29),
#                       species.harw = c(30, 31),
#                       fun.final.felling = "harv.prob",
#                       fun.thinning      = "thin.prob",
#                       'BN2009',
#                       'BBG2008', 'SBA.m2.ha', 'spp','pr.spru.ba', 'QMD.cm',
#                       per.vol.harv = 0.83
#                       )
#
#      names(vol) <- paste0("t", 0:5)
#      for (i.period in 0:5){
#          sa <- prep.common.vars.fun (
#              tr = res1\$live, fl= res1\$plot.data,
#              i.period, this.period = paste0("t", i.period),
#              common.vars = NULL, vars.required = "vuprha.m3.ha",
#              period.length = 5 )
#
#          vol.1[i, (i.period +1)] <- sum(sa\$res\$vuprha.m3.ha * sa\$fl\$ha2total )
#
#      }
#  }

## ----eval = FALSE--------------------------------------------------------
#
#  vol.2 <- data.frame(matrix(NA, ncol = 6, nrow = i.max))
#  for (i in (1:i.max)){
#      print(i)
#
#      res2 <- sitree (tree.df   = tr,
#                       stand.df  = fl,
#                       functions = list(
#                           fn.growth     = 'grow.dbhinc.hgtinc',
#                           fn.dbh.inc = "dbhi.BN2009",
#                           fn.hgt.inc =  "height.korf",
#                           fn.mort       = 'ET2001',
#                           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,
#                       species.spruce = c(1, 2, 3),
#                       species.pine = c(10, 11, 20, 21, 29),
#                       species.harw = c(30, 31),
#                       fun.final.felling = "harv.prob",
#                       fun.thinning      = "thin.prob",
#                       'BN2009','PBAL.m2.ha',
#                       'BBG2008', 'SBA.m2.ha', 'spp','pr.spru.ba', 'QMD.cm',
#                       per.vol.harv = 0.83
#                       )
#      names(vol) <- paste0("t", 0:5)
#      for (i.period in 0:5){
#          sa <- prep.common.vars.fun (
#              tr = res1\$live, fl= res1\$plot.data,
#              i.period, this.period = paste0("t", i.period),
#              common.vars = NULL, vars.required = "vuprha.m3.ha",
#              period.length = 5 )
#
#          vol.2[i, (i.period +1)] <- sum(sa\$res\$vuprha.m3.ha * sa\$fl\$ha2total )
#
#      }
#  }
#  ## During the five periods, there is actually not that much of a difference
#  colMeans(vol.1)/1e6
#  colMeans(vol.2)/1e6
#
#
```

## Try the sitree package in your browser

Any scripts or data that you put into this service are public.

sitree documentation built on June 1, 2018, 5:06 p.m.