## -----------------------------------------------------------------------------
## 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
## 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)
square$sides
## 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
square.to.triangle(2)
## the object remains unchanged
square$sides
## -----------------------------------------------------------------------------
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
new.dead.trees$last.measurement$found.dead
## 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
## ant it was found dead in "t3"
new.dead.trees$last.measurement$found.dead
## -----------------------------------------------------------------------------
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
)
## -----------------------------------------------------------------------------
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)
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()
## -----------------------------------------------------------------------------
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()
## -----------------------------------------------------------------------------
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()
## -----------------------------------------------------------------------------
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)
## 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()
## -----------------------------------------------------------------------------
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)
ggplot(age.short, aes(x = t)) + geom_histogram() + ylab('dbh.mm') +
facet_wrap(~ period) + theme_minimal()
## -----------------------------------------------------------------------------
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)
}
## -----------------------------------------------------------------------------
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)
## -----------------------------------------------------------------------------
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"
## -----------------------------------------------------------------------------
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())