#' Produce polygon plots of time series for aggregated loglike indicators
#'
#' This function produces all accumulated time series over month to compare scenario outcomes
#' All the plots are being stored in a polygons_plots folder that can further be found in the output folder.
#' (compare up to 5 scenarios simultaneously)
#' @param fname First name
#' @param lname Last name
#' @export
#' @examples
#' \dontrun{
#' general <- setGeneralOverallVariable (pathToRawInputs =file.path("C:", "Users", "fbas",
#' "Documents", "GitHub", paste0("DISPLACE_input_gis_",
#' "DanishFleet")),
#' pathToDisplaceInputs = file.path("C:", "Users", "fbas",
#' "Documents", "GitHub", paste0("DISPLACE_input_", "DanishFleet")),
#' pathToOutputs =file.path("C:","DISPLACE_outputs"),
#' caseStudy="DanishFleet",
#' iGraph=41,
#' iYear="2015",
#' iCountry="DEN",
#' nbPops=39,
#' nbSzgroup=14,
#' theScenarios= c("svana_baseline",
#' "svana_sub1mx20",
#' "svana_sub4mx20",
#' "svana_sub4mx5ns20bt",
#' "svana_sub4mx20ns5bt",
#' "svana_sub4mx5ns5bt" ),
#' nbSimus=20,
#' useSQLite=FALSE
#' )
#'
#'
#'
#'
#' loadLoglikeFiles(general=general, use_port_info=FALSE)
#'
#'
#'
#' polygonPlotsFromAggLoglikeFiles (general=general,
#' the_baseline="svana_baseline",
#' a_width=3500,
#' a_height=1000,
#' selected_scenarios=general$namefolderoutput[1:3],
#' the_scenario_names=general$namefolderoutput[1:3],
#' nby=5
#' )
#'
#' }
polygonPlotsFromAggLoglikeFiles <- function(general=general,
the_baseline="svana_baseline",
a_width=3500,
a_height=1000,
selected_scenarios=general$namefolderoutput[1:3],
the_scenario_names=general$namefolderoutput[1:3],
nby=5
)
{
sce1 <- selected_scenarios[1] # init
do_polygon_plot <- function(
a_variable="gradva",
nby=5,
a_set_of_scenarios=general$namefolderoutput[1:3],
the_scenario_names= general$namefolderoutput[1:3],
name_set_of_sces= "setA",
selected=selected,
export=TRUE,
a_xlab="# months",
a_ylab="Accumulated Gross Added Value (millions Euro)",
add_legend=FALSE,
color_legend= c(rgb(94/255,79/255,162/255,0.5), rgb (158/255,1/255,66/255,0.5),
rgb(140/255,81/255,10/255,0.4), rgb(1,0,0,0.5), rgb(0,0.5,1,0.5), rgb(0,1,0.5,0.5), rgb(1,0,0.5,0.5), rgb(0,0,0,0.2)) ,
a_width=a_width,
a_height=a_height
) {
documsum <- TRUE
if(a_variable %in% c("GVAPerFTE")) documsum <- FALSE
sce1 <- a_set_of_scenarios[1]
sce2 <- a_set_of_scenarios[2]
sce3 <- a_set_of_scenarios[3]
if(length(a_set_of_scenarios)>=4) sce4 <- a_set_of_scenarios[4]
if(length(a_set_of_scenarios)>=5) sce5 <- a_set_of_scenarios[5]
obj1 <- get(paste("lst_loglike_agg_weight", selected, sce1, sep=''), env=.GlobalEnv)
obj2 <- get(paste("lst_loglike_agg_weight", selected, sce2, sep=''), env=.GlobalEnv)
obj3 <- get(paste("lst_loglike_agg_weight", selected, sce3, sep=''), env=.GlobalEnv)
if(length(a_set_of_scenarios)>=4) obj4 <- get(paste("lst_loglike_agg_weight", selected, sce4, sep=''), env=.GlobalEnv)
if(length(a_set_of_scenarios)>=5) obj5 <- get(paste("lst_loglike_agg_weight", selected, sce5, sep=''), env=.GlobalEnv)
# caution: complete missing records if 0s in given year.month
complete_all_year_month <- function (x, years_span=2015:(2015+nby-1)){
allcombi <- expand.grid(month=sprintf("%02d", 1:12), year=years_span)
allcombi$year.month <- paste0(allcombi$year,".",allcombi$month)
allcombi <- cbind.data.frame(year.month=allcombi$year.month, matrix(0, ncol=ncol(x)-1))
colnames(allcombi) <- colnames(x)
allmissingcombi <- allcombi[!allcombi$year.month %in% x$year.month,]
dd <- rbind.data.frame(x, allmissingcombi)
rownames(dd) <- dd$year.month
dd <- dd[as.character(allcombi$year.month),] # get the right order...
return(dd)
}
obj1 <- lapply(obj1, complete_all_year_month)
obj2 <- lapply(obj2, complete_all_year_month)
obj3 <- lapply(obj3, complete_all_year_month)
if(length(a_set_of_scenarios)>=4) obj4 <- lapply(obj4, complete_all_year_month)
if(length(a_set_of_scenarios)>=5) obj5 <- lapply(obj5, complete_all_year_month)
simu_names <- names(obj1)
simu_names <- simu_names[simu_names %in% names(obj2)]
simu_names <- simu_names[simu_names %in% names(obj3)]
if(length(a_set_of_scenarios)>=4) simu_names <- simu_names[simu_names %in% names(obj4)]
if(length(a_set_of_scenarios)>=5) simu_names <- simu_names[simu_names %in% names(obj5)]
## a nice quantile plot for profit
#plot(cumsum(loglike_Scenario1_save[loglike_Scenario1_save$simu=="simu2",]$gradva), type="l", col=2)
mat_sce_base <- matrix(NA, nrow=length(simu_names), ncol=nby*12) # assuming max 11000 records
rownames(mat_sce_base) <- simu_names
mat_sce1 <- matrix(NA, nrow=length(simu_names), ncol=nby*12)
rownames(mat_sce1) <- simu_names
mat_sce2 <- matrix(NA, nrow=length(simu_names), ncol=nby*12)
rownames(mat_sce2) <- simu_names
if(length(a_set_of_scenarios)>=4) mat_sce3 <- matrix(NA, nrow=length(simu_names), ncol=nby*12)
if(length(a_set_of_scenarios)>=4) rownames(mat_sce3) <- simu_names
if(length(a_set_of_scenarios)>=5) mat_sce4 <- matrix(NA, nrow=length(simu_names), ncol=nby*12)
if(length(a_set_of_scenarios)>=5) rownames(mat_sce4) <- simu_names
for (sim in simu_names){
if(documsum){
mat_sce_base[sim, ] <-
cumsum(as.numeric(obj1[[sim]][,a_variable]) )[1:dim(mat_sce_base)[2]]
mat_sce1[sim, ] <-
cumsum(as.numeric(obj2[[sim]][,a_variable]) )[1:dim(mat_sce1)[2]]
mat_sce2[sim, ] <-
cumsum(as.numeric(obj3[[sim]][,a_variable]) )[1:dim(mat_sce2)[2]]
if(length(a_set_of_scenarios)==4) mat_sce3[sim, ] <-
cumsum(as.numeric(obj4[[sim]][,a_variable]) )[1:dim(mat_sce3)[2]]
if(length(a_set_of_scenarios)>=5) mat_sce4[sim, ] <-
cumsum(as.numeric(obj5[[sim]][,a_variable]) )[1:dim(mat_sce4)[2]]
} else{
mat_sce_base[sim, ] <-
as.numeric(obj1[[sim]][,a_variable]) [1:dim(mat_sce_base)[2]]
mat_sce1[sim, ] <-
as.numeric(obj2[[sim]][,a_variable]) [1:dim(mat_sce1)[2]]
mat_sce2[sim, ] <-
as.numeric(obj3[[sim]][,a_variable])[1:dim(mat_sce2)[2]]
if(length(a_set_of_scenarios)==4) mat_sce3[sim, ] <-
as.numeric(obj4[[sim]][,a_variable]) [1:dim(mat_sce3)[2]]
if(length(a_set_of_scenarios)>=5) mat_sce4[sim, ] <-
as.numeric(obj5[[sim]][,a_variable]) [1:dim(mat_sce4)[2]]
}
}
if(export) tiff(file=file.path(general$main.path, general$namefolderoutput,
paste("accumulated",a_variable,"perScenario", selected,"_",name_set_of_sces, ".tiff", sep="")),
width=a_width, height=a_height,
compression="lzw")
sim_ref <- names(which.max (apply(mat_sce_base, 1, function(x) sum(as.numeric(x), na.rm=TRUE))) )
er <- try( {
if(documsum) plot(cumsum(as.numeric(obj1[[sim_ref]][,a_variable]))[1:(nby*12)] /1e6, ylim=c(min(cumsum(as.numeric(obj1[[sim_ref]][,a_variable]))[1:(nby*12)] /1e6), max(cumsum(as.numeric(obj1[[sim_ref]][,a_variable]))[1:(nby*12)] /1e6)),
col=2, type="n", xlab="Nb months", ylab="", cex.lab=1.6, axes=FALSE)
if(!documsum) plot(as.numeric(obj1[[sim_ref]][,a_variable])[1:(nby*12)] /1e6, ylim=c(min(as.numeric(obj1[[sim_ref]][,a_variable])[1:(nby*12)] /1e6), max(as.numeric(obj1[[sim_ref]][,a_variable])[1:(nby*12)] /1e6)),
col=2, type="n", xlab="Nb months", ylab="", cex.lab=1.6, axes=FALSE)
polygon(x=c((1:ncol(mat_sce_base))[1:(nby*12)], rev((1:ncol(mat_sce_base))[1:(nby*12)])),
y=c(apply(mat_sce_base[,1:(nby*12)], 2, quantile, probs=c(0.05, 0.95), na.rm=TRUE)["5%",], rev(apply(mat_sce_base[,1:(nby*12)], 2, quantile, probs=c(0.05, 0.95), na.rm=TRUE)["95%",])) /1e6,
col= color_legend[1], border=NA) # blue
polygon(x=c((1:ncol(mat_sce1))[1:(nby*12)], rev((1:ncol(mat_sce1))[1:(nby*12)])),
y=c(apply(mat_sce1[,1:(nby*12)], 2, quantile, probs=c(0.05, 0.95), na.rm=TRUE)["95%",], rev(apply(mat_sce1[,1:(nby*12)], 2, quantile, probs=c(0.05, 0.95), na.rm=TRUE)["5%",])) /1e6,
col= color_legend[2], border=NA) # green
polygon(x=c((1:ncol(mat_sce2))[1:(nby*12)], rev((1:ncol(mat_sce2))[1:(nby*12)])),
y=c(apply(mat_sce2[,1:(nby*12)], 2, quantile, probs=c(0.05, 0.95), na.rm=TRUE)["5%",], rev(apply(mat_sce2[,1:(nby*12)], 2, quantile, probs=c(0.05, 0.95), na.rm=TRUE)["95%",])) /1e6,
col= color_legend[3] , border=NA) # red
if(length(a_set_of_scenarios)>=4) polygon(x=c((1:ncol(mat_sce3))[1:(nby*12)], rev((1:ncol(mat_sce3))[1:(nby*12)])),
y=c(apply(mat_sce3[,1:(nby*12)], 2, quantile, probs=c(0.05, 0.95), na.rm=TRUE)["5%",], rev(apply(mat_sce3[,1:(nby*12)], 2, quantile, probs=c(0.05, 0.95), na.rm=TRUE)["95%",])) /1e6,
col= color_legend[4], border=NA) # grey
if(length(a_set_of_scenarios)>=5) polygon(x=c((1:ncol(mat_sce4))[1:(nby*12)], rev((1:ncol(mat_sce4))[1:(nby*12)])),
y=c(apply(mat_sce4[,1:(nby*12)], 2, quantile, probs=c(0.05, 0.95), na.rm=TRUE)["5%",], rev(apply(mat_sce4[,1:(nby*12)], 2, quantile, probs=c(0.05, 0.95), na.rm=TRUE)["95%",])) /1e6,
col= color_legend[5], border=NA) # grey
axis(2, las=2, cex.axis=1.5)
axis(1, cex.axis=1.5)
if(add_legend) legend("topleft", fill=color_legend, border =color_legend, legend=the_scenario_names, cex=1.3, bty="n")
box()
abline(h=0, lty=2, col=grey(0.9))
mtext(side = 2, text = a_ylab, line = 3.6, cex=1.2)
}, silent=TRUE)
if(class(er)=="try-error"){
print(paste("no data."))
browser()
}
if(export) dev.off()
return()
}
var_names <- colnames(get(paste("lst_loglike_agg_weight_","selected_set1_",sce1, sep=''), env=.GlobalEnv)[[1]])
dir.create(file.path(general$main.path, general$namefolderinput,"polygon_plots"))
for (a_var in var_names[-1]){
graphics.off()
tiff(file=file.path(general$main.path, general$namefolderinput, "polygon_plots", paste("accumulated_per_scenario_polygon_", a_var, ".tiff", sep="")),
width = a_width, height = a_height, compression="lzw",
units = "px", pointsize = 12, res=300)
par(mfrow=c(1,3))
par(mar=c(2,4.4,2,2))
par(oma=c(4,4,1,1))
cat (paste0("plot for ", a_var,"\n"))
do_polygon_plot(
a_variable=a_var,
nby=nby,
a_set_of_scenarios= selected_scenarios,
the_scenario_names= selected_scenarios,
name_set_of_sces= "setA",
selected="_selected_set1_",
export=FALSE,
a_xlab="# months",
if(a_var=="gradva") {a_ylab="Acc. GVA (mio Euro)"} else{ if(a_var=="rev_explicit_from_av_prices"){a_ylab="Income from landings (mio Euro)"} else{a_ylab=a_var}},
add_legend=TRUE,
color_legend= c(rgb(94/255,79/255,162/255,0.5), rgb (158/255,1/255,66/255,0.5), rgb(140/255,81/255,10/255,0.4), rgb(1,0,0,0.5), rgb(0,0.5,1,0.5), rgb(0,1,0.5,0.5), rgb(1,0,0.5,0.5), rgb(0,0,0,0.2)),
a_width=a_width,
a_height=a_height)
cat (paste0("plot for ", a_var," for set1...ok\n"))
do_polygon_plot(
a_variable=a_var,
nby=nby,
a_set_of_scenarios= selected_scenarios,
the_scenario_names= selected_scenarios,
name_set_of_sces= "setA",
selected="_selected_set2_",
export=FALSE,
a_xlab="# months",
if(a_var=="gradva") {a_ylab="Acc. GVA (mio Euro)"} else{ if(a_var=="rev_explicit_from_av_prices"){a_ylab="Income from landings (mio Euro)"} else{a_ylab=a_var}},
add_legend=TRUE,
color_legend= c(rgb(94/255,79/255,162/255,0.5), rgb (158/255,1/255,66/255,0.5), rgb(140/255,81/255,10/255,0.4), rgb(1,0,0,0.5), rgb(0,0.5,1,0.5), rgb(0,1,0.5,0.5), rgb(1,0,0.5,0.5), rgb(0,0,0,0.2)),
a_width=a_width,
a_height=a_height
)
cat (paste0("plot for ", a_var," for set2...ok\n"))
do_polygon_plot(
a_variable=a_var,
nby=nby,
a_set_of_scenarios= selected_scenarios,
the_scenario_names= selected_scenarios,
name_set_of_sces= "setA",
selected="_selected_set3_",
export=FALSE,
a_xlab="# months",
if(a_var=="gradva") {a_ylab="Acc. GVA (mio Euro)"} else{ if(a_var=="rev_explicit_from_av_prices"){a_ylab="Income from landings (mio Euro)"} else{a_ylab=a_var}},
add_legend=TRUE,
color_legend= c(rgb(94/255,79/255,162/255,0.5), rgb (158/255,1/255,66/255,0.5), rgb(140/255,81/255,10/255,0.4), rgb(1,0,0,0.5), rgb(0,0.5,1,0.5), rgb(0,1,0.5,0.5), rgb(1,0,0.5,0.5), rgb(0,0,0,0.2)),
a_width=a_width,
a_height=a_height
)
cat (paste0("plot for ", a_var," for set3...ok\n"))
mtext("# months", 1, line=2, cex=1.5, outer=TRUE)
mtext(side=2,"Indicators",line=2, cex=1.5, outer=TRUE)
dev.off()
} # end a_var
return()
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.