Commit 9ec30b46 authored by Toni's avatar Toni
Browse files

Added plotting of a second y-axis for Rscript (stable version)

parent 13093ce1
......@@ -6,11 +6,12 @@ text[3 ] <- " Expects the GDX file to contain a 3 dimensional parameter, e.g. p
text[4 ] <- " First index x maps to x-axis, e.g. t1,t2,...t3 (converted to numeric)"
text[5 ] <- " Second index y maps to categories, e.g. oil, gas, nuclear (categories)"
text[6 ] <- " Third index contains parameter values, e.g. 9 ,7 ,5 (numeric expected)"
text[7 ] <- " Example: Rscript do_r_plot.r plot4r.gdx -p p -w 50 -a TRUE"
text[8 ] <- ""
text[9 ] <- " version 2021-12-31"
text[10 ] <- ""
text[11] <- " author Toni Lastusilta (VTT) "
text[7 ] <- " The second parameter p2(x,z) maps to a seconday y-axis (optional)"
text[8 ] <- " Example: Rscript do_r_plot.r plot4r.gdx -p p1 -s p2 -w 50 -a TRUE"
text[9 ] <- ""
text[10 ] <-" version 2022-01-10"
text[11] <- ""
text[12] <- " author Toni Lastusilta (VTT) "
program_description <- vector("character")
for (i in text) {
program_description <- invisible(paste(program_description,i, sep="\n" ))
......@@ -25,6 +26,7 @@ library(rlang)
library(forcats)
library(reshape)
suppressMessages(library(rapportools))
suppressMessages(library(plyr))
#CMD Line Option Management and Error Checking
......@@ -33,12 +35,16 @@ option_list = list(
help="Input file name (GDX) [default= %default]", metavar="character"),
make_option(c("-o", "--output"), type="character", default="Rplots.pdf",
help="Output file name (PDF) [default= %default]", metavar="character"),
make_option(c("-p", "--parameter"), type="character", default="p_r_gen",
make_option(c("-p", "--parameter"), type="character", default="p1",
help="parameter name in gdx, expects 3-dimensional, e.g. p(x,y) [default= %default]", metavar="character"),
make_option(c("-s", "--parameter2"), type="character", default="p2",
help="second parameter name in gdx, expects 3-dimensional, e.g. p2(x,y) [default= %default]", metavar="character"),
make_option(c("-x", "--xlabel"), type="character", default="",
help="Optional X-axis label for plot [default= %default]", metavar="character"),
make_option(c("-y", "--ylabel"), type="character", default="",
help="Optional Y-axis label for plot [default= %default]", metavar="character"),
help="Optional left Y-axis label for plot [default= %default]", metavar="character"),
make_option(c("-z", "--zlabel"), type="character", default="",
help="Optional right Y-axis label for plot [default= %default]", metavar="character"),
make_option(c("-w", "--window"), type="integer", default="",
help="do additional plots by splitting x (x-axis) to the specified time step length [default= %default]", metavar="integer"),
make_option(c("-a", "--archive"), type="logical", default=FALSE,
......@@ -59,7 +65,6 @@ tryCatch({
, error = function(ex) { print(ex) ; stop("FAILED: see above message", call.=FALSE) }
)
# check that gdx file and gdx reader is found
if (is.null(opt$input) || isFALSE(file.exists(opt$input))) {
print_help(opt_parser)
......@@ -69,6 +74,7 @@ if (! require(gdxrrw)) stop ("gdxrrw package is not available")
if (0 == igdx(silent=TRUE)) stop ("the gdx shared library has not been loaded")
# read gdx parameter p of p(x,y)
df_p <- NULL
tryCatch({
gdxName <- opt$input
symName <- opt$parameter
......@@ -76,12 +82,31 @@ tryCatch({
}
, error = function(ex) { print(ex) ; stop("FAILED: see above message", call.=FALSE) }
)
# gdx error check, expecting parameter format p(x,y)
if(ncol(df_p)!=3){
print(head(df_p))
stop(sprintf("Symbol %s has %d columns, expected 3 (third column holds values)", symName, ncol(df_p)), call.=FALSE)
}
df_p_org <- df_p
df_p2 <- NULL
if(!is.empty(opt$parameter2)){
# read gdx parameter p2 of p2(x,y)
tryCatch({
gdxName <- opt$input
symName <- opt$parameter2
df_p2 <- rgdx.param(gdxName,symName, ts=TRUE, squeeze=FALSE)
}
, error = function(ex) { print(ex) ; stop("FAILED: see above message", call.=FALSE) }
)
# gdx error check, expecting parameter format p(x,y)
if(ncol(df_p2)!=3){
print(head(df_p2))
warning(sprintf("Symbol %s ignored, it has %d columns, expected 3 (third column holds values)", symName, ncol(df_p2)), call.=FALSE)
df_p2 <- NULL
}
}
df_p2_org <- df_p2
# read set x and y from parameter p(x,y) if not set by swithches -x and -y
xvar_desc <- vector("character")
......@@ -113,6 +138,22 @@ if(is.empty(opt$ylabel)){
yvar_desc = opt$ylabel
}
# read set z from parameter p2(x,z) if not set by switch -z
zvar_desc <- vector("character")
if(is.empty(opt$zlabel)){
tryCatch({
gdxName <- opt$input
symName <- opt$parameter2
zvar=as_string(sym(colnames(df_p2)[2]))
df_z <- rgdx.set(gdxName, zvar, names=zvar, ts=TRUE)
zvar_desc = paste(attr(df_z,"ts"), " ")
}
, error = function(ex) { zvar_desc="" ; print(ex) ; print("FAILED to read set description: see above message", call.=FALSE) }
)
}else{
zvar_desc = opt$zlabel
}
# store data from gdx
xvar_labels <- unique(df_p[1])
window_end <- nrow(xvar_labels)
......@@ -140,22 +181,85 @@ df_p[,1] <- ordered(df_p[,1])
levels(df_p[,1]) <- 1:nrow(unique(df_p[1]))
df_p[,1] <- as.numeric(as.character(df_p[,1]))
# Create Plots
# Prepare plot p(x,y)
xvar <- sym(colnames(df_p)[1])
yvar <- sym(colnames(df_p)[2])
value <- sym(colnames(df_p)[3])
xlabel <- paste(xvar_desc,xvar,"(",xvar_labels[1,1],"=1, ...,",xvar_labels[window_end,1],"=",window_end,")")
ylabel <- paste(yvar_desc,yvar)
nr_colors <- nrow(unique(df_p[2]))
mypalette="Set3"
if(nr_colors<=8){
# Prepare plot p2(x,z)
nr_shapes <- 0
if(!empty(df_p2)){
xvar2 <- sym(colnames(df_p2)[1])
zvar <- sym(colnames(df_p2)[2])
value2 <- sym(colnames(df_p2)[3])
zlabel <- paste(zvar_desc,zvar)
if(!identical(levels(df_p_org[,1]),levels(df_p2_org[,1]))){
stop(paste("Error: index missmatch. Symbols ",value,"and",value2, " must have same indices for factor",xvar))
}
nr_shapes <- nrow(unique(df_p2[2]))
# convert x-axis to numeric
df_p2[,1] <- ordered(df_p2[,1])
levels(df_p2[,1]) <- 1:nrow(unique(df_p2[1]))
df_p2[,1] <- as.numeric(as.character(df_p2[,1]))
}
#Prepare plot p1 and p2
my_title_font_size = 8
my_axis_title_font_size = 7
my_axis_scale_font_size = 6
my_legend_title_font_size <- 7
my_legend_entry_font_size <- 6
nr_cat <- nr_colors + nr_shapes
if(nr_cat<=8){
mypalette="Dark2"
} else if(nr_cat<=12){
mypalette="Set"
my_legend_title_font_size <- 7
my_legend_entry_font_size <- 6
}else if(nr_cat<=20) {
mypalette=NULL
my_legend_title_font_size <- 4
my_legend_entry_font_size <- 3
}else{
mypalette=NULL
my_legend_title_font_size <- 1
my_legend_entry_font_size <- 1
}
# plot p(x,y)
df_p_index2_ordered <- fct_reorder2(df_p[,2],df_p[,1],df_p[,3])
g1 <- ggplot() +
geom_line(data=df_p, aes(x=!!xvar, y=!!value, group=df_p_index2_ordered, color=df_p_index2_ordered)) +
labs(colour = value) + ggtitle(ptitle) + xlab(xlabel) +
theme(plot.title = element_text(size = my_title_font_size), axis.title=element_text(size=my_axis_title_font_size),
axis.text = element_text(size = my_axis_scale_font_size), legend.position="bottom",
legend.title = element_text(size=my_legend_title_font_size), legend.text = element_text(size=my_legend_entry_font_size))
if(!is.empty(mypalette)){
g1 <- g1 +
scale_color_brewer(palette = mypalette)
}
g2 <- g1
# plot p2(x,z)
if(!empty(df_p2)){
scaleFactor <- max(df_p[,3]) / max(df_p2[,3])
df_p2_index2_ordered <- fct_reorder2(df_p2[,2],df_p2[,1],df_p2[,3])
g2 <- g1 +
geom_line(data=df_p2, aes(x=!!xvar2, y=!!value2*scaleFactor, group=df_p2_index2_ordered, linetype=df_p2_index2_ordered)) +
scale_y_continuous(
name = paste(value, "-axis. ", ylabel),
sec.axis = sec_axis(~./scaleFactor, name=paste(value2, "-axis. ",zlabel))) +
scale_linetype_discrete(name = paste(value2))
}
g <- ggplot(df_p, aes(x=!!xvar, y=!!value, group=!!yvar, color=fct_reorder2(df_p[,2],df_p[,1],df_p[,3]))) + labs(colour = yvar) + geom_line() + theme(text = element_text(size=6))+ scale_color_brewer(palette = mypalette) + ggtitle(ptitle) + xlab(xlabel) + ylab(ylabel) + theme(legend.position="bottom", legend.text = element_text(size=6))
#g2_copy<-g2 # to avoid errors we create a copy, g2 changes with switch -w
#plot(g2_copy)
# create PDF
pdf(file = pdf_name, width = 6.25, height = 4, family = "Times", pointsize = 6, onefile = TRUE)
print(g)
print(g2)
cnt_plots=1;
# additional figures to PDF with switch -w (window)
......@@ -166,7 +270,7 @@ if(window_step>=1){
window_curr_end=window_end
}
xlabel <- paste(xvar_desc,xvar,"(",xvar_labels[window_curr_start,1],"=",window_curr_start,", ...,",xvar_labels[window_curr_end,1],"=",window_curr_end,")")
print(g + xlab(xlabel) + coord_cartesian(xlim = c(window_curr_start,window_curr_end), expand=0))
print(g2 + xlab(xlabel) + coord_cartesian(xlim = c(window_curr_start,window_curr_end), expand=0))
cnt_plots=cnt_plots+1
}
}
......
......@@ -15,15 +15,20 @@ $batinclude backbone_definitions_4plots.gms
* Script to call R for making plots
$onechoV > runR.inc
$$log runR.inc creates line graphs with R from a two-dimensional parameter in a GDX file
$$log The file uses compile time variables %Ylabel% and %Xlabel% in case they are set ($set)
$$log Usage: $batInclude runR.gms GDXfile GDXparam
$$log Example: $batInclude runR.gms plot4r.gdx p value time
$$log Now: $batInclude runR.gms %1 %2 %3 %4
$$log It can also plot a second parameter
$$log The file uses compile time variables %Ylabel%, %Ylabel2% and %Xlabel% in case they are set ($set)
$$log Usage: $batInclude runR.gms GDXfile GDXparam GDXparam2
$$log Example: $batInclude runR.gms plot4r.gdx p p2
$$log Now: $batInclude runR.gms %1 %2 %3
$$set r_param2
$$if not "%3"=="" $set r_param2 -s %3
$$set r_ylabel
$$if set Ylabel $set r_ylabel -y "%Ylabel%"
$$set r_ylabel2
$$if set Ylabel2 $set r_ylabel2 -z "%Ylabel2%"
$$set r_xlabel
$$if set Xlabel $set r_xlabel -x "%Xlabel%"
execute 'Rscript "%gams.wDir%plots\do_r_plot.r" -i "%gams.wDir%plots/%1" -o "%gams.wDir%output/%2.pdf" -p %2 -w 720 -a FALSE %r_ylabel% %r_xlabel% ';
execute 'Rscript "%gams.wDir%plots\do_r_plot.r" -i "%gams.wDir%plots/%1" -o "%gams.wDir%output/%2.pdf" -p %2 %r_param2% -w 168 -a FALSE %r_ylabel% %r_ylabel2% %r_xlabel% ';
myerrorlevel = errorlevel;
if(myerrorlevel=0,
Execute.ASyncNC 'SumatraPDF.exe "%gams.wDir%output/%2.pdf"';
......@@ -32,22 +37,35 @@ $onechoV > runR.inc
)
$offecho
* Define source GDX for plotting data
$set backbone_output_GDX "output\archive\out_alt_obj.gdx"
* Preparation
$set GDXparam2 priceElec
$set GDXparam2set2 price_2013
* Define sets and parameters to be used to prepare data for plots or in ggplots in R
Set
tparam(t) "Temporary: Time steps used for a parameter"
tparam(t) "Time steps (in hours)"
%GDXparam2set2% "Electricity price (EUR/MWh)" /price/
;
Parameter
myerrorlevel Stop execution if error in external call / 0 /
w2kw Convert whatt to kilowhatt / 0.001 /
sign Switch sign / -1 /
myerrorlevel "Save Rscript callback error level" / 0 /
w2kw "Convert whatt to kilowhatt" / 0.001 /
sign "Switch sign" / -1 /
%GDXparam2%(t,%GDXparam2set2%) "Electricity spot price (EUR/MWh)"
;
* Define source GDX for plotting data
$set backbone_output_GDX "output\archive\out_alt_obj.gdx"
* Prepare priceElspot (continued)
execute_loaddc "%backbone_output_GDX%", tparam=t_realized, ts_priceElspot;
%GDXparam2%(tparam,%GDXparam2set2%)=ts_priceElspot(tparam) +eps ;
$set Ylabel2 "Price of electricity (EUR/MWh) :"
* CREATE PLOTS:
* %GDXfile% for R, %GDXparam% parameter to plot, %TMPParam% two-dimensional counterpart to %GDXparam%
* Common for all plots
* None
* %GDXfile% for R, %GDXparam% parameter to plot, %TMPParam% two-dimensional counterpart to %GDXparam%
* Plots 1 : ts_influx
* The influence of external power (weather) on building, namely envelope_mass due outside temperature changes and interior air due to ventilation
* Plot ts_influx_absolute *******************************************************
......@@ -59,8 +77,8 @@ Parameter %TMPparam%(t,node) "External power (weather) influencing building" ;
execute_loaddc "%backbone_output_GDX%", tparam=t_realized, %GDXparam%;
* Below line is adjusted for each paramGDX and paramPlot pair
%TMPparam%(tparam(t),node)=sum((grid,f)$(sameas(grid,"building") and sameas(f,'f00')),w2kw*%GDXparam%(grid,node,f,t) + eps);
execute_unload 'plots/%GDXfile%', tparam=t, node, %TMPparam%;
$batInclude runR.inc %GDXfile% %TMPparam%
execute_unload 'plots/%GDXfile%', tparam=t, node, %TMPparam%, %GDXparam2set2%, %GDXparam2% ;
$batInclude runR.inc %GDXfile% %TMPparam% %GDXparam2%
* Plot ts_influx_normalized *******************************************************
$set GDXfile plot4r.gdx
$set GDXparam ts_influx
......@@ -70,8 +88,8 @@ Parameter %TMPparam%(t,node) "External power (weather) influencing building -
execute_loaddc "%backbone_output_GDX%", tparam=t_realized, %GDXparam%;
* Below line is adjusted for each paramGDX and paramPlot pair
%TMPparam%(tparam(t),node)=sum((grid,node_building2node(node_building,node),f)$(sameas(grid,"building") and sameas(f,'f00')), w2kw*%GDXparam%(grid,node,f,t)/building_squares(node_building) + eps);
execute_unload 'plots/%GDXfile%', tparam=t, node, %TMPparam%;
$batInclude runR.inc %GDXfile% %TMPparam%
execute_unload 'plots/%GDXfile%', tparam=t, node, %TMPparam%, %GDXparam2set2%, %GDXparam2% ;
$batInclude runR.inc %GDXfile% %TMPparam% %GDXparam2%
* Plots 2 : r_gen
* Electric power consumption for heating and cooling of air and heating of domestic hot water
......@@ -84,8 +102,8 @@ Parameter %TMPparam%(t,unit) "Electric power consumption for heating and cooli
execute_loaddc "%backbone_output_GDX%", tparam=t_realized, %GDXparam%;
* Below line is adjusted for each paramGDX and paramPlot pair
%TMPparam%(tparam(t),unit_heat_and_cool(unit))=sum((grid,node,f)$(sameas(grid,"elec") and sameas(node,'link_node') and sameas(f,'f00')),sign* w2kw* %GDXparam%(grid,node,unit,f,t) + eps);
execute_unload 'plots/%GDXfile%', tparam=t, unit, %TMPparam%;
$batInclude runR.inc %GDXfile% %TMPparam%
execute_unload 'plots/%GDXfile%', tparam=t, unit, %TMPparam%, %GDXparam2set2%, %GDXparam2% ;
$batInclude runR.inc %GDXfile% %TMPparam% %GDXparam2%
* Plot r_gen_heat_normalized *******************************************************
$set GDXfile plot4r.gdx
$set GDXparam r_gen
......@@ -95,8 +113,8 @@ Parameter %TMPparam%(t,unit) "Electric power consumption for heating and cooli
execute_loaddc "%backbone_output_GDX%", tparam=t_realized, %GDXparam%;
* Below line is adjusted for each paramGDX and paramPlot pair
%TMPparam%(tparam(t),unit_heat_and_cool(unit))=sum((grid,node,node_building2unit(node_building,unit),f)$(sameas(grid,"elec") and sameas(node,'link_node') and sameas(f,'f00')),sign* w2kw* %GDXparam%(grid,node,unit,f,t)/building_squares(node_building) + eps);
execute_unload 'plots/%GDXfile%', tparam=t, unit, %TMPparam%;
$batInclude runR.inc %GDXfile% %TMPparam%
execute_unload 'plots/%GDXfile%', tparam=t, unit, %TMPparam%, %GDXparam2set2%, %GDXparam2% ;
$batInclude runR.inc %GDXfile% %TMPparam% %GDXparam2%
* Plot r_gen DHW absolute *******************************************************
$set GDXfile plot4r.gdx
$set GDXparam r_gen
......@@ -106,8 +124,8 @@ Parameter %TMPparam%(t,unit) "Electric power consumption for domestic hot wate
execute_loaddc "%backbone_output_GDX%", tparam=t_realized, %GDXparam%;
* Below line is adjusted for each paramGDX and paramPlot pair
%TMPparam%(tparam(t),unit_DHW(unit))=sum((grid,node,f)$(sameas(grid,"elec") and sameas(node,'link_node') and sameas(f,'f00')),sign* w2kw* %GDXparam%(grid,node,unit,f,t) + eps);
execute_unload 'plots/%GDXfile%', tparam=t, unit, %TMPparam%;
$batInclude runR.inc %GDXfile% %TMPparam%
execute_unload 'plots/%GDXfile%', tparam=t, unit, %TMPparam%, %GDXparam2set2%, %GDXparam2% ;
$batInclude runR.inc %GDXfile% %TMPparam% %GDXparam2%
* Plot r_gen DHW normalized *******************************************************
$set GDXfile plot4r.gdx
$set GDXparam r_gen
......@@ -117,8 +135,8 @@ Parameter %TMPparam%(t,unit) "Electric power consumption for domestic hot wate
execute_loaddc "%backbone_output_GDX%", tparam=t_realized, %GDXparam%;
* Below line is adjusted for each paramGDX and paramPlot pair
%TMPparam%(tparam(t),unit_DHW(unit))=sum((grid,node,node_building2unit(node_building,unit),f)$(sameas(grid,"elec") and sameas(node,'link_node') and sameas(f,'f00')),sign* w2kw* %GDXparam%(grid,node,unit,f,t)/building_squares(node_building) + eps);
execute_unload 'plots/%GDXfile%', tparam=t, unit, %TMPparam%;
$batInclude runR.inc %GDXfile% %TMPparam%
execute_unload 'plots/%GDXfile%', tparam=t, unit, %TMPparam%, %GDXparam2set2%, %GDXparam2% ;
$batInclude runR.inc %GDXfile% %TMPparam% %GDXparam2%
* Plots 3 r_state
* Current building temperature for interior_air_and_furniture, internal_mass, envelope_mass and DHWT
......@@ -131,8 +149,8 @@ Parameter %TMPparam%(t,node) "Temperature at interior_air_and_furniture (Celci
execute_loaddc "%backbone_output_GDX%", tparam=t_realized, %GDXparam%;
* Below line is adjusted for each paramGDX and paramPlot pair
%TMPparam%(tparam(t),node_building_interior_air_and_furniture(node))=sum((grid,f)$(sameas(grid,"building") and sameas(f,'f00')),%GDXparam%(grid,node,f,t) - c2k + eps);
execute_unload 'plots/%GDXfile%', tparam=t, node, %TMPparam%;
$batInclude runR.inc %GDXfile% %TMPparam%
execute_unload 'plots/%GDXfile%', tparam=t, node, %TMPparam%, %GDXparam2set2%, %GDXparam2% ;
$batInclude runR.inc %GDXfile% %TMPparam% %GDXparam2%
* Plot r_state(internal_mass) *******************************************************
$set GDXfile plot4r.gdx
$set GDXparam r_state
......@@ -142,8 +160,8 @@ Parameter %TMPparam%(t,node) "Temperature at internal_mass: inside walls (Cel
execute_loaddc "%backbone_output_GDX%", tparam=t_realized, %GDXparam%;
* Below line is adjusted for each paramGDX and paramPlot pair
%TMPparam%(tparam(t),node_building_internal_mass(node))=sum((grid,f)$(sameas(grid,"building") and sameas(f,'f00')),%GDXparam%(grid,node,f,t) - c2k + eps);
execute_unload 'plots/%GDXfile%', tparam=t, node, %TMPparam%;
$batInclude runR.inc %GDXfile% %TMPparam%
execute_unload 'plots/%GDXfile%', tparam=t, node, %TMPparam%, %GDXparam2set2%, %GDXparam2% ;
$batInclude runR.inc %GDXfile% %TMPparam% %GDXparam2%
* Plot r_state(envelope_mass) *******************************************************
$set GDXfile plot4r.gdx
$set GDXparam r_state
......@@ -153,8 +171,8 @@ Parameter %TMPparam%(t,node) "Temperature at envelope_mass: outside walls (Cel
execute_loaddc "%backbone_output_GDX%", tparam=t_realized, %GDXparam%;
* Below line is adjusted for each paramGDX and paramPlot pair
%TMPparam%(tparam(t),node_building_envelope_mass(node))=sum((grid,f)$(sameas(grid,"building") and sameas(f,'f00')),%GDXparam%(grid,node,f,t) - c2k + eps);
execute_unload 'plots/%GDXfile%', tparam=t, node, %TMPparam%;
$batInclude runR.inc %GDXfile% %TMPparam%
execute_unload 'plots/%GDXfile%', tparam=t, node, %TMPparam%, %GDXparam2set2%, %GDXparam2% ;
$batInclude runR.inc %GDXfile% %TMPparam% %GDXparam2%
* Plot r_state(DHWT) *******************************************************
$set GDXfile plot4r.gdx
$set GDXparam r_state
......@@ -164,7 +182,5 @@ Parameter %TMPparam%(t,node) "Temperature at Domestic Hot Water Tank (Celcius)
execute_loaddc "%backbone_output_GDX%", tparam=t_realized, %GDXparam%;
* Below line is adjusted for each paramGDX and paramPlot pair
%TMPparam%(tparam(t),node_building_DHWT(node))=sum((grid,f)$(sameas(grid,"building") and sameas(f,'f00')),%GDXparam%(grid,node,f,t) - c2k + eps);
execute_unload 'plots/%GDXfile%', tparam=t, node, %TMPparam%;
$batInclude runR.inc %GDXfile% %TMPparam%
execute_unload 'plots/%GDXfile%', tparam=t, node, %TMPparam%, %GDXparam2set2%, %GDXparam2% ;
$batInclude runR.inc %GDXfile% %TMPparam% %GDXparam2%
*Create data for simple Rscript create_simple_r_plot.r
$title Create GDX files for testing RScript plotting
$ontext
Create data for simple Rscript "create_simple_r_plot.r"
Addiotnally we create data for "do_r_plot.r"
$offtext
set
i "Set i " /i1*i5/
j "Set j " /j1*j3/
k "Set k " /k1*k2/
;
parameter
p(i,j) "parameter p"
p1(i,j) "parameter to plot on first y-axis"
p2(i,k) "parameter to plot on second y-axis"
;
p(i,j)=ord(i)*ord(j);
p1(i,j)=ord(i)*ord(j);
p2(i,k)=(1+card(i)-ord(i))*10*ord(k);
execute_unload "%gams.wDir%\plots\simple_r_plot_data.gdx";
*execute 'Rscript "%gams.curDir%plots\create_simple_r_plot.r"';
* Create data for do_r_plot.r
execute_unload "%gams.wDir%\plots\plot4r.gdx";
# A simple example to plot a GDX file with R
library(reshape2)
library(plyr)
library(gdxrrw)
library(ggplot2)
r_dir = "/Users/tltoni/GIT/FlexiB/backbone_building_v2/plots"
# we use gams backbone project directory to avoid path issues
gams_project_dir = "/Users/tltoni/GIT/FlexiB/backbone_building_v2"
gdx_filename = "simple_r_plot_data.gdx"
setwd(r_dir)
if (! require(gdxrrw)) stop ("gdxrrw package is not available")
if (0 == igdx(silent=TRUE)) stop ("the gdx shared library has not been loaded")
if(!file.exists(gdx_filename)){
exec_string= sprintf("gams plots/simple_r_plot_create_data.gms wdir=%s",gams_project_dir)
system(exec_string)
}
#gdxInfo("gdx_filename")
tryCatch({
gdxName <- gdx_filename
symName <- "p1"
df1 <- rgdx.param(gdxName,symName)
# ix <- rgdx.set(gdxName,'i',names="I")
# jx <- rgdx.set(gdxName,'j',names="J")
}
, error = function(ex) { stop(ex) ; FALSE }
)
tryCatch({
gdxName <- gdx_filename
symName2 <- "p2"
df2 <- rgdx.param(gdxName,symName2)
# ix <- rgdx.set(gdxName,'i',names="I")
# kx <- rgdx.set(gdxName,'k',names="k")
}
, error = function(ex) { stop(ex) ; FALSE }
)
# df3 <- merge(df1,df2,by="i")
g1 <- ggplot() +
geom_line(data=df1, aes(x=i, y=p1, group=j, color=j)) +
scale_color_discrete(name = paste(symName)) +
guides(color=guide_legend(override.aes=list(fill=NA)))
g2 <- g1
#scaleFactor <- 1
scaleFactor <- max(df1$p1) / max(df2$p2)
if(identical(levels(df1[,1]),levels(df2[,1]))){
g2 <- g1 +
geom_line(data=df2, aes(x=i, y=p2*scaleFactor, group=k, linetype=k )) +
scale_y_continuous(
name = paste(symName, "-axis"),
sec.axis = sec_axis(~./scaleFactor, name=paste(symName2, "-axis"))) +
scale_linetype_discrete(name=paste(symName2))
}else{
print(paste("Warning! Parameter ", symName2, " is ignored"))
print(paste(symName, " differs from ", symName2," in set ", sym(colnames(df1)[1])))
}
plot(g2)
......@@ -21,14 +21,14 @@ if(!file.exists(gdx_filename)){
#gdxInfo("gdx_filename")
tryCatch({
gdxName <- gdx_filename
symName <- "p"
df <- rgdx.param(gdxName,symName)
symName <- "p1"
df1 <- rgdx.param(gdxName,symName)
# ix <- rgdx.set(gdxName,'i',names="I")
# jx <- rgdx.set(gdxName,'j',names="J")
}
, error = function(ex) { stop(ex) ; FALSE }
)
g <- ggplot(df, aes(x=i, y=p, group=j, color=j)) + geom_line()
plot(g)
g1 <- ggplot(df1, aes(x=i, y=p1, group=j, color=j)) + geom_line()
plot(g1)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment