add new SA simulations without deadwood extraction

Main.R

 #####
 #
-# CB
-# 2019-05-22
+# Load simulated Watershed data of SIMO and create CSV.files for further analysis
+# 2020-02-21
 #
 #####
 
 
+
 library(dplyr)
 
 
@@ -14,22 +15,27 @@ path <- paste0(getwd(),"/")
 
 
 ### load some parameter
-# Management regimes and their abbreviation 
+# Management regimes and their abbreviation, they are merged to the data by the branching group (script loadDB.R)
 regime <- read.csv(paste0(path, "params/regimes.csv"), sep = ",", stringsAsFactors = FALSE)
 
 
-### Simulation variant
-# specify which one:  "CC45" climate change with RCP scenraio 4.5 
-#                     "CC85" climate change with RCP scenario 8.5 
-#                     "without" no climate change
-sim_variant <- "without"  
+### Simulation variant that will be load
+# specify which one:  "CC45"        climate change with RCP scenraio 4.5 
+#                     "CC45_SA"     RCP 4.5 and set aside without deadwood extraction (no other management regimes!)
+#                     "CC85"        climate change with RCP scenario 8.5 
+#                     "CC85_SA"     RCP 8.5 and set aside without deadwood extraction (no other management regimes!)
+#                     "without"     no climate change
+#                     "without_SA"  no climate change and set aside without deadwood extraction (no other management regimes!)
+
 
+sim_variant <- "CC45"  
 
