---
title: "Overview simulated data"
author: "CB"
date: "`r format(Sys.time(), '%d %m, %Y')`"
output:  html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

```


```{r, include=FALSE}
path <- paste0(getwd(),"/")

# source(paste0(path,"Main.R"))

### DEFINE simulation variant ###

# only the following three, special set aside simulations are bind
#
# sim_variant:  "Without"
#               "CC45"
#               "CC85"


sim_variant = "CC85"


library(tidyr)
library(dplyr)
library(ggplot2)
library(knitr)
library(kableExtra)


#### Load the data

# Simulation results of all management regimes, SA includes DW extraction
rslt <- read.csv(paste0(path, "output/rslt_", sim_variant, "_all.csv" ), sep = ";", header = TRUE, stringsAsFactors = FALSE)

# Simulation results only for set aside without DW extraction
rslt_SA <- read.csv(paste0(path, "output/rslt_", sim_variant,"_SA_all.csv" ), sep = ";", header = TRUE, stringsAsFactors = FALSE)

rslt <- rslt %>%  rbind(rslt_SA)


```
### Simulated Climate Change: `r sim_variant `

### Names of the restructured GPKGs
```{r, echo=TRUE}
gpkg <- unique(rslt$gpkg)
print(gpkg)

```

### Simulated branching_groups

All simulated branching_groups must be listed in the file /params/regimes.csv !!
Based on this the regime names are merged to the results.
```{r, echo = FALSE}
# simulated_regimes <- unique(rslt$regime)
# print(simulated_regimes)

sim_branching_group <- unique(rslt$branching_group)
print(sim_branching_group)

print(unique(rslt$regime))

```


### Number of simulated stands per GPKG and their size
```{r, echo = TRUE}
stands <- rslt %>% 
  group_by(gpkg) %>% 
  summarise(simulated_stands = n_distinct(id),
            min_stand_size = min(AREA),
            max_stand_size = max(AREA),
            mean_size = mean(AREA))

kable(stands) %>% kable_styling()
```

### Development of average stand volume V (m3/ha) 
#### for certain regimes

Volume development of "SA" and "SA_DWextract" should be identical
```{r}
meanV <- rslt[rslt$regime %in% c("BAU", "SA", "CCF_2", "BAUwGTR", "SA_DWextract"), ] %>% 
   group_by(year, regime, gpkg) %>% 
   summarise(meanV = mean(V)) %>% 
   ggplot(aes(year, meanV)) + 
   geom_line( aes(color = regime)) +
   theme(axis.text.x = element_text(angle = 90)) +
   scale_x_continuous(breaks=c(2016, 2026, 2036, 2046, 2056, 2066, 2076, 2086, 2096, 2106)) +
   scale_y_continuous(limit = c(0,700)) +
   facet_wrap(. ~gpkg)
 plot(meanV)

```

#### only for CCF regimes
```{r}
meanV_CCF <- rslt[rslt$regime %in% c("CCF_1", "CCF_2", "CCF_3", "CCF_4"), ] %>% 
   group_by(year, regime, gpkg) %>% 
   summarise(meanV = mean(V)) %>% 
   ggplot(aes(year, meanV)) + 
   geom_line( aes(color = regime)) +
   theme(axis.text.x = element_text(angle = 90)) +
   scale_x_continuous(breaks=c(2016, 2026, 2036, 2046, 2056, 2066, 2076, 2086, 2096, 2106)) +
   # scale_y_continuous(limit = c(0,700)) +
   facet_wrap(. ~gpkg)
 plot(meanV_CCF)

```

### Development of average h_dom for certain regimes
```{r}
meanH_dom <- rslt[rslt$regime %in% c("BAU", "SA", "CCF_2", "BAUwGTR", "SA_DWextract") , ] %>% 
   group_by(year, regime, gpkg) %>% 
   summarise(meanH_dom = mean(H_dom, na.rm = TRUE )) %>% 
   ggplot(aes(year, meanH_dom)) + 
   geom_line( aes(color = regime)) +
   theme(axis.text.x = element_text(angle = 90)) +
   scale_x_continuous(breaks=c(2016, 2026, 2036, 2046, 2056, 2066, 2076, 2086, 2096, 2106)) +
   facet_wrap(. ~gpkg)
 plot(meanH_dom)
