Merge pull request #154 from piklev/develop

*refactoring* changed the efficiency plots in the calibration reporting

scripts/output/single/reportCEScalib.R

@@ -14,16 +14,17 @@ library(quitte)
 require(lucode)
 require(colorspace)
 
+gdx_name     <- "fulldata.gdx"        # name of the gdx  
+gdx_ref_name <- "input_ref.gdx"       # name of the reference gdx (for policy cost calculation)
 
 if(!exists("source_include")) {
   #Define arguments that can be read from command line
-  output_folder <- "dummy"
-  readArgs("output_folder")
-  scenario<-"remind17_6013_SSP2-tax20-Noaff-CT-rem-5"
-  output_folder <- "remind17_6013_SSP2-tax20-Noaff-CT-rem-5"
-} else {
-  output_folder <- outputdir
-}
+   outputdir <- "output/R17IH_SSP2_postIIASA-26_2016-12-23_16.03.23"     # path to the output folder
+   readArgs("outputdir","gdx_name","gdx_ref_name")
+} 
+gdx      <- path(outputdir,gdx_name)
+gdx_ref  <- path(outputdir,gdx_ref_name)
+if(!file.exists(gdx_ref)) { gdx_ref <- NULL }
 scenario <- getScenNames(outputdir)
 
 #---------------------------------------------------------------------------
@@ -83,9 +84,9 @@ cat("Reading CES calibration output from ",filename,"\n")
 if (file.exists(filename)) {
   CES.cal.report <- read.table(filename, header = TRUE, sep = ",", quote = "\"") %>% 
     as.data.frame()
-} else if (file.exists(path(output_folder,filename))) {
+} else if (file.exists(path(outputdir,filename))) {
   
-  CES.cal.report <- read.table(path(output_folder,filename), header = TRUE, sep = ",", quote = "\"") %>% 
+  CES.cal.report <- read.table(path(outputdir,filename), header = TRUE, sep = ",", quote = "\"") %>% 
     as.data.frame() 
 } else {
   stop("No CES_calibration.csv file found. CES_calibration.csv is normally produced during calibration runs")
@@ -97,7 +98,7 @@ if (file.exists(filename)) {
 #---------------------------------------------------------------------------
 
 
-in_set = readGDX(path(outputdir,"fulldata.gdx"), "in", "sets")
+in_set = readGDX(gdx, "in", "sets")
 
 itr <- getColValues(CES.cal.report,"iteration")
 itr_num <- sort(as.double(setdiff(itr, c("origin","target"))))
@@ -108,10 +109,11 @@ col <- c("#fc0000", "#000000",
          "#bc80bd"#,
          #"#808080"
          )
-
+names(col) <- c("origin", "target", seq(1,length(itr_num)))
 
 
 lns <- c(rep("solid", 2), rep("longdash", length(itr_num)))
+names(lns) <-  c("origin", "target", seq(1,length(itr_num)))
 
 .pf <- list("TE" = c("gastr", "refliq", "biotr", "coaltr","hydro", "ngcc","ngt","pc", "apCarDiT","apCarPeT","apCarElT","dot","gaschp","wind","tnrs"))
 
@@ -202,8 +204,8 @@ for (s in levels(CES.cal.report$scenario)) {
       geom_line() + 
       facet_wrap(~ pf, scales = "free", as.table = FALSE) + 
       expand_limits(y = 0) + 
-      scale_colour_manual(values = col[-2]) + 
-      scale_linetype_manual(values = lns[-2]) + 
+      scale_colour_manual(values = col) + 
+      scale_linetype_manual(values = lns) + 
       ggtitle(paste("prices", r, s)) -> p
     plot(p)
     
@@ -212,17 +214,22 @@ for (s in levels(CES.cal.report$scenario)) {
       filter(scenario == s, 
              t        <= 2100, 
              regi     == r, 
-             variable == "total efficiency") %>% 
+             variable == "total efficiency",
+             iteration != "origin") %>% 
+      group_by(scenario,t,regi,pf,variable) %>%
+      mutate(value = value / value[as.character(iteration) == "1"]) %>% 
+      ungroup() %>% 
       order.levels(pf = getElement(.pf,"structure" )) %>% 
       ggplot(aes(x = t, y = value, colour = iteration, 
                  linetype = iteration)) + 
       geom_line() + 
       facet_wrap(~ pf, scales = "free", as.table = FALSE) + 
-      expand_limits(y = 0) + 
-      scale_colour_manual(values = col[-2]) + 
-      scale_linetype_manual(values = lns[-2]) + 
-      ggtitle(paste("total efficiency", r, s)) -> p
+      scale_colour_manual(values = col) + 
+      scale_linetype_manual(values = lns) + 
+      ggtitle(paste("total efficiency (1 = iteration 1)", r, s)) -> p
     plot(p)
+    
+    
     # plot Putty quantities
     if ( dim(CES.cal.report %>% filter(variable == "quantity_putty"))[1] > 0){
     CES.cal.report %>% 
@@ -254,8 +261,8 @@ for (s in levels(CES.cal.report$scenario)) {
       geom_line() + 
       facet_wrap(~ pf, scales = "free", as.table = FALSE) + 
       expand_limits(y = 0) + 
-      scale_colour_manual(values = col[-2]) + 
-      scale_linetype_manual(values = lns[-2]) + 
+      scale_colour_manual(values = col) + 
+      scale_linetype_manual(values = lns) + 
       ggtitle(paste("prices", r, s)) -> p
     plot(p)
     
@@ -264,16 +271,20 @@ for (s in levels(CES.cal.report$scenario)) {
       filter(scenario == s, 
              t        <= 2100, 
              regi     == r, 
-             variable == "total efficiency putty") %>% 
+             variable == "total efficiency putty",
+             iteration != "origin") %>% 
+      group_by(scenario,t,regi,pf,variable) %>%
+      mutate(value = value / value[as.character(iteration) == "1"]) %>% 
+      ungroup() %>% 
       order.levels(pf = getElement(.pf,"structure" )) %>% 
       ggplot(aes(x = t, y = value, colour = iteration, 
                  linetype = iteration)) + 
       geom_line() + 
       facet_wrap(~ pf, scales = "free", as.table = FALSE) + 
       expand_limits(y = 0) + 
-      scale_colour_manual(values = col[-2]) + 
-      scale_linetype_manual(values = lns[-2]) + 
-      ggtitle(paste("total efficiency", r, s)) -> p
+      scale_colour_manual(values = col) + 
+      scale_linetype_manual(values = lns) + 
+      ggtitle(paste("total efficiency (1 = iteration 1)", r, s)) -> p
     plot(p)
     
     }
@@ -297,8 +308,8 @@ for (s in levels(CES.cal.report$scenario)) {
       geom_line() + 
       facet_wrap(~ pf, scales = "free", as.table = FALSE) + 
       expand_limits(y = 0) + 
-      scale_colour_manual(values = col[-2]) + 
-      scale_linetype_manual(values = lns[-2]) + 
+      scale_colour_manual(values = col) + 
+      scale_linetype_manual(values = lns) + 
       geom_vline(xintercept = 2005) +
       ggtitle(paste("vm_deltaCap", r, s)) -> p
     plot(p)