-### Define the names of the SQL databases (SIMO-output) that will imported
+
+### Define the names of the databases (SIMO-output for 10 watersheds) that will be imported
 # It is used in the skripts "structure_SIMO_rslDB_FBE.R" and "loadDB.R"
 db_names <- c("MV_Hartola", 
               "MV_Kitee", 
-              "MV_Korsnas", 
+              "MV_Korsnas",
               "MV_Parikkala",
               "MV_Pori",
               "MV_Pyhtaa",
@@ -39,12 +45,35 @@ db_names <- c("MV_Hartola",
               "MV_Voyri")
 
 
-### Restructure the SQL database: querry that creates the final UNIT table, which contains indicators over time under regimes
+### Restructure the SQL database. 
+# The query creates a table called UNIT, which contains indicators over time and under management regimes
 #
-# !!! Only needed if the DB is read/loaded for the first time. THis may take some time  !!!
-# !!! For your downloaded watershed data this is already done. !!!
+# !!! Only needed if the DB is loaded for the first time (this may take some time):     "first_load = TRUE"
+# !!! For the downloaded watershed data this is already done:                           "first_load = FALSE"
+
+
+first_load = FALSE
+
+
+### the following lines do not need any changes ###
+if(first_load == TRUE){
 
-# source(paste0(path, "structure_SIMO_rslDB_FBE.R")) # Afterwards, you can exclude this line with: #
+  # If one of the follwing simulation variants is read ...
+  if(sim_variant %in% c("CC45", "CC85", "without")) {
+    
+    # Run the script with the SQL query for all management regimes
+    source(paste0(path, "structure_SIMO_rslDB_FBE.R"))  
+  
+    } else {
+    
+      # Otherwise, run the script witht the SQL query for Set aside without deadwood extraction (...SA)
+      # This needs a different SQL query, since the database does not contain harvest information
+      source(paste0(path, "structure_SIMO_rslDB_SA.R")) # Query 
+      
+      }
+
+}
+###
 
 
 ### Import the restructured SIMO data (from UNIT table) in the R-environment
@@ -72,13 +101,31 @@ columns <-  paste0("id,
                    CARBON_STORAGE")
 
 
-### !!! Needed if indicator columns are selected for the first time, since dataframes are also stored as csv-file under "output/".  !!!
+### !!! CSV files are stored unter output
+#
+# If data/columns have already been loaded and saved as csv file  under ../output:                csv_exist = TRUE
+# If data is loaded for the firest time OR "new columns have to be loaded" (takes some time):     csv_exits = FALSE 
+
+
+csv_exist = FALSE
+
+
+### the following lines do not need any changes ###
+if(csv_exist == FALSE){
+  
+  # If FALSE run the following script 
+  source(paste0(path, "loadDB.R")) 
+ 
+  } else {
+  
+    # If TRUE read the CSV file that contains all watersheds (GPKGs)
+    rslt <- read.csv(paste0(path, "output/rslt_", sim_variant, "_all.csv" ), sep = ";", header = TRUE, stringsAsFactors = FALSE)
+
+    }
+###
 
-# source(paste0(path, "loadDB.R")) # Afterwards, you can also exclude this line with: #
 
-# ... and the single csv can later also be loaded into the R-environment with:
 
-rslt <- read.csv(paste0(path, "output/rslt_", sim_variant, "_all.csv" ), sep = ";", header = TRUE, stringsAsFactors = FALSE)
 
 
 

loadDB.R

@@ -21,11 +21,13 @@ library(dplyr)
 
 for (name in db_names){
   
+  # name = "MV_Hartola"
   db <- dbConnect(dbDriver("SQLite"), dbname = paste0(path,"input_data/simulated_", sim_variant, "_", name, "_rsu.db"))
   rsl  <- dbGetQuery( db, paste0("select ", columns, " from UNIT"))
   rsl$gpkg <- name
   dbDisconnect(db)
   
+  ### Add the abbreviation of the regimes
   rsl <- rsl %>% 
     left_join(regime, by = "branching_group", all.x = TRUE)
   
@@ -37,6 +39,15 @@ for (name in db_names){
   rsl <- rsl %>%
     anti_join(error_stands, by = c("id", "gpkg"))
   
+  
+  ### Rename set aside scenario if it considers deadwood extraction
+  if(sim_variant %in% c("CC45", "CC85", "without")) {
+    
+   rsl <- rsl %>% mutate(regime = ifelse(regime %in% "SA", "SA_DWextract", regime))
+    
+  }
+  
+  
   write.table(rsl, paste0(path, "output/rsl_",sim_variant,"_",name,".csv" ), sep = ";", row.names = F, col.names = TRUE)
   
   assign( paste("rsl", name, sep="_"), rsl)
@@ -59,10 +70,19 @@ for (name in db_names){
   rm(db, rsl) 
 }  
 
+### Add the abbreviation of the regimes
 rslt <- rslt %>% 
   left_join(regime, by= "branching_group", all.X = TRUE)
 
 
+### Rename set aside scenario if it considers deadwood extraction
+if(sim_variant %in% c("CC45", "CC85", "without")) {
+  
+  rslt <- rslt %>% mutate(regime = ifelse(regime %in% "SA", "SA_DWextract", regime))
+  
+}
+
+
 ### Filter those stands that caused errors during the simulation with SIMO
 # import csv file that contains the error stands
 error_stands <- read.csv(paste0(path, "params/errors_watersheds.csv"), sep = ",", header = TRUE, stringsAsFactors = FALSE)
@@ -77,19 +97,92 @@ write.table(rslt, paste0(path, "output/rslt_", sim_variant, "_all.csv" ), sep =
 
 
 
+##################################################################
+##################################################################
 
 
 
+# # First extraction of wind simulated Korsnas data for MSc
+# 
+# sim_variant <- "without_SA"
+# db_names <- c("MV_Korsnas_Wind_1") #  , "MV_Korsnas_Wind_2", "rsu_example2"
+# 
+# 
+# columns <-  paste0("id,
+#                    year,
+#                    branch,
+#                    branch_desc,
+#                    branching_group,
+#                    Age,
+#                    area,
+#                    cash_flow,
+#                    V_total_deadwood,
+#                    BA,
+#                    V,
+#                    N,
+#                    H_dom,
+#                    D_gm,
+#                    Harvested_V,
+#                    Biomass,
+#                    income_biomass,
+#                    CARBON_STORAGE")
+#                     
+# 
+#                     # Harvested_V_log,
+#                     # Harvested_V_log_under_bark,
+#                     # Harvested_V_pulp,
+#                     # Harvested_V_pulp_under_bark
+#                     # V_log,   
+#                     # H_gm,
+#                     # SINCE_DRAINAGE,
+#                     # DRAINAGE_STATUS,
+#                     # regen_age,
+#                     # SINCE_DRAINAGE_ORIG,
+#                     # V_pulp,
+#                     # SOIL_CLASS
+#                     # Carb_flux_nat_wd_nrg,
+#                     # Carbon_flux_natural_rm_wind")
+# 
+# 
+# 
+# rslt <- NULL
+# 
+# for (name in db_names){
+#   db <- dbConnect(dbDriver("SQLite"), dbname = paste0(path,"input_data/without_SA_5_stands/simulated_", sim_variant, "_" , name, ".db"))
+#   rsl  <- dbGetQuery( db, paste0("select ", columns, " from UNIT"))
+#   rsl$gpkg <- name
+#   dbDisconnect(db)
+#   rslt <- rbind(rslt, rsl)
+#   rm(db, rsl)
+# }
+# 
+# rslt <- rslt %>%
+#   left_join(regime, by= "branching_group", all.X = TRUE)
+# 
+# ### write table
+# write.table(rslt, paste0(path, "output/rsl_without_SA_MV_Kitee.csv" ), sep = ";", row.names = F, col.names = TRUE)
+# 
+# rslt <- read.csv(paste0(path, "output/rsl_without_SA_MV_Kitee.csv" ), sep = ";", header = TRUE, stringsAsFactors = FALSE)
 
-# ##### Getting all Colnames from a SQL querry
+
+
+
+# ##################################################################
+# ##################################################################
+# 
 # 
+# ##### Getting all Colnames from a SQL querry
+#
+#
 # # was used to create List of SIMO output -> What indicators can be calculate for ES
-#  
-# db <- dbConnect(dbDriver("SQLite"), dbname = paste0(path,"input_data/JENNI.db"))
+# 
+# db <- dbConnect(dbDriver("SQLite"), dbname = paste0(path,"input_data/test/simulated_without_rsu_example2.db"))
 # rsl  <- dbGetQuery(db, 'select * from UNIT')
 # dbDisconnect(db)
 # 
-# colnames_FBE <- colnames(rsl) 
-# noquote(colnames_FBE)
+# colnames_without <- colnames(rsl)
+# noquote(colnames_without)
 # write.table(colnames_FBE, paste0(path,"temp/colnames_FBE.csv"), row.names = FALSE)
-
+# 
+# 
+# onlyinCC <- colnames_CC45[!colnames_CC45 %in% colnames_without]

params/SQL_Maiju_SA.sql

+DROP TABLE IF EXISTS opers2;
+DROP TABLE IF EXISTS opers3;
+DROP TABLE IF EXISTS max_v;
+
+CREATE TABLE max_v AS SELECT comp_unit.id AS id,MAX(comp_unit.V) AS max_v FROM comp_unit GROUP BY comp_unit.id;
+ 
+DROP TABLE IF EXISTS unit;
+
+Create Table UNIT AS SELECT u.*,(select max(stratum.H_dom) From stratum where stratum.data_id = u.data_id) as H_dom, 
+(select max(stratum.D_gm) From stratum where stratum.data_id = u.data_id) as D_gm, 
+(select sum(stratum.N) From stratum where stratum.data_id = u.data_id and D_gm >40) as N_where_D_gt_40,
+(select sum(stratum.N) From stratum where stratum.data_id = u.data_id and D_gm <=40 and D_gm > 35) as N_where_D_gt_35_lt_40,
+(select sum(stratum.N) From stratum where stratum.data_id = u.data_id and D_gm <=35 and D_gm > 30) as N_where_D_gt_30_lt_35,
+(select sum(stratum.V) From stratum where stratum.data_id = u.data_id and D_gm >40) as V_where_D_gt_40,
+(select sum(stratum.V) From stratum where stratum.data_id = u.data_id and D_gm <=40 and D_gm > 35) as V_where_D_gt_35_lt_40,
+(select sum(stratum.V) From stratum where stratum.data_id = u.data_id and D_gm <=35 and D_gm > 30) as V_where_D_gt_30_lt_35,
+(select sum(stratum.V) From stratum where stratum.data_id = u.data_id and SP = 5) as V_populus,
+(select sum(stratum.V) From stratum where stratum.data_id = u.data_id and SP = 6) as V_Alnus_incana,
+(select sum(stratum.V) From stratum where stratum.data_id = u.data_id and SP = 7) as V_Alnus_glutinosa,
+(select sum(stratum.V) From stratum where stratum.data_id = u.data_id and SP = 8) as V_o_coniferous,
+(select sum(stratum.V) From stratum where stratum.data_id = u.data_id and SP = 9) as V_o_decidious,
+ l.data_date, b.branch, b.branch_desc, b.branching_group, 0 as income, 0 as cash_flow,  0 as clearcut, 0 as Harvested_V_log, 
+ 0 as Harvested_V_pulp,  0 as Harvested_V,  m.max_v, 0 as income_log, 0 as income_pulp,
+0 as income_log_change,  0 as income_pulp_change, 0 as income_biomass, 0 as Biomass
+ FROM comp_unit u, data_link l
+left outer join branch_desc b on l.branch = b.branch and l.id = b.id 
+cross join max_v m on l.id = m.id
+WHERE u.data_id=l.data_id and max_v <1000
+ORDER BY u.id, l.branch, l.data_date;
+
+DROP TABLE IF EXISTS BAU;
+DROP TABLE IF EXISTS BAU_10;DROP TABLE IF EXISTS BAU_15;DROP TABLE IF EXISTS BAU_30;DROP TABLE IF EXISTS BAU_5;DROP TABLE IF EXISTS BAU_m5;DROP TABLE IF EXISTS BAUwGTR;DROP TABLE IF EXISTS BAUwoT;DROP TABLE IF EXISTS BAUwoT_10;DROP TABLE IF EXISTS BAUwoT_m20;DROP TABLE IF EXISTS BAUwT;DROP TABLE IF EXISTS BAUwT_10;DROP TABLE IF EXISTS BAUwT_15;DROP TABLE IF EXISTS BAUwT_30;DROP TABLE IF EXISTS BAUwT_5;DROP TABLE IF EXISTS BAUwT_GTR;DROP TABLE IF EXISTS BAUwt_m5;DROP TABLE IF EXISTS CCF_1;
+DROP TABLE IF EXISTS CCF_2;
+DROP TABLE IF EXISTS CCF_3;
+DROP TABLE IF EXISTS CCF_4;DROP TABLE IF EXISTS SA;
+
+CREATE TABLE SA AS
+SELECT u.*
+FROM unit u
+where u.branching_group is Null;
+

quick_summary.Rmd

@@ -14,13 +14,38 @@ knitr::opts_chunk$set(echo = TRUE)
 ```{r, include=FALSE}
 path <- paste0(getwd(),"/")
 
-source(paste0(path,"Main.R"))
+# 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 = "without"
+
 
 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 `
 
@@ -42,6 +67,8 @@ Based on this the regime names are merged to the results.
 sim_branching_group <- unique(rslt$branching_group)
 print(sim_branching_group)
 
+print(unique(rslt$regime))
+
 ```
 
 
@@ -59,8 +86,10 @@ 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"), ] %>% 
+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)) + 
@@ -82,7 +111,7 @@ meanV_CCF <- rslt[rslt$regime %in% c("CCF_1", "CCF_2", "CCF_3", "CCF_4"), ] %>%
    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)) +
+   # scale_y_continuous(limit = c(0,700)) +
    facet_wrap(. ~gpkg)
  plot(meanV_CCF)
 
@@ -91,7 +120,7 @@ meanV_CCF <- rslt[rslt$regime %in% c("CCF_1", "CCF_2", "CCF_3", "CCF_4"), ] %>%
 
 ### Development of average V_total_deadwood (m3/ha) for certain regimes
 ```{r}
-meanDW <- rslt[rslt$regime %in% c("BAU", "SA", "CCF_2", "BAUwGTR") , ] %>% 
+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)) + 

structure_SIMO_rslDB_FBE.R

@@ -121,7 +121,11 @@ create_table_UNIT <- 'Create Table UNIT AS SELECT u.*,
                       ORDER BY u.id, l.branch, l.data_date'
 
 
-# For each "db_names", defined in main.R ...
+
+
+##### For each "db_names", defined in main.R ...
+
+
 
 for (name in db_names){
   # Connect to the database
@@ -148,3 +152,48 @@ for (name in db_names){
 }
 
 rm(create_table_max_v, create_table_OPERS2, create_table_OPERS3, create_table_UNIT)
+
+
+
+
+##################################################################
+##################################################################
+
+
+
+
+# # First extraction of wind simulated Korsnas data for MSc
+# 
+# sim_variant <- "without_SA"
+# db_names <- c("MV_Korsnas_Wind_1") #  , "MV_Korsnas_Wind_2", "rsu_example2"
+#
+# 
+# 
+# for (name in db_names){
+#   # Connect to the database
+#   con <- dbConnect(dbDriver("SQLite"), dbname = paste0(path,"input_data/simulated_", sim_variant,"_" ,name , "_rsu.db"))
+#   # con <- dbConnect(dbDriver("SQLite"), dbname = paste0(path,"input_data/test/simulated_", sim_variant,"_" ,name , ".db"))
+#   
+#   # If the following tables already exist, for which the query is defined, remove them
+#   tab_to_delete <- c("OPERS2", "OPERS3", "max_v", "UNIT")
+#   
+#   for(i in tab_to_delete){
+#     if(dbExistsTable(con, i)) {dbRemoveTable(con, i)}
+#   }
+#   
+#   # Run the Queries and create the final table "UNIT", 
+#   # which contains the development of all stands and indicators under the simulated management regimes 
+#   query_to_run <- c(create_table_max_v, create_table_OPERS2, create_table_OPERS3, create_table_UNIT)
+#   
+#   for(i in query_to_run){
+#     dbExecute(con, i)
+#   }
+#   
+#   dbDisconnect(con)
+#   
+#   rm(query_to_run, tab_to_delete, con)
+# }
+# 
+# rm(create_table_max_v, create_table_OPERS2, create_table_OPERS3, create_table_UNIT)
+
+

structure_SIMO_rslDB_SA.R

+#####
+#
+# Restructer the SIMO output in the SQLite database (final "UNIT" table)
+# 
+# !!! ONLY FOR SetAside simulations
+# 
+# 2020-02-19
+#
+#####
+
+
+# Define the SQL queries:
+# The following queries are based on a SQL-script from K. Eyvindson (see folder params)
+#  -  SQL_Maiju.sql
+# They create a final table "UNIT" in the SIMO output database. 
+# Table UNIT contains the development of indicators over the time under the different management regimes.
+
+
+
+## load libraries
+library(RSQLite)
+
+
+
+# Queries:
+create_table_max_v <- 'CREATE TABLE max_v AS SELECT comp_unit.id AS id,
+                       MAX(comp_unit.V) AS max_v FROM comp_unit GROUP BY comp_unit.id'
+
+
+create_table_UNIT <-    'Create Table UNIT AS SELECT u.*,(select max(stratum.H_dom) From stratum where stratum.data_id = u.data_id) as H_dom, 
+                        (select max(stratum.D_gm) From stratum where stratum.data_id = u.data_id) as D_gm, 
+                        (select sum(stratum.N) From stratum where stratum.data_id = u.data_id and D_gm >40) as N_where_D_gt_40,
+                        (select sum(stratum.N) From stratum where stratum.data_id = u.data_id and D_gm <=40 and D_gm > 35) as N_where_D_gt_35_lt_40,
+                        (select sum(stratum.N) From stratum where stratum.data_id = u.data_id and D_gm <=35 and D_gm > 30) as N_where_D_gt_30_lt_35,
+                        (select sum(stratum.V) From stratum where stratum.data_id = u.data_id and D_gm >40) as V_where_D_gt_40,
+                        (select sum(stratum.V) From stratum where stratum.data_id = u.data_id and D_gm <=40 and D_gm > 35) as V_where_D_gt_35_lt_40,
+                        (select sum(stratum.V) From stratum where stratum.data_id = u.data_id and D_gm <=35 and D_gm > 30) as V_where_D_gt_30_lt_35,
+                        (select sum(stratum.V) From stratum where stratum.data_id = u.data_id and SP = 5) as V_populus,
+                        (select sum(stratum.V) From stratum where stratum.data_id = u.data_id and SP = 6) as V_Alnus_incana,
+                        (select sum(stratum.V) From stratum where stratum.data_id = u.data_id and SP = 7) as V_Alnus_glutinosa,
+                        (select sum(stratum.V) From stratum where stratum.data_id = u.data_id and SP = 8) as V_o_coniferous,
+                        (select sum(stratum.V) From stratum where stratum.data_id = u.data_id and SP = 9) as V_o_decidious,
+                        l.data_date, b.branch, b.branch_desc, b.branching_group, 0 as income, 0 as cash_flow,  0 as clearcut, 0 as Harvested_V_log, 
+                        0 as Harvested_V_pulp,  0 as Harvested_V,  m.max_v, 0 as income_log, 0 as income_pulp,
+                        0 as income_log_change,  0 as income_pulp_change, 0 as income_biomass, 0 as Biomass
+                        FROM comp_unit u, data_link l
+                        left outer join branch_desc b on l.branch = b.branch and l.id = b.id 
+                        cross join max_v m on l.id = m.id
+                        WHERE u.data_id=l.data_id and max_v <1000
+                        ORDER BY u.id, l.branch, l.data_date'
+
+
+##### For each "db_names", defined in main.R ...
+
+
+# sim_variant <- "without_SA"
+# db_names <- c("MV_Hartola", 
+#               "MV_Kitee", 
+#               "MV_Korsnas", 
+#               "MV_Parikkala",
+#               "MV_Pori",
+#               "MV_Pyhtaa",
+#               "MV_Raasepori",
+#               "MV_Simo",
+#               "MV_Vaala",
+#               "MV_Voyri")
+
+
+for (name in db_names){
+  
+  # name = "MV_Hartola"
+  
+  # Connect to the database
+  con <- dbConnect(dbDriver("SQLite"), dbname = paste0(path,"input_data/simulated_", sim_variant,"_" ,name , "_rsu.db"))
+  
+  # If the following tables already exist, for which the query is defined, remove them
+  tab_to_delete <- c("OPERS2", "OPERS3", "max_v", "UNIT")
+  
+  for(i in tab_to_delete){
+    if(dbExistsTable(con, i)) {dbRemoveTable(con, i)}
+  }
+  
+  # Run the Queries and create the final table "UNIT", 
+  # which contains the development of all stands and indicators under the simulated management regimes 
+  query_to_run <- c(create_table_max_v, create_table_UNIT)
+  
+  for(i in query_to_run){
+    dbExecute(con, i)
+  }
+  
+  dbDisconnect(con)
+  
+  rm(query_to_run, tab_to_delete, con)
+}
+
+rm(create_table_max_v, create_table_UNIT)
+
+
+
+