sitree.summary <- function(sitrees.res, plots, by.stand = TRUE, plot = FALSE,
plot.all.together = FALSE, with.legend = FALSE){
## check if plot.size.m2
if(! 'plot.size.m2' %in% names(sitrees.res$plot.data))
stop('plot.data in sitrees.res$plot.data should have a field named plot.size.m2')
plot.data <- data.frame(plot.id = sitrees.res$plot.data$plot.id,
plot.size.m2 = sitrees.res$plot.data$plot.size.m2)
plot.data$tree2ha <- 10000/plot.data$plot.size.m2
nplots <- nrow(plot.data)
res <- list()
############################
## PREPARING THE DATA
##############################
## plot1 -- BA per ha
res$i.stand <- match(sitrees.res$live$data$plot.id,
plot.data$plot.id)
## NA mean no trees in that plot
res$i.tree <- match(plot.data$plot.id,
sitrees.res$liv$data$plot.id)
res$tree2ha <- 10000 / plot.data$plot.size.m2
res$SBA.m2.ha <- res$stems.ha <- res$heights.10 <-
data.frame(t0 = matrix(NA, nrow = nplots))
## initialize my.plots
my.plots <- list(plot1 = NULL, plot2 = NULL, plot3 = NULL, plot4 = NULL,
plot5 = NULL)
data.summary <- list()
for (i in 0:sitrees.res$live$nperiods){
this.period <- paste0("t", i)
res$tree.BA.m2 <- pi * (sitrees.res$live$data$dbh.mm[, this.period]/1000/2)^2
## SBA.m2.ha
sa <-
aggregate(res$tree.BA.m2 * res$tree2ha[res$i.stand],
by = list(plot.id = sitrees.res$live$data$plot.id),
FUN = sum)
res$SBA.m2.ha[, this.period] <- sa$x[ match(plot.data$plot.id,
sa$plot.id)]
## stems.ha
living.trees <- (sitrees.res$live$data$dbh.mm[, this.period] > 0)
sa <-
aggregate(res$tree2ha[res$i.stand][living.trees],
by = list(plot.id =
sitrees.res$live$data$plot.id[living.trees]),
FUN = sum)
res$stems.ha[, this.period] <- sa$x[ match(plot.data$plot.id,
sa$plot.id)]
## height of 10 tallest trees
## stems.ha
heights <- data.frame(
heights = sitrees.res$live$data$height.dm[, this.period],
plot.id = sitrees.res$live$data$plot.id
)
heights <- heights[heights$heights != 0,]
heights <- heights[order(heights$plot.id, heights$heights),]
heights <- aggregate(heights ~ plot.id, data = heights,
FUN = function(x) ifelse (length(x) > 9,
mean(x[1:10]),
mean(x[1:length(x)])
)
)
res$heights.10[, this.period] <-
heights$heights [ match(plot.data$plot.id,
heights$plot.id)]
}
## number of dead trees
dead.trees <- data.frame(
period = sitrees.res$dead$last.measurement$found.dead,
plot.id = sitrees.res$dead$data$plot.id)
num.dead.trees.ha <-
aggregate(dead.trees$period,
by = list(plot.id = dead.trees$plot.id,
period = dead.trees$period),
FUN = length)
num.dead.trees.ha$x <- num.dead.trees.ha$x *
res$tree2ha[match(num.dead.trees.ha$plot.id, plot.data$plot.id)]
## number of removed trees
if (!is.null(sitrees.res$removed)){
removed.trees <- data.frame(
period = sitrees.res$removed$last.measurement$found.removed,
plot.id = sitrees.res$removed$data$plot.id)
num.removed.trees.ha <-
aggregate(removed.trees$period,
by = list(plot.id = removed.trees$plot.id,
period = removed.trees$period),
FUN = length)
num.removed.trees.ha$x <- num.removed.trees.ha$x *
res$tree2ha[match(num.removed.trees.ha$plot.id, plot.data$plot.id)]
} else {
num.removed.trees.ha <- expand.grid(plot.id = plot.data$plot.id,
period = colnames(sitrees.res$live$data$dbh.mm))
num.removed.trees.ha$x <- 0
}
my.period.levels <- paste0("t", 0:sitrees.res$live$nperiods)
############################
## PREPARING THE DATA -- DIFFERENTIATE BETWEEN BY PLOT OR GLOBAL
##############################
## plot1 -- SBA.m2.ha
if (1 %in% plots){
## we need to make the data long
SBA.m2.ha <- reshape( res$SBA.m2.ha,
varying = names(res$SBA.m2.ha),
timevar = "period",
direction = "long", sep = "")
SBA.m2.ha$period <- factor(paste0("t", SBA.m2.ha$period),
levels = my.period.levels)
SBA.m2.ha <- data.table(SBA.m2.ha)
SBA.m2.ha$plot.id <- SBA.m2.ha$id
SBA.m2.ha[plot.data, plot.size.m2 := i.plot.size.m2, on = 'plot.id']
SBA.m2.ha$plot.id <- as.factor(SBA.m2.ha$id)
if (by.stand){
plot1 <-
ggplot(SBA.m2.ha,
aes(x = period, y = t, color = plot.id, group=plot.id)) +
geom_line() + ylab("SBA (m2/ha)")+ ylim(0,NA)+ xlab('Simulation period') +
guides(color = guide_legend(title = "Plot ID"))
} else{
SBA.m2.ha <- SBA.m2.ha[, list(t = weighted.mean(t, w = plot.size.m2)),
by = period]
plot1 <-
ggplot(SBA.m2.ha, aes(period, t, group = 1)) +
ylim(0,NA)+ xlab('Simulation period') +
geom_line() + ylab("Average SBA (m2/ha)")
}
if (with.legend == FALSE){
plot1 <- plot1 + theme(legend.position="none")
}
my.plots$plot1 <- plot1
data.summary$SBA.m2.ha <- SBA.m2.ha
}
## plot2 -- stems per ha
## we need to make the data long
if (2 %in% plots){
stems.ha <- reshape( res$stems.ha,
varying = names(res$stems.ha),
timevar = "period",
direction = "long", sep = "")
stems.ha$period <- factor(paste0("t", stems.ha$period),
levels = my.period.levels)
stems.ha <- stems.ha[!is.na(stems.ha$t),]
stems.ha <- data.table(stems.ha)
stems.ha$plot.id <- stems.ha$id
stems.ha[plot.data, plot.size.m2 := i.plot.size.m2, on = 'plot.id']
stems.ha$plot.id <- as.factor(stems.ha$id)
if (by.stand){
my.plots$plot2 <-
ggplot(stems.ha, aes(x = period, y = t,
color = plot.id, group = plot.id)) + xlab('Simulation period') +
geom_line() + ylab("stems/ha")+ ylim(0,NA)+
guides(color = guide_legend(title = "Plot ID"))
} else {
stems.ha <- stems.ha[, list(t = weighted.mean(t, w = plot.size.m2)),
by = period]
my.plots$plot2 <-
ggplot(stems.ha , aes(period, t, group = 1)) +
geom_line() +ylim(0,NA)+ xlab('Simulation period') +
ylab("Average stems/ha")
}
if (with.legend == FALSE){
my.plots$plot2 <- my.plots$plot2 + theme(legend.position="none")
}
data.summary$stems.ha <- stems.ha
}
## plot3 -- height of 10 tallest trees
if (3 %in% plots){
heights.10 <- reshape( res$heights.10,
varying = names(res$heights.10),
timevar = "period",
direction = "long", sep = "")
heights.10$period <- factor(paste0("t", heights.10$period),
levels = my.period.levels)
heights.10 <- data.table(heights.10)
heights.10$plot.id <- heights.10$id
heights.10[plot.data, plot.size.m2 := i.plot.size.m2, on = 'plot.id']
heights.10$plot.id <- as.factor(heights.10$id)
## if there are no trees, we don't plot anything
heights.10 <- heights.10[!is.na(heights.10$t),]
if (by.stand){
my.plots$plot3 <-
ggplot(heights.10, aes(x = period, y = t, color = plot.id, group = plot.id)) +
geom_line() + xlab('Simulation period') +
ylab("Average height of the 10 tallest trees (dm)")+ylim(0,NA)+
guides(color = guide_legend(title = "Plot ID"))
} else {
heights.10 <- heights.10[, list(t = weighted.mean(t, w = plot.size.m2)),
by = period]
my.plots$plot3 <-
ggplot(heights.10, aes(x = period, y = t, group = 1)) +
geom_line() +ylim(0,NA)+ xlab('Simulation period') +
ylab("Average height of the 10 tallest trees per plot (dm)")
}
if (with.legend == FALSE){
my.plots$plot3 <- my.plots$plot3 + theme(legend.position="none")
}
data.summary$heights.10 <- heights.10
}
## plot4 --number of dead
if (4 %in% plots){
dead.trees <- expand.grid(plot.id = plot.data$plot.id,
period = colnames(sitrees.res$live$data$dbh.mm))
dead.trees$num.dead.trees.ha <- 0
dead.trees <- data.table(dead.trees)
dead.trees[num.dead.trees.ha, num.dead.trees.ha := i.x, on = list(plot.id, period)]
dead.trees[plot.data, tree2ha := i.tree2ha, on = 'plot.id']
dead.trees[plot.data, plot.size.m2 := i.plot.size.m2, on = 'plot.id']
dead.trees$plot.id <- as.factor(dead.trees$plot.id)
if (by.stand){
my.plots$plot4 <-
ggplot(dead.trees, aes(x = period, y = num.dead.trees.ha ,
color = plot.id, group = plot.id)) +ylim(0,NA)+
xlab('Simulation period') +
geom_line() + ylab( "dead trees per ha")+
guides(color = guide_legend(title = "Plot ID"))
} else {
num.dead.trees.ha <-
dead.trees[,
list(num.dead.trees.ha =
weighted.mean(num.dead.trees.ha, w = plot.size.m2)),
by = period]
my.plots$plot4 <-
ggplot(num.dead.trees.ha, aes(period, num.dead.trees.ha, group = 1)) +
xlab('Simulation period') +
ylim(0,NA)+
geom_line() + ylab("Average number of dead trees per ha")
}
if (with.legend == FALSE){
my.plots$plot4 <- my.plots$plot4 + theme(legend.position="none")
}
data.summary$num.dead.trees.ha <- num.dead.trees.ha
}
## plot5 --number of trees harvest
if (5 %in% plots){
## add zeros
rem.trees <- expand.grid(plot.id = plot.data$plot.id,
period = colnames(sitrees.res$live$data$dbh.mm))
rem.trees$num.removed.trees.ha <- 0
rem.trees <- data.table(rem.trees)
rem.trees[num.removed.trees.ha, num.removed.trees.ha := i.x,
on = list(plot.id, period)]
rem.trees[plot.data, plot.size.m2 := i.plot.size.m2, on = 'plot.id']
rem.trees$plot.id <- as.factor(rem.trees$plot.id)
if (by.stand){
my.plots$plot5 <-
ggplot(rem.trees, aes(x = period, y = num.removed.trees.ha, color = plot.id,
group = plot.id)) + xlab('Simulation period') +
geom_line() + ylab( "removed trees per ha") +ylim(0,NA)+
guides(color = guide_legend(title = "Plot ID"))
} else {
num.removed.trees.ha <-
rem.trees[, list(num.removed.trees.ha = weighted.mean(num.removed.trees.ha,
w = plot.size.m2)),
by = period]
my.plots$plot5 <-
ggplot(num.removed.trees.ha,
aes(period, num.removed.trees.ha, group = 1)) + geom_line() +ylim(0,NA)+
xlab('Simulation period') +
ylab("Average number of removed trees per ha")
}
if (with.legend == FALSE){
my.plots$plot5 <- my.plots$plot5 + theme(legend.position="none")
}
data.summary$num.removed.trees.ha <- num.removed.trees.ha
}
if (plot) {
if (plot.all.together & length(my.plots) > 1){
if (by.stand){
if(1 %in% plots) SBA.m2.ha <- SBA.m2.ha[,c("period", "t", "id")] else {
SBA.m2.ha <- NA}
if(2 %in% plots) stems.ha <- stems.ha[,c("period", "t", "id")] else {
stems.ha <- NA}
if(3 %in% plots) heights.10 <- heights.10[,c("period", "t", "id")] else {
heights.10 <- NA}
if(4 %in% plots) {
num.dead.trees.ha <- num.dead.trees.ha[,c("period", "x", "plot.id")]
} else num.dead.trees.ha <- NA
if(5 %in% plots) {
num.removed.trees.ha <- num.removed.trees.ha[,c("period", "x",
"plot.id")]
} else num.removed.trees.ha <- NA
xx <- list(
SBA.m2.ha = SBA.m2.ha,
stems.ha = stems.ha,
heights.10 = heights.10,
num.dead.trees.ha = num.dead.trees.ha,
num.removed.trees.ha = num.removed.trees.ha
)
xx[!1:5 %in% plots] <- NULL
xx <- lapply(xx, FUN = function(x) {
names(x) <- c("period", "t", "id")
return(x)
})
xx <- rbindlist(xx, idcol = 'id.plot')
xx$plot.id <- as.factor(xx$id)
xx[id.plot == 'SBA.m2.ha', id.plot := 'SBA~(m^2/ha)']
xx[id.plot == 'stems.ha', id.plot := 'Stems/ha']
xx[id.plot == 'heights.10', id.plot := 'Height~of~the~10~tallest~trees~(dm)']
xx[id.plot == 'num.dead.trees.ha', id.plot := 'Dead~trees~(stems/ha)']
xx[, id.plot := factor(id.plot, levels = c('SBA~(m^2/ha)', 'Stems/ha',
'Height~of~the~10~tallest~trees~(dm)',
'Dead~trees~(stems/ha)'))]
my.plots[['all_plots_figure']] <-
ggplot(xx, aes(period, t, group = plot.id, col = plot.id)) +
geom_line() +
ylim(0,NA) +
facet_wrap(~ id.plot, scales = "free",
labeller = label_parsed) +
ylab("") + xlab ("Simulation period") +
guides(color = guide_legend(title = "Plot ID"))
} else {
if(1 %in% plots) SBA.m2 <- SBA.m2.ha[,c("period", "t")] else {
SBA.m2.ha <- NA}
if(2 %in% plots) stems <- stems.ha[,c("period", "t")] else {
stems.ha <- NA}
if(3 %in% plots) heights.10 <- heights.10[,c("period", "t")] else {
heights.10 <- NA}
if(4 %in% plots) {
num.dead.trees.ha <- num.dead.trees.ha[,c("period", "num.dead.trees.ha")]
} else num.dead.trees.ha <- NA
if(5 %in% plots) {
num.removed.trees.ha <- num.removed.trees.ha[,c("period", "num.removed.trees.ha")]
} else num.removed.trees.ha <- NA
xx <- list(
SBA.m2.ha = SBA.m2,
stems.ha = stems,
heights.10 = heights.10,
num.dead.trees.ha = num.dead.trees.ha,
num.removed.trees.ha = num.removed.trees.ha
)
xx[!1:5 %in% plots] <- NULL
xx <- lapply(xx, FUN = function(x) {
if (ncol(x) == 2) x$id <- 1
names(x) <- c("period", "t", "id")
return(x)
})
xx <- rbindlist(xx, idcol = 'id.plot')
xx$plot.id <- as.factor(xx$id)
xx[id.plot == 'SBA.m2.ha', id.plot := 'Average~SBA~(m^2/ha)']
xx[id.plot == 'stems.ha', id.plot := 'Average~stems/ha']
xx[id.plot == 'heights.10', id.plot := 'Average~height~of~the~10~tallest~trees~(dm)']
xx[id.plot == 'num.dead.trees.ha', id.plot := 'Average~dead~trees~(stems/ha)']
xx[, id.plot := factor(id.plot, levels = c('Average~SBA~(m^2/ha)',
'Average~stems/ha',
'Average~height~of~the~10~tallest~trees~(dm)',
'Average~dead~trees~(stems/ha)'))]
my.plots[['all_plots_figure']] <-
ggplot(xx, aes(period, t, group = plot.id, col = plot.id)) +
geom_line()+ ylim (0,NA) +
facet_wrap(~id.plot, scales = "free",
labeller = label_parsed) +
ylab("") + xlab ("Simulation period") +
theme(legend.position="none")
print(my.plots[['all_plots_figure']])
}## end if by.stand
} else { ## if plot.all.together
for (i in 1:length(my.plots)){
print(my.plots[[i]])
cat ("Press [enter] to continue")
line <- readline()
}
} ## end else (if plot.all.together)
}## end if plot
attr(my.plots,'data') <- data.summary
invisible(my.plots)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.