check_opt_data_rcp0_2_kyle.Rmd

---
title: "check_opt_data"
author: "CB"
date: "7/28/2020"
output: html_document
---


# RCP0 whole Finland - version 7
### (rslt_RCP0_FIN_V7.csv)

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


### Running on cPouta or WIN local
on_cPouta <- TRUE


### Set path to input files

if(on_cPouta == TRUE) {
  # for pouta
  datapath <- "/media/volume/outp_rcp0/"
  path <- "/home/ubuntu/workspace/MultiForest2/" 
  # path <- paste0(getwd(),"/")
  
} else {
  # for Windows
  datapath <- "C:/MyTemp/r_SIMO_output/output/MF_FIN/RCP0_test/"

}



### load libraries
library(dplyr)
library(tidyr)
library(ggplot2)
library(data.table)
library(gridExtra)
library(egg)
library(kableExtra)


```



```{r , echo=FALSE}

poutaxy <- read.csv(paste0(datapath, "rslt_RCP0_FIN_V7_XY.csv"), sep = ";", header = TRUE, stringsAsFactors = FALSE)
pouta0 <- read.csv(paste0(datapath, "rslt_RCP0_FIN_V7.csv"), sep = ";", header = TRUE, stringsAsFactors = FALSE)

# poutaxy <- read.csv(paste0(datapath, "rslt_test_XY.csv"), sep = ";", header = TRUE, stringsAsFactors = FALSE)
# pouta0 <- read.csv(paste0(datapath, "rslt_test.csv"), sep = ";", header = TRUE, stringsAsFactors = FALSE)



pouta <- pouta0 %>%
  dplyr::left_join(poutaxy[,c("id", "year", "regime", "standid", "CARBON_STORAGE")], by = c("id", "year", "regime"))

```


## How many stands?

```{r, echo=TRUE}

length(unique(pouta$id))

```



## Simulated regimes

```{r, echo=TRUE}

unique(pouta$regime)

```



# Check for outliers of certain indicators

Based on summary, the following are looked at due to high max-min values: </br>
- V </br>
- i_Vm3 </br>
- CARBON SINK </br>


## Indicator: V

```{r , echo=TRUE}

boxplot(pouta$V, main="Indicator: V")

V_out <- pouta %>% dplyr::filter(V > 1500) %>% # Standing volume above xx m3/ha
  dplyr::mutate(error = "V")

# how many stands effected
length(unique(V_out$standid))

# ---------------
# which stands
# ---------------
unique(V_out$standid)

# which region
unique(V_out$region)

# which regime
unique(V_out$regime)

```




## Indicator: i_Vm3

```{r , echo=TRUE}

boxplot(pouta$i_Vm3, main="Indicator: i_Vm3")

i_Vm3_out <- pouta %>%  filter(i_Vm3 > 50) %>%  # annual volume increment above xx m3/ha
  mutate(error = "i_Vm3")

# how many stands effected
length(unique(i_Vm3_out$id))


# ---------------
# which stands
# ---------------
unique(i_Vm3_out$standid)

# which region
unique(i_Vm3_out$region)

# which regime
unique(i_Vm3_out$regime)


```



### What are the high values

```{r , echo=TRUE}

i_Vm3_high <- pouta %>% 
  select(c("id", "standid", "year", "regime", "i_Vm3", "V","Harvested_V")) %>%
  filter(id %in% unique(i_Vm3_out$id)) %>% 
  filter(i_Vm3 > 50 ) 

# first 20 rows
i_Vm3_high[1:20,] %>%
  knitr::kable(caption = "FIRST 20 ROWS ONLY !!") %>%
  kable_styling()

```


```{r , echo=TRUE}


## Indicator: CARBON SINK

```{r , echo=TRUE}

plot_carbon <- pouta %>% 
  filter(CARBON_SINK > -5000000)

boxplot(plot_carbon$CARBON_SINK, main="CARBON_SINK without the huge outlier")

carbon_sink_out <- pouta %>% filter(CARBON_SINK < -1000000 | CARBON_SINK > 1000000) %>% # source/sink of more than 1000 t CO2/ha
  mutate(error = "CARBON_SINK")

# how many stands
length(unique(carbon_sink_out$id))

# ---------------
# which stands
# ---------------
unique(carbon_sink_out$standid)

# which regime
unique(carbon_sink_out$regime)

# which region
unique(carbon_sink_out$region)

```



### What are the high values

```{r , echo=TRUE}

carbononly <- pouta  %>%
  select(c("id", "standid", "year", "regime", "V", "Harvested_V" ,"CARBON_STORAGE", "CARBON_SINK")) %>%
  filter(id %in% unique(carbon_sink_out$id)) 

carbononly %>% 
  filter(CARBON_SINK < -1000000 | CARBON_SINK > 1000000) %>% 
  knitr::kable(caption = "ALL ROWS !!") %>%
  kable_styling()

```



### Plot some example stands

```{r , echo=TRUE}

carbonplot <- carbononly %>% 
  filter(regime %in% c("SA","initial_state") & id %in% 99037593 | 
           regime %in% c("BAU_30","initial_state") & id %in% 99034584 |
           regime %in% c("CCF_2", "CCF_1", "CCF_3", "initial_state") & id %in% 99032541) %>% 
  mutate(deltaCarbonstorage = round(CARBON_SINK / 3.67, digits = 2))

# for scenarios BAU_30 and SA
p1 = ggplot(carbonplot[carbonplot$regime %in% c("BAU_30", "SA"),], aes(year, V)) + 
  geom_line() +
  facet_grid(.~regime) +
  theme(axis.text.x = element_text(angle = 90))

p2 = ggplot(carbonplot[carbonplot$regime %in% c("BAU_30", "SA"),], aes(year, CARBON_STORAGE)) + 
  geom_line() +
  facet_grid(.~regime) +
  theme(axis.text.x = element_text(angle = 90)) +
  scale_y_continuous(labels = function(x) format(x, scientific = FALSE))

p3 = ggplot(carbonplot[carbonplot$regime %in% c("BAU_30", "SA"),], aes(year, deltaCarbonstorage)) + 
  geom_line() +
  facet_grid(.~regime) +
  theme(axis.text.x = element_text(angle = 90)) +
  scale_y_continuous(labels = function(x) format(x, scientific = FALSE))

ggarrange(p1, p2, p3,
          top = "BAU_30 99034584 standid: 9932854 - SA 99037593	standid: 10628286")

```


!!! HIGH Carbon sinks due to extremly high simulated Volumes !!!


```{r , echo=TRUE}

# for scenarios CC_2, _1 and _3
p4 = ggplot(carbonplot[carbonplot$regime %in% c("CCF_2", "CCF_1", "CCF_3"),], aes(year, V)) +
  geom_line() +
  facet_grid(.~regime) +
  theme(axis.text.x = element_text(angle = 90))

p5 = ggplot(carbonplot[carbonplot$regime %in% c("CCF_2", "CCF_1", "CCF_3"),], aes(year, CARBON_STORAGE)) +
  geom_line() +
  facet_grid(.~regime) +
  theme(axis.text.x = element_text(angle = 90)) +
  scale_y_continuous(labels = function(x) format(x, scientific = FALSE))

p6 = ggplot(carbonplot[carbonplot$regime %in% c("CCF_2", "CCF_1", "CCF_3"),], aes(year, deltaCarbonstorage)) +
  geom_line() +
  facet_grid(.~regime) +
  theme(axis.text.x = element_text(angle = 90)) +
  scale_y_continuous(labels = function(x) format(x, scientific = FALSE))

ggarrange(p4, p5, p6,
          top = "99032541 - standid 23271031")

```