```


 
### What is the maximum H_dom per watershed, and which regime is causing it

```{r, warning=FALSE, message=FALSE}
require(data.table)
dt <- data.table(rslt)
max <- dt[ , max(H_dom, na.rm = TRUE ), by = gpkg]
max <- max %>% rename(H_dom = V1 )

maxH_dom <- rslt %>% 
   semi_join(max, by = c("gpkg","H_dom")) %>% 
   select(gpkg, id, year, regime, H_dom)

kable(maxH_dom) %>% kable_styling()
 
```


### Development of average V_total_deadwood (m3/ha) for certain regimes
```{r}
meanDW <- rslt[rslt$regime %in% c("BAU", "SA", "CCF_2", "BAUwGTR", "SA_DWextract") , ] %>% 
   group_by(year, regime, gpkg) %>% 
   summarise(mean_DW = mean(V_total_deadwood)) %>% 
   ggplot(aes(year, mean_DW)) + 
   geom_line( aes(color = regime)) +
   theme(axis.text.x = element_text(angle = 90)) +
   scale_x_continuous(breaks=c(2016, 2026, 2036, 2046, 2056, 2066, 2076, 2086, 2096, 2106)) +
   facet_wrap(. ~gpkg)
 plot(meanDW)

```


### Development of average CARBON_STORAGE (kg/ha) for certain regimes
```{r}
meanCS <- rslt[rslt$regime %in% c("BAU", "SA", "CCF_2", "BAUwGTR")  , ] %>% 
   group_by(year, regime, gpkg) %>% 
   summarise(mean_CS = mean(CARBON_STORAGE, na.rm = TRUE )/1000) %>% 
   ggplot(aes(year, mean_CS)) + 
   geom_line( aes(color = regime)) +
   theme(axis.text.x = element_text(angle = 90)) +
   scale_x_continuous(breaks=c(2016, 2026, 2036, 2046, 2056, 2066, 2076, 2086, 2096, 2106)) +
   facet_wrap(. ~gpkg)
 plot(meanCS)

```


### Average harvested timber volume (m3/ha) for certain regimes 

Cash flow = The sum of all revenues and costs for a specific forest stand
```{r}
meanHarvested_V <- rslt[rslt$regime %in%c("BAU", "SA", "CCF_2", "BAUwGTR") , ] %>% 
   group_by(year, regime, gpkg) %>% 
   mutate(Harvested_V = ifelse(is.na(Harvested_V), 0, Harvested_V)) %>% 
   summarise(meanHarvested_V = mean(Harvested_V), na.rm = TRUE) %>% 
   ggplot(aes(year, meanHarvested_V)) + 
   geom_line( aes(color = regime)) +
   theme(axis.text.x = element_text(angle = 90)) +
   scale_x_continuous(breaks=c(2016, 2026, 2036, 2046, 2056, 2066, 2076, 2086, 2096, 2106)) +
   facet_wrap(. ~gpkg)
 plot(meanHarvested_V)

```


### Average cash flow (Euro/ha) for certain regimes 

Cash flow = The sum of all revenues and costs for a specific forest stand
```{r}
meanCash <- rslt[rslt$regime %in%c("BAU", "SA", "CCF_2", "BAUwGTR") , ] %>% 
   group_by(year, regime, gpkg) %>% 
   mutate(cash_flow = ifelse(is.na(cash_flow), 0, cash_flow)) %>% 
   summarise(meanCash = mean(cash_flow), na.rm = TRUE) %>% 
   ggplot(aes(year, meanCash)) + 
   geom_line( aes(color = regime)) +
   theme(axis.text.x = element_text(angle = 90)) +
   scale_x_continuous(breaks=c(2016, 2026, 2036, 2046, 2056, 2066, 2076, 2086, 2096, 2106)) +
   facet_wrap(. ~gpkg)
 plot(meanCash)

```
