diff --git a/.buildlibrary b/.buildlibrary index 429fbd24..e319e051 100644 --- a/.buildlibrary +++ b/.buildlibrary @@ -1,4 +1,4 @@ -ValidationKey: '231860058' +ValidationKey: '232170880' AcceptedWarnings: - .*following variables are expected in the piamInterfaces.* - Summation checks have revealed some gaps.* diff --git a/CITATION.cff b/CITATION.cff index f7327740..0e488d5a 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -2,8 +2,8 @@ cff-version: 1.2.0 message: If you use this software, please cite it using the metadata from this file. type: software title: 'remind2: The REMIND R package (2nd generation)' -version: 1.158.2 -date-released: '2024-10-23' +version: 1.159.0 +date-released: '2024-11-05' abstract: Contains the REMIND-specific routines for data and model output manipulation. authors: - family-names: Rodrigues diff --git a/DESCRIPTION b/DESCRIPTION index 4b8b3c05..4568ae9f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Type: Package Package: remind2 Title: The REMIND R package (2nd generation) -Version: 1.158.2 -Date: 2024-10-23 +Version: 1.159.0 +Date: 2024-11-05 Authors@R: c( person("Renato", "Rodrigues", , "renato.rodrigues@pik-potsdam.de", role = c("aut", "cre")), person("Lavinia", "Baumstark", role = "aut"), diff --git a/NAMESPACE b/NAMESPACE index 8dcf8a57..ce4df492 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -136,7 +136,6 @@ importFrom(ggplot2,ggtitle) importFrom(ggplot2,labs) importFrom(ggplot2,scale_color_manual) importFrom(ggplot2,scale_fill_brewer) -importFrom(ggplot2,scale_fill_discrete) importFrom(ggplot2,scale_fill_manual) importFrom(ggplot2,scale_linetype_identity) importFrom(ggplot2,scale_x_continuous) @@ -148,7 +147,6 @@ importFrom(ggplot2,theme) importFrom(ggplot2,theme_bw) importFrom(ggplot2,theme_minimal) importFrom(ggplot2,unit) -importFrom(ggplot2,xlab) importFrom(gms,getLine) importFrom(gms,readDefaultConfig) importFrom(lucode2,getScenNames) @@ -179,6 +177,7 @@ importFrom(magclass,getNames) importFrom(magclass,getRegions) importFrom(magclass,getSets) importFrom(magclass,getYears) +importFrom(magclass,is.magpie) importFrom(magclass,lowpass) importFrom(magclass,magpie_expand) importFrom(magclass,mbind) @@ -222,7 +221,6 @@ importFrom(quitte,overwrite) importFrom(quitte,read.gdx) importFrom(quitte,read.quitte) importFrom(quitte,revalue.levels) -importFrom(quitte,sum_total) importFrom(readr,read_csv) importFrom(readr,write_rds) importFrom(remulator,calc_supplycurve) diff --git a/R/convGDX2MIF.R b/R/convGDX2MIF.R index 50718201..c675813b 100644 --- a/R/convGDX2MIF.R +++ b/R/convGDX2MIF.R @@ -16,8 +16,9 @@ #' @param testthat boolean whether called by tests, turns some messages into warnings #' @author Lavinia Baumstark #' @examples -#' -#' \dontrun{convGDX2MIF(gdx,gdx_refpolicycost,file="REMIND_generic_default.csv",scenario="default")} +#' \dontrun{ +#' convGDX2MIF(gdx, gdx_refpolicycost, file = "REMIND_generic_default.csv", scenario = "default") +#' } #' #' @export #' @importFrom dplyr %>% bind_rows filter @@ -30,13 +31,12 @@ convGDX2MIF <- function(gdx, gdx_ref = NULL, file = NULL, scenario = "default", t = c(seq(2005, 2060, 5), seq(2070, 2110, 10), 2130, 2150), gdx_refpolicycost = gdx_ref, testthat = FALSE) { - # Define region subsets ---- regionSubsetList <- toolRegionSubsets(gdx) # ADD EU-27 region aggregation if possible - if("EUR" %in% names(regionSubsetList)){ - regionSubsetList <- c(regionSubsetList,list( - "EU27"=c("ENC","EWN","ECS","ESC","ECE","FRA","DEU","ESW") + if ("EUR" %in% names(regionSubsetList)) { + regionSubsetList <- c(regionSubsetList, list( + "EU27" = c("ENC", "EWN", "ECS", "ESC", "ECE", "FRA", "DEU", "ESW") )) } @@ -44,26 +44,26 @@ convGDX2MIF <- function(gdx, gdx_ref = NULL, file = NULL, scenario = "default", output <- NULL message("running reportMacroEconomy...") - output <- mbind(output,reportMacroEconomy(gdx,regionSubsetList,t)[,t,]) + output <- mbind(output, reportMacroEconomy(gdx, regionSubsetList, t)[, t, ]) message("running reportTrade...") - output <- mbind(output,reportTrade(gdx,regionSubsetList,t)[,t,]) + output <- mbind(output, reportTrade(gdx, regionSubsetList, t)[, t, ]) message("running reportPE...") - output <- mbind(output,reportPE(gdx,regionSubsetList,t)[,t,]) + output <- mbind(output, reportPE(gdx, regionSubsetList, t)[, t, ]) message("running reportSE...") - output <- mbind(output,reportSE(gdx,regionSubsetList,t)[,t,]) + output <- mbind(output, reportSE(gdx, regionSubsetList, t)[, t, ]) message("running reportFE...") - output <- mbind(output,reportFE(gdx,regionSubsetList,t)) + output <- mbind(output, reportFE(gdx, regionSubsetList, t)) message("running reportExtraction...") - output <- mbind(output,reportExtraction(gdx,regionSubsetList,t)[,t,]) + output <- mbind(output, reportExtraction(gdx, regionSubsetList, t)[, t, ]) message("running reportCapacity...") - output <- mbind(output,reportCapacity(gdx,regionSubsetList,t)[,t,]) - #output <- mbind(output,reportLCOE(gdx)[,t,]) now moved to additional LCOE.mif file because many variables + output <- mbind(output, reportCapacity(gdx, regionSubsetList, t, gdx_ref = gdx_ref)[, t, ]) + # output <- mbind(output,reportLCOE(gdx)[,t,]) now moved to additional LCOE.mif file because many variables message("running reportCapitalStock...") - output <- mbind(output,reportCapitalStock(gdx,regionSubsetList,t)[,t,]) + output <- mbind(output, reportCapitalStock(gdx, regionSubsetList, t, gdx_ref = gdx_ref)[, t, ]) message("running reportEnergyInvestment...") - output <- mbind(output,reportEnergyInvestment(gdx,regionSubsetList,t)[,t,]) + output <- mbind(output, reportEnergyInvestment(gdx, regionSubsetList, t, gdx_ref = gdx_ref)[, t, ]) message("running reportEmiAirPol...") - tmp <- try(reportEmiAirPol(gdx,regionSubsetList,t)) # test whether reportEmiAirPol works + tmp <- try(reportEmiAirPol(gdx, regionSubsetList, t)) # test whether reportEmiAirPol works if (!inherits(tmp, "try-error")) { if (!is.null(tmp)) output <- mbind(output, tmp[, t, ]) } else { @@ -71,21 +71,28 @@ convGDX2MIF <- function(gdx, gdx_ref = NULL, file = NULL, scenario = "default", } # reporting of variables that need variables from different other report functions - message("running reportEmi...") - output <- mbind(output,reportEmi(gdx,output,regionSubsetList,t)[,t,]) # needs output from reportFE + message("running reportEmi...") # needs output from reportFE + output <- mbind(output, reportEmi(gdx, output, regionSubsetList, t)[, t, ]) + message("running reportTechnology...") - output <- mbind(output,reportTechnology(gdx,output,regionSubsetList,t)[,t,]) # needs output from reportSE + # needs output from reportSE + output <- mbind(output, reportTechnology(gdx, output, regionSubsetList, t)[, t, ]) + message("running reportPrices...") - output <- mbind(output,reportPrices(gdx,output,regionSubsetList,t,gdx_ref = gdx_ref)[,t,]) # needs output from reportSE, reportFE, reportEmi, reportExtraction, reportMacroEconomy + # needs output from reportSE, reportFE, reportEmi, reportExtraction, reportMacroEconomy + output <- mbind(output, reportPrices(gdx, output, regionSubsetList, t, gdx_ref = gdx_ref)[, t, ]) + message("running reportCosts...") - output <- mbind(output,reportCosts(gdx,output,regionSubsetList,t)[,t,]) # needs output from reportEnergyInvestment, reportPrices, reportEnergyInvestments + # needs output from reportEnergyInvestment, reportPrices, reportEnergyInvestments + output <- mbind(output, reportCosts(gdx, output, regionSubsetList, t)[, t, ]) + message("running reportTax...") - output <- mbind(output,reportTax(gdx,output,regionSubsetList,t)[,t,]) + output <- mbind(output, reportTax(gdx, output, regionSubsetList, t)[, t, ]) # cross variables ---- # needs variables from different other report* functions message("running reportCrossVariables...") - output <- mbind(output,reportCrossVariables(gdx,output,regionSubsetList,t)[,t,]) + output <- mbind(output, reportCrossVariables(gdx, output, regionSubsetList, t)[, t, ]) # policy costs, if possible and sensible ---- if (is.null(gdx_refpolicycost)) { @@ -94,7 +101,7 @@ convGDX2MIF <- function(gdx, gdx_ref = NULL, file = NULL, scenario = "default", if (file.exists(gdx_refpolicycost)) { gdp_scen <- try(readGDX(gdx, "cm_GDPscen", react = "error"), silent = TRUE) gdp_scen_ref <- try(readGDX(gdx_refpolicycost, "cm_GDPscen", react = "error"), silent = TRUE) - if (! inherits(gdp_scen, "try-error") && ! inherits(gdp_scen_ref, "try-error")) { + if (!inherits(gdp_scen, "try-error") && !inherits(gdp_scen_ref, "try-error")) { if (gdp_scen[1] == gdp_scen_ref[1]) { if (gdx == gdx_refpolicycost) { msg_refpc <- "reporting 0 everywhere" @@ -102,7 +109,7 @@ convGDX2MIF <- function(gdx, gdx_ref = NULL, file = NULL, scenario = "default", msg_refpc <- paste0("comparing to ", basename(dirname(gdx_refpolicycost)), "/", basename(gdx_refpolicycost), "...") } message("running reportPolicyCosts, ", msg_refpc) - output <- mbind(output, reportPolicyCosts(gdx, gdx_refpolicycost, regionSubsetList, t)[,t,]) + output <- mbind(output, reportPolicyCosts(gdx, gdx_refpolicycost, regionSubsetList, t)[, t, ]) } else { warning("The GDP scenario differs from that of the reference run. Did not execute 'reportPolicyCosts'! ", "If a policy costs reporting is desired, please use the 'policyCosts' output.R script.") @@ -119,9 +126,9 @@ convGDX2MIF <- function(gdx, gdx_ref = NULL, file = NULL, scenario = "default", # SDP variables ---- message("running reportSDPVariables...") - tmp <- try(reportSDPVariables(gdx,output,t)) # test whether reportSDPVariables works + tmp <- try(reportSDPVariables(gdx, output, t)) # test whether reportSDPVariables works if (!inherits(tmp, "try-error")) { - if(!is.null(tmp)) output <- tmp + if (!is.null(tmp)) output <- tmp } else { message("function reportSDPVariables does not work and is skipped") } @@ -133,16 +140,16 @@ convGDX2MIF <- function(gdx, gdx_ref = NULL, file = NULL, scenario = "default", # clean and test output ---- # Add dimension names "scenario.model.variable" getSets(output)[3] <- "variable" - output <- add_dimension(output,dim=3.1,add = "model",nm = "REMIND") - output <- add_dimension(output,dim=3.1,add = "scenario",nm = scenario) + output <- add_dimension(output, dim = 3.1, add = "model", nm = "REMIND") + output <- add_dimension(output, dim = 3.1, add = "scenario", nm = scenario) ## check variable names ---- checkVarNames(getNames(output, dim = 3)) ## summation checks ---- .reportSummationErrors <- function(msg, testthat) { - if (!any(grepl('All summation checks were fine', msg))) { - msgtext <- paste(msg, collapse = '\n') + if (!any(grepl("All summation checks were fine", msg))) { + msgtext <- paste(msg, collapse = "\n") if (isTRUE(testthat)) warning("### Analyzing ", basename(gdx), ":\n", msgtext) else message(msgtext) } } @@ -153,16 +160,16 @@ convGDX2MIF <- function(gdx, gdx_ref = NULL, file = NULL, scenario = "default", summationsFile = "extractVariableGroups", absDiff = 0.01, relDiff = 0.02, roundDiff = TRUE ), - type = 'message') %>% + type = "message") %>% .reportSummationErrors(testthat = testthat) capture.output(sumChecks <- checkSummations( mifFile = output, dataDumpFile = NULL, outputDirectory = NULL, - summationsFile = system.file('extdata/additional_summation_checks.csv', - package = 'remind2'), + summationsFile = system.file("extdata/additional_summation_checks.csv", + package = "remind2"), absDiff = 0.01, relDiff = 0.02, roundDiff = TRUE) %>% - bind_rows(sumChecks), - type = 'message' + bind_rows(sumChecks), + type = "message" ) %>% .reportSummationErrors(testthat = testthat) @@ -174,7 +181,7 @@ convGDX2MIF <- function(gdx, gdx_ref = NULL, file = NULL, scenario = "default", "^Emi\\|CO2\\|Energy\\|Demand\\|Industry\\|.*Fossil \\(Mt CO2/yr\\)$", low = 0), list("Share.*\\((%|Percent)\\)$", low = 0, up = 100)), - reaction = 'warning') + reaction = "warning") # write or return output ---- if (!is.null(file)) { @@ -184,18 +191,17 @@ convGDX2MIF <- function(gdx, gdx_ref = NULL, file = NULL, scenario = "default", # write additional file on summation errors if needed if (nrow(sumChecks) > 0) { - summation_errors_file <- sub('(\\.[^.]+)$', '_summation_errors.csv', file) + summation_errors_file <- sub("(\\.[^.]+)$", "_summation_errors.csv", file) warning("Summation checks have revealed some gaps! See file ", summation_errors_file) write.csv(sumChecks, summation_errors_file, quote = FALSE, row.names = FALSE) } - } - else { + } else { # return summation errors as attribute if (nrow(sumChecks) > 0) { warning("Summation checks have revealed some gaps! ", "See `summation_errors` attribute on output for details.") - attr(output, 'summation_errors') <- sumChecks + attr(output, "summation_errors") <- sumChecks } return(output) } diff --git a/R/convGDX2MIF_LCOE.R b/R/convGDX2MIF_LCOE.R index 3d651b98..87a2b728 100644 --- a/R/convGDX2MIF_LCOE.R +++ b/R/convGDX2MIF_LCOE.R @@ -5,7 +5,6 @@ #' #' #' @param gdx a GDX as created by readGDX, or the file name of a gdx -#' @param gdx_ref a GDX as created by readGDX of the reference run #' @param file name of the mif file which will be written, if no name is #' provided a magpie object containing all the reporting information is #' returned @@ -15,7 +14,7 @@ #' @author Lavinia Baumstark #' @examples #' \dontrun{ -#' convGDX2MIF(gdx, gdx_ref, file = "REMIND_generic_LCOE.csv", scenario = "default") +#' convGDX2MIF_LCOE(gdx, file = "REMIND_generic_LCOE.csv", scenario = "default") #' } #' #' @export @@ -23,11 +22,11 @@ #' @importFrom magclass mbind write.report #' @importFrom utils write.table -convGDX2MIF_LCOE <- function(gdx, gdx_ref, file = NULL, scenario = "default", +convGDX2MIF_LCOE <- function(gdx, file = NULL, scenario = "default", t = c(seq(2005, 2060, 5), seq(2070, 2110, 10), 2130, 2150)) { # make the reporting output <- NULL - output <- mbind(output, reportLCOE(gdx)[, t, ]) + output <- mbind(output, reportLCOE(gdx = gdx)[, t, ]) # write the LCOE.mif or give back the magpie object output if (!is.null(file)) { diff --git a/R/modifyInvestmentVariables.R b/R/modifyInvestmentVariables.R new file mode 100644 index 00000000..f1c088b5 --- /dev/null +++ b/R/modifyInvestmentVariables.R @@ -0,0 +1,106 @@ +#' Modify Investment Variables +#' +#' This function transforms investment variables to the normal reporting convention. +#' +#' Years represented by investment variables in the energy system +#' ('vm_deltaCap', 'vm_costInvTeDir' and 'vm_costInvTeAdj') are different from the normal +#' reporting convention. In the current REMIND version, vm_deltacap(t) represents +#' the average of the years t-4..t, while in the general reporting +#' convention it represents the average of t-2.5..t+2.5 (for 5 year time steps). +#' +#' See also: https://github.com/remindmodel/remind/pull/1238. +#' +#' +#' @param x a magclass object to be manipulated, must have timesteps in 'ttot' +#' @param ref an optional magclass object to be used for fixing values before 'startYear' +#' @param startYear years before will be overwritten with values from 'ref' +#' +#' +#' @author Falk Benke +#' +#' +modifyInvestmentVariables <- function(x, ref = NULL, startYear = NULL) { + + ttot <- c(seq(1900, 2060, 5), seq(2070, 2110, 10), 2130, 2150) + + if (!setequal(getYears(x, as.integer = TRUE), ttot)) { + stop("Timesteps must equal to 'ttot'") + } + + # generate mapping from REMIND timesteps to yearly timesteps for REMIND investment variables (vm_deltaCap etc.) + # REMIND investment variables are defined to cover years between the last and the current REMIND time step. + # Example: vm_deltaCap(2020) refers to annual capacity additions from 2016 to 2020. + investTs <- rbind( + # first year is 1898, as 1990 represents 1898 to 1902 in REMIND + data.frame(year = seq(1898, 2110, 1)) %>% + mutate("period" = ifelse( + .data$year <= 2060, + ifelse(.data$year %% 5 == 0, .data$year, .data$year + 5 - (.data$year %% 5)), + ifelse(.data$year %% 10 == 0, .data$year, .data$year + 10 - (.data$year %% 10)) + )), + # highest year is 2167, as 2150 represents 2140 - 2167 in REMIND + data.frame(year = seq(2111, 2167, 1)) %>% + mutate("period" = ifelse(.data$year <= 2130, 2130, 2150)) + ) + + # generate mapping from yearly timesteps to REMIND reporting timesteps. + # Reported variables are defined as averages between t-2 and t+2 around the central year t. + # Example: New Cap|XYZ(2020) refers to average annual capacity additions of technology XYZ in 2018 to 2022. + remindTs <- + rbind( + # mapping for 2005 - 2150 from quitte + quitte::remind_timesteps, + # enhance the logic to cover all years starting from 1898 (e.g. period 2000 represents 1998 - 2002) + data.frame(year = seq(1898, 2002, 1), weight = 1) %>% + mutate("period" = ifelse(.data$year %% 10 %in% c(1, 2, 8, 9), round(.data$year, digits = -1), + ifelse(.data$year %% 10 %in% c(6, 7), .data$year - (.data$year %% 5), + ifelse(.data$year %% 10 %in% c(3, 4), .data$year + 5 - (.data$year %% 5), .data$year) + ) + )) + ) + + investTs <- investTs %>% + mutate( + "year" = paste0("y", .data$year), + "period" = paste0("y", .data$period), + ) + + remindTs <- remindTs %>% + mutate( + "year" = paste0("y", .data$year), + "period" = paste0("y", .data$period), + ) + + # map variables from model timesteps to yearly timesteps (e.g. 2020 -> 2016-2020) + x <- toolAggregate(x, dim = 2, rel = investTs, from = "period", to = "year", verbosity = 2) + + w <- remindTs %>% + select(-"period") %>% + # period x weight combos should be distinct for this to work + distinct() %>% + as.magpie() + + # Average variables with yearly timesteps to 5-year reporting time steps defined + # around center year (e.g. 2018-2022 average -> 2020) + x <- toolAggregate(x, dim = 2, rel = remindTs, weight = w, from = "year", to = "period") + + if (!is.null(ref)) { + if (!all( + setequal(getYears(x), getYears(ref)), + setequal(getItems(x, dim = 1), getItems(ref, dim = 1)), + setequal(getNames(x), getNames(ref)) + )) { + stop("ref does not match the dimensions of x") + } + + fixedYears <- getYears(x, as.integer = TRUE)[getYears(x, as.integer = TRUE) < startYear] + if (length(fixedYears) == 0) { + return(x) + } + + ref <- modifyInvestmentVariables(ref) + x[, fixedYears, ] <- ref[, fixedYears, ] + } + + return(x) +} diff --git a/R/plotLCOE.R b/R/plotLCOE.R index 68684856..80890b1d 100644 --- a/R/plotLCOE.R +++ b/R/plotLCOE.R @@ -1,105 +1,43 @@ #' Read in LCOE mif and write LCOE_plots.pdf -#' +#' #' Read in all information from LCOE mif file and create #' the LCOE_plots.pdf -#' +#' #' @param LCOEfile a path to the LCOE reporting csv file of the scenario to be plotted #' @param gdx gdx file of run -#' @param y time span for the data in line plots, default: y=c(2015, 2020,2030,2040,2050,2060) +#' @param y time span for the data in line plots, default: y=c(2015, 2020,2030,2040,2050,2060) #' @param reg region(s) in focus, reg ="all_regi" shows all regions if the mifs contain different regions -#' @param fileName name of the pdf, default = "LCOE_plots.pdf" -#' +#' @param fileName name of the pdf, default = "LCOE_plots.pdf" +#' #' @author Felix Schreyer #' @export -#' @importFrom magclass read.report mbind +#' @importFrom magclass read.report mbind #' @importFrom lusweave swopen swlatex swfigure swclose -#' @importFrom ggplot2 ggplot facet_wrap geom_errorbar ggtitle xlab scale_y_continuous scale_fill_discrete geom_col aes element_text theme_bw sec_axis scale_x_discrete scale_linetype_identity unit +#' @importFrom ggplot2 ggplot facet_wrap geom_errorbar ggtitle scale_y_continuous +#' geom_col aes element_text theme_bw sec_axis scale_x_discrete scale_linetype_identity unit #' @importFrom dplyr left_join -#' @importFrom quitte order.levels sum_total getRegs revalue.levels getScenarios +#' @importFrom quitte order.levels getRegs revalue.levels getScenarios #' @importFrom tidyr spread gather #' @importFrom data.table frollmean #' @importFrom gdx readGDX -plotLCOE <- function(LCOEfile, gdx, y=c(2015,2020,2030,2040,2050,2060),reg="all_regi",fileName="LCOE_plots.pdf") { - - - df.LCOE.in <- read.csv(LCOEfile, sep = ";") - - # # read LCOE reporting files - # df.LCOE.fromFiles <- NULL - # for (i in 1:length(reportFiles)) { - # df.LCOE.in <- read.csv(reportFiles[i], sep = ";") - # df.LCOE.fromFiles <- rbind(df.LCOE.fromFiles, df.LCOE.in) - # } - - - # # read model results - # data <- NULL - # for(i in 1:length(mif)){ - # data_new <- read.report(mif[i],as.list=FALSE) - # if (magclass::getNames(data_new,fulldim = TRUE)[["scenario"]] %in% magclass::getNames(data,fulldim = TRUE)[["scenario"]]) magclass::getNames(data_new) <- gsub(magclass::getNames(data_new,fulldim = TRUE)["scenario"],paste0(magclass::getNames(data_new,fulldim = TRUE)["scenario"],i),magclass::getNames(data_new)) - # if(all(getRegions(data) %in% getRegions(data_new))) { - # data <- mbind(data,data_new) - # } else { - # if(is.null(reg)){ - # stop("the regional aggregation of the results are different, you might use reg='all_regi'") - # } else if(reg=="all_regi"){ - # if(all(getRegions(data_new) %in% getRegions(data))) { - # # expand data_new by old regions from data - # oldReg <- getRegions(data)[-which(getRegions(data) %in% getRegions(data_new))] - # dummy_data_new <- new.magpie(oldReg,getYears(data_new),getNames(data_new),fill=NA) - # data_new <- mbind(data_new,dummy_data_new) - # # compine old and new data - # data <- mbind(data,data_new) - # } else { - # # expand data by new regions from data_new - # newReg <- getRegions(data_new)[-which(getRegions(data_new) %in% getRegions(data))] - # dummy_data <- new.magpie(newReg,getYears(data),getNames(data),fill=NA) - # data <- mbind(data,dummy_data) - # # expand data_new by old regions from data - # oldReg <- getRegions(data)[-which(getRegions(data) %in% getRegions(data_new))] - # dummy_data_new <- new.magpie(oldReg,getYears(data_new),getNames(data_new),fill=NA) - # data_new <- mbind(data_new,dummy_data_new) - # # compine old and new data - # data <- mbind(data,data_new) - # } - # - # } else { - # stop("the regional aggregation of the results are different, you might use reg='all_regi'") - # } - # } - # } - # - # if (!(is.null(reg))) { - # if (!reg=="all_regi") { - # data <- data[reg,y,] - # } else { - # data <- data[,y,] - # } - # } else { - # data <- data[,y,] - # } - # - - +plotLCOE <- function(LCOEfile, gdx, y = c(2015, 2020, 2030, 2040, 2050, 2060), + reg = "all_regi", fileName = "LCOE_plots.pdf") { + df.LCOE.in <- read.csv(LCOEfile, sep = ";") + # initialize variables used in dplyr pipes and ggplot below - scenario <- NULL + tech <- NULL period <- NULL region <- NULL variable <- NULL value <- NULL - Total <- NULL type <- NULL - vm_deltaCap <- NULL LCOE <- NULL all_te <- NULL output <- NULL - aes <- NULL - element_text <- NULL cost <- NULL - Total <- NULL `CO2 Tax Cost` <- NULL `Second Fuel Cost` <- NULL TooHigh <- NULL @@ -107,238 +45,230 @@ plotLCOE <- function(LCOEfile, gdx, y=c(2015,2020,2030,2040,2050,2060),reg="all_ all_enty <- NULL Price <- NULL plot.tech <- NULL - - - # ---- plot settings ---- - template <- c("\\documentclass[a4paper,landscape]{article}", - "\\setlength{\\oddsidemargin}{-0.8in}", - "\\setlength{\\evensidemargin}{-0.5in}", - "\\setlength{\\topmargin}{-0.8in}", - "\\setlength{\\parindent}{0in}", - "\\setlength{\\headheight}{0in}", - "\\setlength{\\topskip}{0in}", - "\\setlength{\\headsep}{0in}", - "\\setlength{\\footskip}{0.2in}", - "\\setlength\\textheight{0.95\\paperheight}", - "\\setlength\\textwidth{0.95\\paperwidth}", - "\\setlength{\\parindent}{0in}", - "\\usepackage{float}", - "\\usepackage[bookmarksopenlevel=section,colorlinks=true,linkbordercolor={0.9882353 0.8352941 0.7098039}]{hyperref}", - "\\hypersetup{bookmarks=true,pdfauthor={GES group, PIK}}", - "\\usepackage{graphicx}", - "\\usepackage[strings]{underscore}", - "\\usepackage{Sweave}", - "\\begin{document}", - "<>=", - "options(width=110)", - "@") - - + template <- c( + "\\documentclass[a4paper,landscape]{article}", + "\\setlength{\\oddsidemargin}{-0.8in}", + "\\setlength{\\evensidemargin}{-0.5in}", + "\\setlength{\\topmargin}{-0.8in}", + "\\setlength{\\parindent}{0in}", + "\\setlength{\\headheight}{0in}", + "\\setlength{\\topskip}{0in}", + "\\setlength{\\headsep}{0in}", + "\\setlength{\\footskip}{0.2in}", + "\\setlength\\textheight{0.95\\paperheight}", + "\\setlength\\textwidth{0.95\\paperwidth}", + "\\setlength{\\parindent}{0in}", + "\\usepackage{float}", + "\\usepackage[bookmarksopenlevel=section,colorlinks=true,linkbordercolor={0.9882353 0.8352941 0.7098039}]{hyperref}", + "\\hypersetup{bookmarks=true,pdfauthor={GES group, PIK}}", + "\\usepackage{graphicx}", + "\\usepackage[strings]{underscore}", + "\\usepackage{Sweave}", + "\\begin{document}", + "<>=", + "options(width=110)", + "@" + ) + + ylimit.up <- 250 # y axis max ylimit.lo <- -150 # y axes min - plot.cost <- c( "Flex Tax","Second Fuel Cost", "CO2 Tax Cost", "Fuel Cost", - "OMV Cost" , "OMF Cost", "Investment Cost","Total LCOE") + plot.cost <- c( + "Flex Tax", "Second Fuel Cost", "CO2 Tax Cost", "Fuel Cost", + "OMV Cost", "OMF Cost", "Investment Cost", "Total LCOE" + ) plot.period <- y plot.scen <- getScenarios(df.LCOE.in) - #plot.period <- c(2020,2030,2040,2050) - plot_theme <- theme_bw() + theme(axis.text.x = element_text(angle=90, size=40, vjust=0.3), - text = element_text(size=50), - axis.text.y = element_text(size=60), - legend.text = element_text(size=50), - legend.key.size = unit(6,"line")) - - - # plot_theme <- theme_bw() + theme(axis.text.x = element_text(angle=90, size=20, vjust=0.3), - # text = element_text(size=20), - # axis.text.y = element_text(size=20), - # legend.text = element_text(size=20), - # legend.key.size = unit(2,"line")) - # colors of cost components - # cost.colors <- c("Investment Cost" = "azure4", "OMF Cost" = "darkcyan", "OMV Cost" = "cyan", - # "Fuel Cost" = "orange3", "CO2 Tax Cost" = "indianred", "CO2 Provision Cost" = "slategray1", - # "Second Fuel Cost" = "sandybrown", "VRE Storage Cost" = "mediumpurple3", - # "Grid Cost" = "darkseagreen3", "CCS Tax Cost" = "chartreuse4", - # "CCS Cost" = "darkblue", "Flex Tax" = "yellow1") - - cost.colors <- c("Investment Cost" = "azure4", "OMF Cost" = "darkcyan", "OMV Cost" = "cyan", - "Fuel Cost" = "orange3", "CO2 Tax Cost" = "indianred", - "Second Fuel Cost" = "sandybrown", "Flex Tax" = "yellow1") - - + + plot_theme <- theme_bw() + theme( + axis.text.x = element_text(angle = 90, size = 40, vjust = 0.3), + text = element_text(size = 50), + axis.text.y = element_text(size = 60), + legend.text = element_text(size = 50), + legend.key.size = unit(6, "line") + ) + + + cost.colors <- c( + "Investment Cost" = "azure4", "OMF Cost" = "darkcyan", "OMV Cost" = "cyan", + "Fuel Cost" = "orange3", "CO2 Tax Cost" = "indianred", + "Second Fuel Cost" = "sandybrown", "Flex Tax" = "yellow1" + ) + + # relabel outputs and only plot most important technologies per output - relabel.outputs <- c("seliqfos" = "seliq", "seliqbio" = "seliq", - "segafos" = "segas", "segabio" = "segas") - - plot.outputs <- c("seel","seliq","segas","seh2") - plot.techs <- c("pc", "igcc", "ngt","ngcc", "ngccc","tnrs","hydro","spv","wind","windoff","csp", - "biochp","bioigccc", - "refliq","biodiesel","bioftrec","bioftcrec", "MeOH", - "gastr", "biogas","h22ch4", - "elh2","elh2VRE", "bioh2", "bioh2c", "gash2c", "coalh2c") - + relabel.outputs <- c( + "seliqfos" = "seliq", "seliqbio" = "seliq", + "segafos" = "segas", "segabio" = "segas" + ) + + plot.outputs <- c("seel", "seliq", "segas", "seh2") + plot.techs <- c( + "pc", "igcc", "ngt", "ngcc", "ngccc", "tnrs", "hydro", "spv", "wind", "windoff", "csp", + "biochp", "bioigccc", + "refliq", "biodiesel", "bioftrec", "bioftcrec", "MeOH", + "gastr", "biogas", "h22ch4", + "elh2", "elh2VRE", "bioh2", "bioh2c", "gash2c", "coalh2c" + ) + ### end plot settings - + # read in capacity additions from .mif, for second y axes of plot - vm_deltaCap <- readGDX(gdx, "vm_deltaCap", field = "l", restore_zeros = F) - - + vm_deltaCap <- readGDX(gdx, "vm_deltaCap", field = "l", restore_zeros = FALSE) %>% + modifyInvestmentVariables() + # calculate 15-year moving average capacity additions - df.dC <- as.quitte(dimSums(vm_deltaCap, dim=3.2)*1e3) %>% - rename(tech=all_te, vm_deltaCap = value) %>% - filter(period >= 2005) %>% - select(region, period, tech, vm_deltaCap) %>% - group_by(region, tech) %>% - mutate( vm_deltaCap = frollmean(vm_deltaCap, 3, align = "center", fill = NA)) %>% - ungroup() - + df.dC <- as.quitte(dimSums(vm_deltaCap, dim = 3.2) * 1e3) %>% + rename(tech = all_te, vm_deltaCap = value) %>% + filter(period >= 2005) %>% + select(region, period, tech, vm_deltaCap) %>% + group_by(region, tech) %>% + mutate(vm_deltaCap = frollmean(vm_deltaCap, 3, align = "center", fill = NA)) %>% + ungroup() + # calculate prices - tdptwyr2dpgj <- 31.71 #TerraDollar per TWyear to Dollar per GJ - qBalSe.m <- readGDX(gdx, "q_balSe", field = "m", restore_zeros = F) - qBalSeel.m <- readGDX(gdx,"q32_balSe",field="m", restore_zeros = F) - budget.m <- readGDX(gdx,'qm_budget',field = "m")[,getYears(qBalSe.m),] - vm_prodSe <- readGDX(gdx, "vm_prodSe", field = "l", restore_zeros = F) - - - - vm_prodSe_output <- dimSums(vm_prodSe, dim=c(3.1,3.3), na.rm = T) # sum SE over all technologies per output type + tdptwyr2dpgj <- 31.71 # TerraDollar per TWyear to Dollar per GJ + qBalSe.m <- readGDX(gdx, "q_balSe", field = "m", restore_zeros = FALSE) + qBalSeel.m <- readGDX(gdx, "q32_balSe", field = "m", restore_zeros = FALSE) + budget.m <- readGDX(gdx, "qm_budget", field = "m")[, getYears(qBalSe.m), ] + vm_prodSe <- readGDX(gdx, "vm_prodSe", field = "l", restore_zeros = FALSE) + + vm_prodSe_output <- dimSums(vm_prodSe, dim = c(3.1, 3.3), na.rm = TRUE) # sum SE over all technologies per output type qBalSe.m <- mbind(qBalSeel.m, qBalSe.m) # bind marginals from seel balanace to marginal from other se balance eqation - SePrice <- qBalSe.m / (budget.m+1e-10) * tdptwyr2dpgj * 1.2 * 3.66 # calculate SE Prices in USD2015/Mwh - + SePrice <- qBalSe.m / (budget.m + 1e-10) * tdptwyr2dpgj * 1.2 * 3.66 # calculate SE Prices in USD2015/Mwh + # calculate gases and liquids prices as demand-weighted average of biofuel and fossil fuel - SePrice <- mbind(SePrice, - setNames( - SePrice[,,"seliqfos"]*vm_prodSe_output[,,"seliqfos"]+ - SePrice[,,"seliqbio"]*vm_prodSe_output[,,"seliqbio"] / - (vm_prodSe_output[,,"seliqfos"]+vm_prodSe_output[,,"seliqbio"]), - "seliq"), - setNames( - SePrice[,,"segafos"]*vm_prodSe_output[,,"segafos"]+ - SePrice[,,"segabio"]*vm_prodSe_output[,,"segabio"] / - (vm_prodSe_output[,,"segafos"]+vm_prodSe_output[,,"segabio"]), - "segas")) - + SePrice <- mbind( + SePrice, + setNames( + SePrice[, , "seliqfos"] * vm_prodSe_output[, , "seliqfos"] + + SePrice[, , "seliqbio"] * vm_prodSe_output[, , "seliqbio"] / + (vm_prodSe_output[, , "seliqfos"] + vm_prodSe_output[, , "seliqbio"]), + "seliq" + ), + setNames( + SePrice[, , "segafos"] * vm_prodSe_output[, , "segafos"] + + SePrice[, , "segabio"] * vm_prodSe_output[, , "segabio"] / + (vm_prodSe_output[, , "segafos"] + vm_prodSe_output[, , "segabio"]), + "segas" + ) + ) + # convert SE prices to dataframe, calculate 15-year moving average - df.SePrice <- as.quitte(SePrice) %>% - rename(output = all_enty, Price = value) %>% - select(region, period, output, Price) %>% - # moving average - group_by(region, output) %>% - mutate( Price = frollmean(Price, 3, align = "center", fill = NA)) %>% - ungroup() - + df.SePrice <- as.quitte(SePrice) %>% + rename(output = all_enty, Price = value) %>% + select(region, period, output, Price) %>% + # moving average + group_by(region, output) %>% + mutate(Price = frollmean(Price, 3, align = "center", fill = NA)) %>% + ungroup() + # join LCOE with capacity additions and SE prices df.LCOE.dC.join <- df.LCOE.in %>% - filter( type == "marginal") %>% - revalue.levels(output = relabel.outputs) %>% - rename(LCOE = value) %>% - left_join(df.dC) %>% - left_join(df.SePrice) %>% - gather(variable, value, LCOE, vm_deltaCap, Price) %>% - # do away with cost dimension for non LCOE variables - filter( cost == "Investment Cost" | variable == "LCOE", - period %in% plot.period, - output %in% plot.outputs, - tech %in% plot.techs, - cost %in% plot.cost) %>% - mutate( tech = as.factor(tech)) %>% - order.levels(plot.tech) - - - df.LCOE.dC.join <- df.LCOE.dC.join %>% - filter(region %in% c("USA")) + filter(type == "marginal") %>% + revalue.levels(output = relabel.outputs) %>% + rename(LCOE = value) %>% + left_join(df.dC) %>% + left_join(df.SePrice) %>% + gather(variable, value, LCOE, vm_deltaCap, Price) %>% + # do away with cost dimension for non LCOE variables + filter( + cost == "Investment Cost" | variable == "LCOE", + period %in% plot.period, + output %in% plot.outputs, + tech %in% plot.techs, + cost %in% plot.cost + ) %>% + mutate(tech = as.factor(tech)) %>% + order.levels(plot.tech) + + + df.LCOE.dC.join <- df.LCOE.dC.join %>% + filter(region %in% c("USA")) ### loop to plot LCOE and capacity additions per region and energy output - - # have to do plotting loop with lapply such that variable scale of second y axes works + + # have to do plotting loop with lapply such that variable scale of second y axes works # (normal loop cannot do this) LCOE.grouped <- split(df.LCOE.dC.join, list(df.LCOE.dC.join$output, df.LCOE.dC.join$region)) make_plot <- function(df.LCOE.grouped) { - # region and output plot.reg <- getRegs(df.LCOE.grouped) plot.output <- unique(df.LCOE.grouped$output) - + print(plot.reg) print(plot.output) - + # remove LCOE that are higher than y axis limit - df.LCOE.plot <- df.LCOE.grouped %>% - select(period, tech, cost, variable, value) %>% - spread(cost, value) %>% - mutate( TooHigh = ifelse(variable == "LCOE", ifelse((`Total LCOE` - ifelse(`CO2 Tax Cost` < 0, `CO2 Tax Cost`, 0) - - ifelse(`Second Fuel Cost` < 0, `Second Fuel Cost`, 0)) > ylimit.up,T,F), F)) %>% - gather(cost, value, -period, -tech, -variable, -TooHigh) %>% - filter(cost == "Investment Cost" | variable == "LCOE", TooHigh == FALSE) - - - - # calculate lower y axis limit, limit of second y axes - df.LCOE.min <- df.LCOE.plot %>% - # the two cost components that can be negative - filter(cost %in% c("CO2 Tax Cost","Second Fuel Cost")) - - #ylimit.lo <- min(df.LCOE.min$value, na.rm = T)-50 # y axis min. - - + df.LCOE.plot <- df.LCOE.grouped %>% + select(period, tech, cost, variable, value) %>% + spread(cost, value) %>% + mutate(TooHigh = ifelse(variable == "LCOE", ifelse((`Total LCOE` - ifelse(`CO2 Tax Cost` < 0, `CO2 Tax Cost`, 0) + - ifelse(`Second Fuel Cost` < 0, `Second Fuel Cost`, 0)) > ylimit.up, TRUE, FALSE), FALSE)) %>% + gather(cost, value, -period, -tech, -variable, -TooHigh) %>% + filter(cost == "Investment Cost" | variable == "LCOE", TooHigh == FALSE) + # maximum capacity addition in plot - df.dC.max <- df.LCOE.plot %>% - filter(variable == "vm_deltaCap") - - sec.axis.limit <- max(df.dC.max$value, na.rm = T) + 1 # second axis max - - - - df.LCOE.plot.out <- df.LCOE.plot %>% - # scale second axis variable s.t it fits into plot - mutate(value = ifelse(variable == "vm_deltaCap", - value / sec.axis.limit * ylimit.up, - value)) %>% - # remove Total LCOE - filter( cost != "Total LCOE") %>% - order.levels(cost = plot.cost) - - df.LCOE.total <- df.LCOE.plot %>% - filter(cost == "Total LCOE") - - + df.dC.max <- df.LCOE.plot %>% + filter(variable == "vm_deltaCap") + + sec.axis.limit <- max(df.dC.max$value, na.rm = TRUE) + 1 # second axis max + + df.LCOE.plot.out <- df.LCOE.plot %>% + # scale second axis variable s.t it fits into plot + mutate(value = ifelse(variable == "vm_deltaCap", + value / sec.axis.limit * ylimit.up, + value + )) %>% + # remove Total LCOE + filter(cost != "Total LCOE") %>% + order.levels(cost = plot.cost) + + df.LCOE.total <- df.LCOE.plot %>% + filter(cost == "Total LCOE") + + p.LCOE <- ggplot() + - geom_col(data=df.LCOE.plot.out %>% filter(variable == "LCOE"), - aes(tech, value, fill=cost), alpha=0.7) + - geom_point(data=df.LCOE.plot.out %>% filter(variable == "vm_deltaCap"), - aes(tech, value, color=variable), - alpha=0.8, size=9) + - # geom_point(data=df.LCOE.plot.out %>% filter(variable == "vm_deltaCap"), - # aes(tech, value), - # alpha=0.9, size=5, color="white") + - geom_hline(data=df.LCOE.plot.out %>% filter(variable == "Price"), - aes(yintercept=value, color=variable), size=2, alpha=0.7) + - geom_errorbar(data=df.LCOE.total, aes(x=tech, ymin=value, ymax=value, linetype="solid"), size=2) + - scale_linetype_identity(name = '', guide = 'legend',labels = c('Total LCOE')) + - #scale_color_identity(name = '', guide = 'legend',labels = c('REMIND Price')) + - scale_color_manual(values=c("Price"="blue", "vm_deltaCap" = "firebrick")) + + geom_col( + data = df.LCOE.plot.out %>% filter(variable == "LCOE"), + aes(tech, value, fill = cost), alpha = 0.7 + ) + + geom_point( + data = df.LCOE.plot.out %>% filter(variable == "vm_deltaCap"), + aes(tech, value, color = variable), + alpha = 0.8, size = 9 + ) + + geom_hline( + data = df.LCOE.plot.out %>% filter(variable == "Price"), + aes(yintercept = value, color = variable), size = 2, alpha = 0.7 + ) + + geom_errorbar(data = df.LCOE.total, aes(x = tech, ymin = value, ymax = value, linetype = "solid"), size = 2) + + scale_linetype_identity(name = "", guide = "legend", labels = c("Total LCOE")) + + scale_color_manual(values = c("Price" = "blue", "vm_deltaCap" = "firebrick")) + facet_wrap(~period) + plot_theme + - scale_fill_manual(values=cost.colors, name="LCOE components") + + scale_fill_manual(values = cost.colors, name = "LCOE components") + scale_x_discrete("Technology") + scale_y_continuous("LCOE and REMIND Price (USD2015/MWh)", - limits = c(ylimit.lo,ylimit.up), - breaks = seq(round(ylimit.lo,digits = -1), - round(ylimit.up, digits = -1), - 50), - sec.axis = sec_axis(~ . * sec.axis.limit / ylimit.up, - name = paste0("Capacity Additions in GW/yr\n(15-year average)")))+ - ggtitle(paste0(plot.scen,", ",plot.reg,":\nMarginal LCOE and capacity additions of ", plot.output," technologies")) - - - - - swfigure(sw,print,p.LCOE, sw_option="height=20,width=35", dpi=1800) - + limits = c(ylimit.lo, ylimit.up), + breaks = seq( + round(ylimit.lo, digits = -1), + round(ylimit.up, digits = -1), + 50 + ), + sec.axis = sec_axis(~ . * sec.axis.limit / ylimit.up, + name = paste0("Capacity Additions in GW/yr\n(15-year average)") + ) + ) + + ggtitle(paste0(plot.scen, ", ", plot.reg, ":\nMarginal LCOE and capacity additions of ", + plot.output, " technologies")) + + + swfigure(sw, print, p.LCOE, sw_option = "height=20,width=35", dpi = 1800) } - - sw <- swopen(fileName,template = template) - fig_list <- lapply(LCOE.grouped, make_plot) - swclose(sw) + sw <- swopen(fileName, template = template) + lapply(LCOE.grouped, make_plot) + swclose(sw) } diff --git a/R/reportCapacity.R b/R/reportCapacity.R index b2c159c1..1b6fb251 100644 --- a/R/reportCapacity.R +++ b/R/reportCapacity.R @@ -10,6 +10,9 @@ #' be created. #' @param t temporal resolution of the reporting, default: #' t=c(seq(2005,2060,5),seq(2070,2110,10),2130,2150) +#' @param gdx_ref a GDX object as created by readGDX, or the path to a gdx of the reference run. +#' It is used to guarantee consistency before 'cm_startyear' for capacity variables +#' using time averaging. #' #' @return MAgPIE object - contains the capacity variables #' @author Lavinia Baumstark, Christoph Bertram @@ -24,7 +27,9 @@ #' @importFrom magclass mbind setNames getSets getSets<- as.magpie #' @importFrom dplyr %>% filter mutate -reportCapacity <- function(gdx, regionSubsetList = NULL, t = c(seq(2005, 2060, 5), seq(2070, 2110, 10), 2130, 2150)) { +reportCapacity <- function(gdx, regionSubsetList = NULL, + t = c(seq(2005, 2060, 5), seq(2070, 2110, 10), 2130, 2150), + gdx_ref = gdx_ref) { # read sets teall2rlf <- readGDX(gdx, name = c("te2rlf", "teall2rlf"), format = "first_found") possibleRefineries <- c("refped", "refdip", "refliq") @@ -43,13 +48,27 @@ reportCapacity <- function(gdx, regionSubsetList = NULL, t = c(seq(2005, 2060, 5 # read scalars sm_c_2_co2 <- as.vector(readGDX(gdx, "sm_c_2_co2")) - # data preparation ttot <- as.numeric(as.vector(ttot)) - vm_cap <- vm_cap[teall2rlf] - vm_cap <- vm_cap[, ttot, ] + vm_cap <- vm_cap[teall2rlf] + vm_cap <- vm_cap[, ttot, ] + vm_deltaCap <- vm_deltaCap[teall2rlf] vm_deltaCap <- vm_deltaCap[, ttot, ] + + # apply 'modifyInvestmentVariables' to shift from the model-internal time coverage (deltacap and investment + # variables for step t represent the average of the years from t-4years to t) to the general convention for + # the reporting template (all variables represent the average of the years from t-2.5years to t+2.5years) + if (!is.null(gdx_ref)) { + cm_startyear <- as.integer(readGDX(gdx, name = "cm_startyear", format = "simplest")) + vm_deltaCapRef <- readGDX(gdx_ref, name = c("vm_deltaCap"), field = "l", format = "first_found") * 1000 + vm_deltaCapRef <- vm_deltaCapRef[teall2rlf] + vm_deltaCapRef <- vm_deltaCapRef[, ttot, ] + vm_deltaCap <- modifyInvestmentVariables(vm_deltaCap, vm_deltaCapRef, cm_startyear) + } else { + vm_deltaCap <- modifyInvestmentVariables(vm_deltaCap) + } + v_earlyreti <- v_earlyreti[, ttot, ] t2005 <- ttot[ttot > 2004] @@ -154,9 +173,8 @@ reportCapacity <- function(gdx, regionSubsetList = NULL, t = c(seq(2005, 2060, 5 # carbon management if ("dac" %in% magclass::getNames(vm_cap, dim = 1)) { - tmp <- mbind(tmp, setNames(dimSums(vm_cap[, , c("dac")], dim = 3) * sm_c_2_co2, "Cap|Carbon Management|DAC (Mt CO2/yr)")) + tmp <- mbind(tmp, setNames(dimSums(vm_cap[, , c("dac")], dim = 3) * sm_c_2_co2, "Cap|Carbon Management|DAC (Mt CO2/yr)")) } - # Newly built capacities electricity (Should all go into tmp2, so that this can be used for calculating cumulated values in tmp5 below) tmp2 <- NULL tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , c("ngcc", "ngt", "gaschp", "ngccc")], dim = 3), "New Cap|Electricity|Gas (GW/yr)")) @@ -225,40 +243,38 @@ reportCapacity <- function(gdx, regionSubsetList = NULL, t = c(seq(2005, 2060, 5 tmp_chpAdd <- mbind(tmp_chpAdd, setNames(dimSums(tmp_chpAdd, dim = 3), "New Cap|Heat (GW/yr)")) tmp2 <- mbind(tmp2, tmp_chpAdd) - # Newly built capacities liquids - tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , c(refineries, "coalftrec", "coalftcrec", "bioftrec", "bioftcrec", "biodiesel", "bioeths", "bioethl", "MeOH")], dim = 3), - "New Cap|Liquids (GW/yr)")) - tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , c(refineries, "coalftrec", "coalftcrec")], dim = 3), - "New Cap|Liquids|Fossil (GW/yr)")) - tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , c("coalftrec", "coalftcrec")], dim = 3), - "New Cap|Liquids|Coal (GW/yr)")) - tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , refineries], dim = 3), - "New Cap|Liquids|Oil (GW/yr)")) - tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , c("bioftrec", "bioftcrec", "biodiesel", "bioeths", "bioethl")], dim = 3), - "New Cap|Liquids|Biomass (GW/yr)")) - tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , c("MeOH")], dim = 3), - "New Cap|Liquids|Hydrogen (GW/yr)")) - - - # Newly built capacities gases - tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , c("gastr", "coalgas", "biogas", "h22ch4")], dim = 3), - "New Cap|Gases (GW/yr)")) - tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , c("gastr", "coalgas")], dim = 3), - "New Cap|Gases|Fossil (GW/yr)")) - tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , "gastr"], dim = 3), - "New Cap|Gases|Gas (GW/yr)")) - tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , "coalgas"], dim = 3), - "New Cap|Gases|Coal (GW/yr)")) - tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , c("biogas")], dim = 3), - "New Cap|Gases|Biomass (GW/yr)")) - tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , c("h22ch4")], dim = 3), - "New Cap|Gases|Hydrogen (GW/yr)")) - - # carbon management - if ("dac" %in% magclass::getNames(vm_cap, dim = 1)) { - tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , c("dac")], dim = 3) * sm_c_2_co2, "New Cap|Carbon Management|DAC (Mt CO2/yr/yr)")) - } + # Newly built capacities liquids + tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , c(refineries, "coalftrec", "coalftcrec", "bioftrec", "bioftcrec", "biodiesel", "bioeths", "bioethl", "MeOH")], dim = 3), + "New Cap|Liquids (GW/yr)")) + tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , c(refineries, "coalftrec", "coalftcrec")], dim = 3), + "New Cap|Liquids|Fossil (GW/yr)")) + tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , c("coalftrec", "coalftcrec")], dim = 3), + "New Cap|Liquids|Coal (GW/yr)")) + tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , refineries], dim = 3), + "New Cap|Liquids|Oil (GW/yr)")) + tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , c("bioftrec", "bioftcrec", "biodiesel", "bioeths", "bioethl")], dim = 3), + "New Cap|Liquids|Biomass (GW/yr)")) + tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , c("MeOH")], dim = 3), + "New Cap|Liquids|Hydrogen (GW/yr)")) + + # Newly built capacities gases + tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , c("gastr", "coalgas", "biogas", "h22ch4")], dim = 3), + "New Cap|Gases (GW/yr)")) + tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , c("gastr", "coalgas")], dim = 3), + "New Cap|Gases|Fossil (GW/yr)")) + tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , "gastr"], dim = 3), + "New Cap|Gases|Gas (GW/yr)")) + tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , "coalgas"], dim = 3), + "New Cap|Gases|Coal (GW/yr)")) + tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , c("biogas")], dim = 3), + "New Cap|Gases|Biomass (GW/yr)")) + tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , c("h22ch4")], dim = 3), + "New Cap|Gases|Hydrogen (GW/yr)")) + # carbon management + if ("dac" %in% magclass::getNames(vm_cap, dim = 1)) { + tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , c("dac")], dim = 3) * sm_c_2_co2, "New Cap|Carbon Management|DAC (Mt CO2/yr/yr)")) + } # add terms calculated from previously calculated capacity values tmp_aux <- NULL @@ -276,7 +292,7 @@ reportCapacity <- function(gdx, regionSubsetList = NULL, t = c(seq(2005, 2060, 5 names_capacities <- intersect(names_capacities, getNames(tmp1)) tmp_aux <- mbind(tmp_aux, setNames(dimSums(tmp1[, , names_capacities], dim = 3) + 0.6 * tmp1[, , "Cap|Electricity|Hydro (GW)"] + tmp[, , "Cap|Electricity|Storage|Battery (GW)"], - "Cap|Electricity|Estimated firm capacity counting hydro at 0p6 (GW)")) + "Cap|Electricity|Estimated firm capacity counting hydro at 0p6 (GW)")) tmp1 <- mbind(tmp1, tmp_aux) @@ -287,19 +303,19 @@ reportCapacity <- function(gdx, regionSubsetList = NULL, t = c(seq(2005, 2060, 5 # Idle capacities and Total (sum of operating and idle) tmp4 <- NULL tmp4 <- mbind(tmp4, setNames(dimSums(vm_cap[, , "igcc"], dim = 3) * v_earlyreti[, , "igcc"] / (1 - v_earlyreti[, , "igcc"]) + - dimSums(vm_cap[, , "coalchp"], dim = 3) * v_earlyreti[, , "coalchp"] / (1 - v_earlyreti[, , "coalchp"]) + - dimSums(vm_cap[, , "pc"], dim = 3) * v_earlyreti[, , "pc"] / (1 - v_earlyreti[, , "pc"]), - "Idle Cap|Electricity|Coal|w/o CC (GW)")) + dimSums(vm_cap[, , "coalchp"], dim = 3) * v_earlyreti[, , "coalchp"] / (1 - v_earlyreti[, , "coalchp"]) + + dimSums(vm_cap[, , "pc"], dim = 3) * v_earlyreti[, , "pc"] / (1 - v_earlyreti[, , "pc"]), + "Idle Cap|Electricity|Coal|w/o CC (GW)")) tmp4 <- mbind(tmp4, setNames(dimSums(vm_cap[, , "ngcc"], dim = 3) * v_earlyreti[, , "ngcc"] / (1 - v_earlyreti[, , "ngcc"]) + - dimSums(vm_cap[, , "gaschp"], dim = 3) * v_earlyreti[, , "gaschp"] / (1 - v_earlyreti[, , "gaschp"]) + - dimSums(vm_cap[, , "ngt"], dim = 3) * v_earlyreti[, , "ngt"] / (1 - v_earlyreti[, , "ngt"]), - "Idle Cap|Electricity|Gas|w/o CC (GW)")) + dimSums(vm_cap[, , "gaschp"], dim = 3) * v_earlyreti[, , "gaschp"] / (1 - v_earlyreti[, , "gaschp"]) + + dimSums(vm_cap[, , "ngt"], dim = 3) * v_earlyreti[, , "ngt"] / (1 - v_earlyreti[, , "ngt"]), + "Idle Cap|Electricity|Gas|w/o CC (GW)")) tmp4 <- mbind(tmp4, setNames(dimSums(vm_cap[, , "dot"], dim = 3) * v_earlyreti[, , "dot"] / (1 - v_earlyreti[, , "dot"]), - "Idle Cap|Electricity|Oil|w/o CC (GW)")) + "Idle Cap|Electricity|Oil|w/o CC (GW)")) tmp4 <- mbind(tmp4, setNames(tmp4[, , "Idle Cap|Electricity|Coal|w/o CC (GW)"] + tmp[, , "Cap|Electricity|Coal|w/o CC (GW)"], - "Total Cap|Electricity|Coal|w/o CC (GW)")) + "Total Cap|Electricity|Coal|w/o CC (GW)")) tmp4 <- mbind(tmp4, setNames(tmp4[, , "Idle Cap|Electricity|Gas|w/o CC (GW)"] + tmp[, , "Cap|Electricity|Gas|w/o CC (GW)"], - "Total Cap|Electricity|Gas|w/o CC (GW)")) + "Total Cap|Electricity|Gas|w/o CC (GW)")) # Cumulate things on extensive time set tmp <- mbind(tmp, tmp7, tmp1, tmp2, tmp4) @@ -309,9 +325,9 @@ reportCapacity <- function(gdx, regionSubsetList = NULL, t = c(seq(2005, 2060, 5 tmp6 <- quitte::as.quitte(tmp6) mylist <- lapply(levels(tmp6$variable), function(x) { calcCumulatedDiscount(data = tmp6 %>% - filter(.data$variable == x), - nameVar = x, - discount = 0.0) %>% + filter(.data$variable == x), + nameVar = x, + discount = 0.0) %>% mutate(variable = gsub("New", replacement = "Cumulative", x)) }) @@ -325,7 +341,6 @@ reportCapacity <- function(gdx, regionSubsetList = NULL, t = c(seq(2005, 2060, 5 magclass::getNames(tmp6)[grep("Cumulative Cap\\|Carbon Management", getNames(tmp6))] <- gsub("GW", "Mt CO2/yr", magclass::getNames(tmp6)[grep("Cumulative Cap\\|Carbon Management", getNames(tmp6))]) - tmp <- mbind(tmp[, t2005, ], tmp6) # add global values tmp <- mbind(tmp, dimSums(tmp, dim = 1)) @@ -334,5 +349,6 @@ reportCapacity <- function(gdx, regionSubsetList = NULL, t = c(seq(2005, 2060, 5 tmp <- mbind(tmp, calc_regionSubset_sums(tmp, regionSubsetList)) getSets(tmp)[3] <- "variable" + return(tmp) } diff --git a/R/reportCapitalStock.R b/R/reportCapitalStock.R index b8c67d59..a29b4937 100644 --- a/R/reportCapitalStock.R +++ b/R/reportCapitalStock.R @@ -3,136 +3,155 @@ #' Read in capital stock information from GDX file, information used in convGDX2MIF.R #' for the reporting #' -#' #' @param gdx a GDX object as created by readGDX, or the path to a gdx #' @param regionSubsetList a list containing regions to create report variables region #' aggregations. If NULL (default value) only the global region aggregation "GLO" will #' be created. #' @param t temporal resolution of the reporting, default: #' t=c(seq(2005,2060,5),seq(2070,2110,10),2130,2150) +#' @param gdx_ref a GDX object as created by readGDX, or the path to a gdx of the reference run. +#' It is used to guarantee consistency before 'cm_startyear' for capacity variables +#' using time averaging. #' #' @return MAgPIE object - contains the capital stock variables #' @author Lavinia Baumstark; Michaja Pehl #' @seealso \code{\link{convGDX2MIF}} #' @examples -#' -#' \dontrun{reportCapitalStock(gdx)} +#' \dontrun{ +#' reportCapitalStock(gdx) +#' } #' #' @export #' @importFrom gdx readGDX #' @importFrom magclass getYears mbind setNames #' @importFrom dplyr tribble -reportCapitalStock <- function(gdx,regionSubsetList=NULL,t=c(seq(2005,2060,5),seq(2070,2110,10),2130,2150)) { - - module2realisation <- readGDX(gdx, "module2realisation", react = "silent") - tran_mod = module2realisation[module2realisation$modules == "transport", 2] - - pm_conv_cap_2_MioLDV <- 650 # The world has ~715million cars in 2005 (IEA TECO2) - - # read sets - teall2rlf <- readGDX(gdx,name=c("te2rlf","teall2rlf"),format="first_found") - teue2rlf <- readGDX(gdx,name=c("teue2rlf", "tees2rlf"),format="first_found") - te <- readGDX(gdx,name=c("te"),format="first_found") - # read variables - vm_cap <- readGDX(gdx,name=c("vm_cap"),field="l",format="first_found") - vm_deltaCap <- readGDX(gdx,name=c("vm_deltaCap"),field="l",format="first_found") - v_investcost <- readGDX(gdx,name=c("vm_costTeCapital","v_costTeCapital","v_investcost"),field="l",format="first_found") - vm_cesIO <- readGDX(gdx, name = 'vm_cesIO', field = 'l') - # read parameters - ppfKap_Ind <- readGDX(gdx, name = 'ppfkap_industry_dyn37', react = 'silent') - steel_process_based <- "steel" %in% readGDX(gdx, "secInd37Prc", react='silent') - - # calculate maximal temporal resolution - y <- Reduce(intersect,list(getYears(vm_cap),getYears(v_investcost))) - vm_cap <- vm_cap[,y,] - vm_deltaCap <- vm_deltaCap[,y,] - v_investcost <- v_investcost[,y,] - - tmp <- NULL - - # ---- report transport capital stocks ---- - if (tran_mod == "complex"){ - LDV35 <- readGDX(gdx,name=c("LDV35"),format="first_found") - tmp <- mbind(tmp,setNames(dimSums( (vm_cap * v_investcost)[teue2rlf] - ,dim=c(3.1,3.2)) * 1000, "Est Capital Stock|ESM|Transp vehic (billion US$2017)")) - tmp <- mbind(tmp,setNames(dimSums( (vm_cap * v_investcost)[teall2rlf][,,LDV35] - ,dim=c(3.1,3.2)) * 1000, "Est Capital Stock|ESM|Pet/EV LDV (billion US$2017)")) - - tmp <- mbind(tmp, - setNames( - dimSums(mbind(vm_cap * v_investcost), - dim = 3) * 1000, - "Estimated Capital Stock|ESM (billion US$2017)" - ) - ) - - # prepare variables - vm_cap <- vm_cap[teall2rlf] - vm_deltaCap <- vm_deltaCap[teall2rlf] - - # build reporting - tmp <- mbind(tmp,setNames(dimSums(vm_cap[,,LDV35] * pm_conv_cap_2_MioLDV,dim=c(3.1,3.2)),"Est LDV Stock (million vehicles)")) - tmp <- mbind(tmp,setNames(dimSums(vm_cap[,,"apCarElT"] * pm_conv_cap_2_MioLDV,dim=c(3.1,3.2)),"Est EV LDV Stock (million vehicles)")) - tmp <- mbind(tmp,setNames(dimSums(vm_cap[,,"apCarH2T"] * pm_conv_cap_2_MioLDV,dim=c(3.1,3.2)),"Est H2 LDV Stock (million vehicles)")) - tmp <- mbind(tmp,setNames(dimSums(vm_cap[,,"apCarPeT"] * pm_conv_cap_2_MioLDV,dim=c(3.1,3.2)),"Est ICE LDV Stock (million vehicles)")) - - tmp <- mbind(tmp,setNames(dimSums(vm_deltaCap[,,LDV35] * pm_conv_cap_2_MioLDV,dim=c(3.1,3.2)),"Est LDV Sales (million vehicles)")) - tmp <- mbind(tmp,setNames(dimSums(vm_deltaCap[,,"apCarElT"] * pm_conv_cap_2_MioLDV,dim=c(3.1,3.2)),"Est EV LDV Sales (million vehicles)")) - tmp <- mbind(tmp,setNames(dimSums(vm_deltaCap[,,"apCarH2T"] * pm_conv_cap_2_MioLDV,dim=c(3.1,3.2)),"Est H2 LDV Sales (million vehicles)")) - tmp <- mbind(tmp,setNames(dimSums(vm_deltaCap[,,"apCarPeT"] * pm_conv_cap_2_MioLDV,dim=c(3.1,3.2)),"Est ICE LDV Sales (million vehicles)")) - - tmp <- mbind(tmp,setNames(dimSums(vm_cap[,,c("apCarDiT","apcarDiEffT","apcarDiEffH2T")],dim=c(3.1,3.2)), "Services and Products|Transport|non-LDV|Stock|uedit (arbitrary unit)")) - tmp <- mbind(tmp,setNames(dimSums(vm_cap[,,"apCarDiT"] ,dim=c(3.1,3.2)), "Services and Products|Transport|non-LDV|Stock|apCarDiT (arbitrary unit)")) - tmp <- mbind(tmp,setNames(dimSums(vm_cap[,,"apcarDiEffT"] ,dim=c(3.1,3.2)), "Services and Products|Transport|non-LDV|Stock|apcarDiEffT (arbitrary unit)")) - tmp <- mbind(tmp,setNames(dimSums(vm_cap[,,"apcarDiEffH2T"] ,dim=c(3.1,3.2)), "Services and Products|Transport|non-LDV|Stock|apcarDiEffH2T (arbitrary unit)")) - - tmp <- mbind(tmp,setNames(dimSums(vm_deltaCap[,,c("apCarDiT","apcarDiEffT","apcarDiEffH2T")],dim=c(3.1,3.2)), "Services and Products|Transport|non-LDV|Sales|uedit (arbitrary unit)")) - tmp <- mbind(tmp,setNames(dimSums(vm_deltaCap[,,"apCarDiT"] ,dim=c(3.1,3.2)), "Services and Products|Transport|non-LDV|Sales|apCarDiT (arbitrary unit)")) - tmp <- mbind(tmp,setNames(dimSums(vm_deltaCap[,,"apcarDiEffT"] ,dim=c(3.1,3.2)), "Services and Products|Transport|non-LDV|Sales|apcarDiEffT (arbitrary unit)")) - tmp <- mbind(tmp,setNames(dimSums(vm_deltaCap[,,"apcarDiEffH2T"] ,dim=c(3.1,3.2)), "Services and Products|Transport|non-LDV|Sales|apcarDiEffH2T (arbitrary unit)")) - - - ## add global values - tmp <- mbind(tmp,dimSums(tmp,dim=1)) +reportCapitalStock <- function(gdx, regionSubsetList = NULL, + t = c(seq(2005, 2060, 5), seq(2070, 2110, 10), 2130, 2150), + gdx_ref = gdx_ref) { + + module2realisation <- readGDX(gdx, "module2realisation", react = "silent") + tran_mod <- module2realisation[module2realisation$modules == "transport", 2] + + pm_conv_cap_2_MioLDV <- 650 # The world has ~715million cars in 2005 (IEA TECO2) + + # read sets + teall2rlf <- readGDX(gdx, name = c("te2rlf", "teall2rlf"), format = "first_found") + teue2rlf <- readGDX(gdx, name = c("teue2rlf", "tees2rlf"), format = "first_found") + + # read variables + vm_cap <- readGDX(gdx, name = c("vm_cap"), field = "l", format = "first_found") + vm_deltaCap <- readGDX(gdx, name = c("vm_deltaCap"), field = "l", format = "first_found") + + v_investcost <- readGDX(gdx, name = c("vm_costTeCapital", "v_costTeCapital", "v_investcost"), field = "l", format = "first_found") + vm_cesIO <- readGDX(gdx, name = "vm_cesIO", field = "l") + + # read parameters + ppfKap_Ind <- readGDX(gdx, name = "ppfkap_industry_dyn37", react = "silent") + steel_process_based <- "steel" %in% readGDX(gdx, "secInd37Prc", react = "silent") + + # calculate maximal temporal resolution + y <- Reduce(intersect, list(getYears(vm_cap), getYears(v_investcost))) + vm_cap <- vm_cap[, y, ] + vm_deltaCap <- vm_deltaCap[, y, ] + v_investcost <- v_investcost[, y, ] + + if (!is.null(gdx_ref)) { + cm_startyear <- as.integer(readGDX(gdx, name = "cm_startyear", format = "simplest")) + vm_deltaCapRef <- readGDX(gdx_ref, name = c("vm_deltaCap"), field = "l", format = "first_found")[, y, ] + vm_deltaCap <- modifyInvestmentVariables(vm_deltaCap, vm_deltaCapRef, cm_startyear) + } else { + vm_deltaCap <- modifyInvestmentVariables(vm_deltaCap) + } + + tmp <- NULL + + # ---- report transport capital stocks ---- + if (tran_mod == "complex") { + LDV35 <- readGDX(gdx, name = c("LDV35"), format = "first_found") + tmp <- mbind(tmp, setNames(dimSums((vm_cap * v_investcost)[teue2rlf], + dim = c(3.1, 3.2)) * 1000, "Est Capital Stock|ESM|Transp vehic (billion US$2017)")) + tmp <- mbind(tmp, setNames(dimSums((vm_cap * v_investcost)[teall2rlf][, , LDV35], + dim = c(3.1, 3.2)) * 1000, "Est Capital Stock|ESM|Pet/EV LDV (billion US$2017)")) + + tmp <- mbind(tmp, + setNames( + dimSums(mbind(vm_cap * v_investcost), + dim = 3) * 1000, + "Estimated Capital Stock|ESM (billion US$2017)" + ) + ) + + # prepare variables + vm_cap <- vm_cap[teall2rlf] + vm_deltaCap <- vm_deltaCap[teall2rlf] + + # build reporting + tmp <- mbind(tmp, setNames(dimSums(vm_cap[, , LDV35] * pm_conv_cap_2_MioLDV, dim = c(3.1, 3.2)), "Est LDV Stock (million vehicles)")) + tmp <- mbind(tmp, setNames(dimSums(vm_cap[, , "apCarElT"] * pm_conv_cap_2_MioLDV, dim = c(3.1, 3.2)), "Est EV LDV Stock (million vehicles)")) + tmp <- mbind(tmp, setNames(dimSums(vm_cap[, , "apCarH2T"] * pm_conv_cap_2_MioLDV, dim = c(3.1, 3.2)), "Est H2 LDV Stock (million vehicles)")) + tmp <- mbind(tmp, setNames(dimSums(vm_cap[, , "apCarPeT"] * pm_conv_cap_2_MioLDV, dim = c(3.1, 3.2)), "Est ICE LDV Stock (million vehicles)")) + + tmp <- mbind(tmp, setNames(dimSums(vm_cap[, , c("apCarDiT", "apcarDiEffT", "apcarDiEffH2T")], dim = c(3.1, 3.2)), "Services and Products|Transport|non-LDV|Stock|uedit (arbitrary unit)")) + tmp <- mbind(tmp, setNames(dimSums(vm_cap[, , "apCarDiT"], dim = c(3.1, 3.2)), "Services and Products|Transport|non-LDV|Stock|apCarDiT (arbitrary unit)")) + tmp <- mbind(tmp, setNames(dimSums(vm_cap[, , "apcarDiEffT"], dim = c(3.1, 3.2)), "Services and Products|Transport|non-LDV|Stock|apcarDiEffT (arbitrary unit)")) + tmp <- mbind(tmp, setNames(dimSums(vm_cap[, , "apcarDiEffH2T"], dim = c(3.1, 3.2)), "Services and Products|Transport|non-LDV|Stock|apcarDiEffH2T (arbitrary unit)")) + + tmp2 <- NULL + + tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , LDV35] * pm_conv_cap_2_MioLDV, dim = c(3.1, 3.2)), "Est LDV Sales (million vehicles)")) + tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , "apCarElT"] * pm_conv_cap_2_MioLDV, dim = c(3.1, 3.2)), "Est EV LDV Sales (million vehicles)")) + tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , "apCarH2T"] * pm_conv_cap_2_MioLDV, dim = c(3.1, 3.2)), "Est H2 LDV Sales (million vehicles)")) + tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , "apCarPeT"] * pm_conv_cap_2_MioLDV, dim = c(3.1, 3.2)), "Est ICE LDV Sales (million vehicles)")) + + tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , c("apCarDiT", "apcarDiEffT", "apcarDiEffH2T")], dim = c(3.1, 3.2)), "Services and Products|Transport|non-LDV|Sales|uedit (arbitrary unit)")) + tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , "apCarDiT"], dim = c(3.1, 3.2)), "Services and Products|Transport|non-LDV|Sales|apCarDiT (arbitrary unit)")) + tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , "apcarDiEffT"], dim = c(3.1, 3.2)), "Services and Products|Transport|non-LDV|Sales|apcarDiEffT (arbitrary unit)")) + tmp2 <- mbind(tmp2, setNames(dimSums(vm_deltaCap[, , "apcarDiEffH2T"], dim = c(3.1, 3.2)), "Services and Products|Transport|non-LDV|Sales|apcarDiEffH2T (arbitrary unit)")) + + tmp <- mbind(tmp, tmp2) + + ## add global values + tmp <- mbind(tmp, dimSums(tmp, dim = 1)) + } + + # ---- report industry energy efficiency capital stocks ---- + if (!is.null(ppfKap_Ind) && 0 < length(ppfKap_Ind)) { + mixer <- tribble( + ~pf, ~name, + "kap_cement", "Cement", + "kap_chemicals", "Chemicals", + "kap_otherInd", "other") + + if (!steel_process_based) { + mixer <- bind_rows(mixer, tribble( + ~pf, ~name, + "kap_steel_primary", "Primary Steel", + "kap_steel_secondary", "Secondary Steel")) } - # ---- report industry energy efficiency capital stocks ---- - if (!is.null(ppfKap_Ind) & 0 < length(ppfKap_Ind)) { - mixer <- tribble( - ~pf, ~name, - 'kap_cement', 'Cement', - 'kap_chemicals', 'Chemicals', - 'kap_otherInd', 'other') + if (0 != length(setdiff(ppfKap_Ind, mixer$pf))) { + warning(paste("Unknown ppfKap_industry_dyn37 entity.", + "Adjust remind2::reportCapitalStock()")) + } - if (!steel_process_based) { - mixer <- bind_rows(mixer, tribble( - ~pf, ~name, - 'kap_steel_primary', 'Primary Steel', - 'kap_steel_secondary', 'Secondary Steel')) - } - - if (0 != length(setdiff(ppfKap_Ind, mixer$pf))) { - warning(paste('Unknown ppfKap_industry_dyn37 entity.', - 'Adjust remind2::reportCapitalStock()')) - } - - if (0 != length(setdiff(mixer$pf, ppfKap_Ind))) { - warning(paste('Missing ppfKap_industry_dyn37 entity.', - 'Adjust remind2::reportCapitalStock()')) - } - - eek_Ind <- setNames(vm_cesIO[,y,ppfKap_Ind], - paste0('Capital|Energy Efficiency|Industry|', - mixer[mixer$pf %in% ppfKap_Ind,][['name']], - ' (billion US$2017)')) - # add industry EEK and global totals - tmp <- mbind(tmp, mbind(eek_Ind, dimSums(eek_Ind, dim = 1))) + if (0 != length(setdiff(mixer$pf, ppfKap_Ind))) { + warning(paste("Missing ppfKap_industry_dyn37 entity.", + "Adjust remind2::reportCapitalStock()")) } - # ---- add region aggregates ---- - if (!is.null(regionSubsetList)) - tmp <- mbind(tmp, calc_regionSubset_sums(tmp, regionSubsetList)) + eek_Ind <- setNames(vm_cesIO[, y, ppfKap_Ind], + paste0("Capital|Energy Efficiency|Industry|", + mixer[mixer$pf %in% ppfKap_Ind, ][["name"]], + " (billion US$2017)")) + # add industry EEK and global totals + tmp <- mbind(tmp, mbind(eek_Ind, dimSums(eek_Ind, dim = 1))) + } + + # ---- add region aggregates ---- + if (!is.null(regionSubsetList)) + tmp <- mbind(tmp, calc_regionSubset_sums(tmp, regionSubsetList)) + + getSets(tmp)[3] <- "variable" - getSets(tmp)[3] <- "variable" - return(tmp) + return(tmp) } diff --git a/R/reportEnergyInvestment.R b/R/reportEnergyInvestment.R index 5cff9225..6b6632ee 100644 --- a/R/reportEnergyInvestment.R +++ b/R/reportEnergyInvestment.R @@ -10,6 +10,9 @@ #' be created. #' @param t temporal resolution of the reporting, default: #' t=c(seq(2005,2060,5),seq(2070,2110,10),2130,2150) +#' @param gdx_ref a GDX object as created by readGDX, or the path to a gdx of the reference run. +#' It is used to guarantee consistency before 'cm_startyear' for investment variables +#' using time averaging. #' #' @return MAgPIE object - contains the price variables #' @author Anastasis Giannousaki @@ -20,11 +23,12 @@ #' } #' #' @export -#' @importFrom magclass mbind +#' @importFrom magclass mbind is.magpie #' @importFrom gdx readGDX reportEnergyInvestment <- function(gdx, regionSubsetList = NULL, - t = c(seq(2005, 2060, 5), seq(2070, 2110, 10), 2130, 2150)) { + t = c(seq(2005, 2060, 5), seq(2070, 2110, 10), 2130, 2150), + gdx_ref = NULL) { # read sets adjte <- readGDX(gdx, name = c("teAdj", "adjte"), format = "first_found") petyf <- readGDX(gdx, c("peFos", "petyf"), format = "first_found") @@ -52,9 +56,13 @@ reportEnergyInvestment <- function(gdx, regionSubsetList = NULL, grid <- readGDX(gdx, name = c("teGrid", "grid"), format = "first_found") pebio <- readGDX(gdx, c("peBio", "pebio"), format = "first_found") ttot <- readGDX(gdx, name = "ttot", format = "first_found") + # read variables + + # investment cost per technology v_directteinv <- readGDX(gdx, name = c("v_costInvTeDir", "vm_costInvTeDir", "v_directteinv"), field = "l", format = "first_found") + # adjustment cost per technology v_adjustteinv <- readGDX(gdx, name = c("v_costInvTeAdj", "vm_costInvTeAdj", "v_adjustteinv"), field = "l", format = "first_found") @@ -62,11 +70,26 @@ reportEnergyInvestment <- function(gdx, regionSubsetList = NULL, pm_data <- readGDX(gdx, c("pm_data"), format = "first_found") # data preparation + ttot <- as.numeric(as.vector(readGDX(gdx, "ttot", format = "first_found"))) - v_directteinv <- v_directteinv[, ttot, ] - v_adjustteinv <- v_adjustteinv[, ttot, ] - costRatioTdelt2Tdels <- pm_data[, , "inco0.tdelt"] / pm_data[, , "inco0.tdels"] + # apply 'modifyInvestmentVariables' to shift from the model-internal time coverage (deltacap and investment + # variables for step t represent the average of the years from t-4years to t) to the general convention for + # the reporting template (all variables represent the average of the years from t-2.5years to t+2.5years) + if (!is.null(gdx_ref)) { + v_directteinv_ref <- readGDX(gdx_ref, name = c("v_costInvTeDir", "vm_costInvTeDir", "v_directteinv"), + field = "l", format = "first_found") + v_adjustteinv_ref <- readGDX(gdx_ref, name = c("v_costInvTeAdj", "vm_costInvTeAdj", "v_adjustteinv"), + field = "l", format = "first_found") + cm_startyear <- as.integer(readGDX(gdx, name = "cm_startyear", format = "simplest")) + v_directteinv <- modifyInvestmentVariables(v_directteinv[, ttot, ], v_directteinv_ref[, ttot, ], cm_startyear) + v_adjustteinv <- modifyInvestmentVariables(v_adjustteinv[, ttot, ], v_adjustteinv_ref[, ttot, ], cm_startyear) + } else { + v_directteinv <- modifyInvestmentVariables(v_directteinv[, ttot, ]) + v_adjustteinv <- modifyInvestmentVariables(v_adjustteinv[, ttot, ]) + } + + costRatioTdelt2Tdels <- pm_data[, , "inco0.tdelt"] / pm_data[, , "inco0.tdels"] ####### internal function for reporting ########### ## "ie" stands for input energy, "oe" for output energy @@ -224,7 +247,6 @@ reportEnergyInvestment <- function(gdx, regionSubsetList = NULL, - tmp[, , "Energy Investments|CO2 Trans&Stor (billion US$2017/yr)"]), "Energy Investments|Other (billion US$2017/yr)")) - # add global values tmp <- mbind(tmp, dimSums(tmp, dim = 1)) # add other region aggregations @@ -232,5 +254,6 @@ reportEnergyInvestment <- function(gdx, regionSubsetList = NULL, tmp <- mbind(tmp, calc_regionSubset_sums(tmp, regionSubsetList)) getSets(tmp)[3] <- "variable" + return(tmp) } diff --git a/R/reportLCOE.R b/R/reportLCOE.R index 290d8f9d..e4fcad90 100644 --- a/R/reportLCOE.R +++ b/R/reportLCOE.R @@ -138,9 +138,16 @@ reportLCOE <- function(gdx, output.type = "both") { pm_SEPrice <- readGDX(gdx, "pm_SEPrice", restore_zeros = FALSE) ## variables - vm_costInvTeDir <- readGDX(gdx, name = c("vm_costInvTeDir", "v_costInvTeDir", "v_directteinv"), field = "l", format = "first_found")[, ttot, ] ## Total direct Investment Cost in Timestep - vm_costInvTeAdj <- readGDX(gdx, name = c("vm_costInvTeAdj", "v_costInvTeAdj"), field = "l", format = "first_found")[, ttot, ] ## total adjustment cost in period - vm_deltaCap <- readGDX(gdx, name = c("vm_deltaCap"), field = "l", format = "first_found")[, ttot, ] + + ## Total direct Investment Cost in Timestep + vm_costInvTeDir <- readGDX(gdx, name = c("vm_costInvTeDir", "v_costInvTeDir", "v_directteinv"), field = "l", format = "first_found")[, ttot, ] + + ## total adjustment cost in period + vm_costInvTeAdj <- readGDX(gdx, name = c("vm_costInvTeAdj", "v_costInvTeAdj"), field = "l", format = "first_found")[, ttot, ] + + # capacity additions per year + vm_deltaCap <- readGDX(gdx, name = c("vm_deltaCap"), field = "l", format = "first_found")[, ttot, ] + vm_demPe <- readGDX(gdx, name = c("vm_demPe", "v_pedem"), field = "l", restore_zeros = FALSE, format = "first_found") v_investcost <- readGDX(gdx, name = c("vm_costTeCapital", "v_costTeCapital", "v_investcost"), field = "l", format = "first_found")[, ttot, ] vm_cap <- readGDX(gdx, name = c("vm_cap"), field = "l", format = "first_found") @@ -1599,13 +1606,12 @@ reportLCOE <- function(gdx, output.type = "both") { LCOE.out.inclGlobal[getRegions(LCOE.out), , ] <- LCOE.out LCOE.out.inclGlobal["GLO", , ] <- dimSums(LCOE.out, dim = 1) / length(getRegions(LCOE.out)) - - - if (output.type %in% c("marginal detail")) { - return(df.LCOE) + if (output.type == "marginal detail") { + out <- df.LCOE } else { - return(LCOE.out.inclGlobal) + out <- LCOE.out.inclGlobal } - return(LCOE.out) + return(out) + } diff --git a/R/reportPrices.R b/R/reportPrices.R index 7f5628af..0b948544 100644 --- a/R/reportPrices.R +++ b/R/reportPrices.R @@ -18,23 +18,24 @@ #' @author Alois Dirnaichner, Felix Schreyer, David Klein, Renato Rodrigues, Falk Benke #' @seealso \code{\link{convGDX2MIF}} #' @examples -#' -#' \dontrun{reportPrices(gdx)} +#' \dontrun{ +#' reportPrices(gdx) +#' } #' #' @importFrom dplyr %>% case_when distinct filter inner_join tibble left_join rename #' @importFrom gdx readGDX -#' @importFrom magclass mbind getYears getRegions setNames dimExists new.magpie lowpass complete_magpie getItems<- getNames unitsplit +#' @importFrom magclass mbind getYears getRegions setNames dimExists new.magpie +#' lowpass complete_magpie getItems<- getNames unitsplit #' @importFrom quitte df.2.named.vector getColValues #' @importFrom readr read_csv #' @importFrom madrat toolAggregate #' - #' @export -reportPrices <- function(gdx, output=NULL, regionSubsetList=NULL, - t=c(seq(2005,2060,5),seq(2070,2110,10),2130,2150),gdx_ref=NULL) { - +reportPrices <- function(gdx, output = NULL, regionSubsetList = NULL, + t = c(seq(2005, 2060, 5), seq(2070, 2110, 10), 2130, 2150), + gdx_ref = NULL) { ## bind to output object - if(is.null(output)){ + if (is.null(output)) { message("reportPrices executes reportPE ", appendLF = FALSE) output <- reportPE(gdx, regionSubsetList = regionSubsetList, t = t) message("- reportSE ", appendLF = FALSE) @@ -44,82 +45,66 @@ reportPrices <- function(gdx, output=NULL, regionSubsetList=NULL, message("- reportEmi ", appendLF = FALSE) output <- mbind(output, reportEmi(gdx, output = output, regionSubsetList = regionSubsetList, t = t)) message("- reportExtraction ", appendLF = FALSE) - output <- mbind(output, reportExtraction(gdx,regionSubsetList = regionSubsetList, t = t)) + output <- mbind(output, reportExtraction(gdx, regionSubsetList = regionSubsetList, t = t)) message("- reportMacroEconomy") - output <- mbind(output, reportMacroEconomy(gdx,regionSubsetList = regionSubsetList, t = t)[, getYears(output), ]) + output <- mbind(output, reportMacroEconomy(gdx, regionSubsetList = regionSubsetList, t = t)[, getYears(output), ]) } output[is.na(output)] <- 0 # substitute na by 0 output <- deletePlus(output) # delete "+" and "++" from variable names # get rid of GLO from output - all_regi_wo_GLO <- getItems(output, dim = "all_regi")[! getItems(output, dim = "all_regi") %in% c("GLO")] + all_regi_wo_GLO <- getItems(output, dim = "all_regi")[!getItems(output, dim = "all_regi") %in% c("GLO")] output_wo_GLO <- output[all_regi_wo_GLO, , ] - ####### get realisations ######### - realisation <- readGDX(gdx, "module2realisation") - ####### conversion factors ########## - s_GWP_CH4 <- readGDX(gdx, c("sm_gwpCH4","s_gwpCH4","s_GWP_CH4"), format="first_found", react = "silent") - s_GWP_N2O <- readGDX(gdx, c("s_gwpN2O","s_GWP_N2O"), format="first_found", react = "silent") + s_GWP_CH4 <- readGDX(gdx, c("sm_gwpCH4", "s_gwpCH4", "s_GWP_CH4"), format = "first_found", react = "silent") + s_GWP_N2O <- readGDX(gdx, c("s_gwpN2O", "s_GWP_N2O"), format = "first_found", react = "silent") s_twa2mwh <- readGDX(gdx, "sm_TWa_2_MWh", format = "first_found", reacht = "silent") - tdptwyr2dpgj <- 31.71 #TerraDollar per TWyear to Dollar per GJ - p80_subset <- c("perm", "good", "peur", "peoil", "pegas", "pecoal", "pebiolc") #TODO: read in from gdx as sets trade + tdptwyr2dpgj <- 31.71 # TerraDollar per TWyear to Dollar per GJ + p80_subset <- c("perm", "good", "peur", "peoil", "pegas", "pecoal", "pebiolc") # TODO: read in from gdx as sets trade ####### read in needed data ######### #---- Functions - find_real_module <- function(module_set, module_name){ - return(module_set[module_set$modules == module_name, 2]) - } - - indu_mod = find_real_module(realisation,"industry") ## parameter - shift_p <- readGDX(gdx,name="p30_pebiolc_pricshift",format="first_found")[, t,] - mult_p <- readGDX(gdx,name="p30_pebiolc_pricmult",format="first_found")[, t,] - pric_mag <- readGDX(gdx,name="p30_pebiolc_pricemag",format="first_found")[, t,] - pric_emu_pre <- readGDX(gdx,name="p30_pebiolc_price_emu_preloop",format="first_found")[, t,] - pric_emu_pre_shifted <- readGDX(gdx,name="p30_pebiolc_price_emu_preloop_shifted",format="first_found")[, t,] - bio_tax_factor <- readGDX(gdx,name="p21_tau_bioenergy_tax",format="first_found",react="silent")[, t,] - if (is.null(bio_tax_factor)) bio_tax_factor <- readGDX(gdx,name="v21_tau_bio",field="l",format="first_found")[, t,] - pm_pvp <- readGDX(gdx,name=c("pm_pvp","p80_pvp"),format="first_found")[, t, p80_subset] - pm_dataemi <- readGDX(gdx,name=c("pm_emifac","pm_dataemi"),format="first_found",restore_zeros=FALSE)[,t, c("pegas.seel.ngt.co2","pecoal.seel.pc.co2")] - pm_taxCO2eq <- readGDX(gdx,name=c("pm_taxCO2eq","pm_tau_CO2_tax"),format="first_found")[, t,] - pm_taxCO2eqSum <- readGDX(gdx,name="pm_taxCO2eqSum",format="first_found")[, t,] - pm_taxCO2eqSCC <- readGDX(gdx,name='pm_taxCO2eqSCC',format="first_found")[, t,] - p21_CO2TaxSectorMarkup <- readGDX(gdx,name=c('p21_CO2TaxSectorMarkup','p21_CO2_tax_sector_markup'),format="first_found",react="silent") - if (dimExists("ttot", p21_CO2TaxSectorMarkup)) p21_CO2TaxSectorMarkup <- p21_CO2TaxSectorMarkup[, t,] - pm_taxemiMkt <- readGDX(gdx,name="pm_taxemiMkt",format="first_found",react="silent")[, t,] - p47_taxCO2eq_AggFE <- readGDX(gdx,name="p47_taxCO2eq_AggFE",format="first_found",react="silent")[, t,] - p47_taxCO2eq_SectorAggFE <- readGDX(gdx,name="p47_taxCO2eq_SectorAggFE",format="first_found",react="silent")[, t,] + shift_p <- readGDX(gdx, name = "p30_pebiolc_pricshift", format = "first_found")[, t, ] + mult_p <- readGDX(gdx, name = "p30_pebiolc_pricmult", format = "first_found")[, t, ] + pric_mag <- readGDX(gdx, name = "p30_pebiolc_pricemag", format = "first_found")[, t, ] + pric_emu_pre <- readGDX(gdx, name = "p30_pebiolc_price_emu_preloop", format = "first_found")[, t, ] + pric_emu_pre_shifted <- readGDX(gdx, name = "p30_pebiolc_price_emu_preloop_shifted", format = "first_found")[, t, ] + bio_tax_factor <- readGDX(gdx, name = "p21_tau_bioenergy_tax", format = "first_found", react = "silent")[, t, ] + if (is.null(bio_tax_factor)) bio_tax_factor <- readGDX(gdx, name = "v21_tau_bio", field = "l", format = "first_found")[, t, ] + pm_pvp <- readGDX(gdx, name = c("pm_pvp", "p80_pvp"), format = "first_found")[, t, p80_subset] + pm_taxCO2eq <- readGDX(gdx, name = c("pm_taxCO2eq", "pm_tau_CO2_tax"), format = "first_found")[, t, ] + pm_taxCO2eqSum <- readGDX(gdx, name = "pm_taxCO2eqSum", format = "first_found")[, t, ] + pm_taxCO2eqSCC <- readGDX(gdx, name = "pm_taxCO2eqSCC", format = "first_found")[, t, ] + p21_CO2TaxSectorMarkup <- readGDX(gdx, name = c("p21_CO2TaxSectorMarkup", "p21_CO2_tax_sector_markup"), format = "first_found", react = "silent") + if (dimExists("ttot", p21_CO2TaxSectorMarkup)) p21_CO2TaxSectorMarkup <- p21_CO2TaxSectorMarkup[, t, ] + pm_taxemiMkt <- readGDX(gdx, name = "pm_taxemiMkt", format = "first_found", react = "silent")[, t, ] + p47_taxCO2eq_AggFE <- readGDX(gdx, name = "p47_taxCO2eq_AggFE", format = "first_found", react = "silent")[, t, ] + p47_taxCO2eq_SectorAggFE <- readGDX(gdx, name = "p47_taxCO2eq_SectorAggFE", format = "first_found", react = "silent")[, t, ] ## variables - pric_emu <- readGDX(gdx,name="vm_pebiolc_price",field="l",format="first_found")[, t,] + pric_emu <- readGDX(gdx, name = "vm_pebiolc_price", field = "l", format = "first_found")[, t, ] ## equations - budget.m <- readGDX(gdx,name='qm_budget',types = "equations",field = "m",format = "first_found")[, t,] # Alternative: calcPrice - balcapture.m <- readGDX(gdx,name=c("q_balcapture", "q12_balcapture"), field = "m", restore_zeros = F)[, t,] - - esm2macro.m <- readGDX(gdx,name='q35_esm2macro',types="equations",field="m",format="first_found", react = "silent")[, t,] - - cm_emiscen <- readGDX(gdx,name='cm_emiscen',format="first_found") - cm_startyear <- as.integer(readGDX(gdx,name='cm_startyear',format='simplest')) - - q37_limit_secondary_steel_share.m <- readGDX( - gdx, name = 'q37_limit_secondary_steel_share', types = 'equation', - field = 'm', react = 'silent')[, t,] + budget.m <- readGDX(gdx, name = "qm_budget", types = "equations", field = "m", format = "first_found")[, t, ] # Alternative: calcPrice + balcapture.m <- readGDX(gdx, name = c("q_balcapture", "q12_balcapture"), field = "m", restore_zeros = FALSE)[, t, ] + esm2macro.m <- readGDX(gdx, name = "q35_esm2macro", types = "equations", field = "m", format = "first_found", react = "silent")[, t, ] + cm_startyear <- as.integer(readGDX(gdx, name = "cm_startyear", format = "simplest")) ##################################### - #choose the CES entries names for transport + # choose the CES entries names for transport if (!is.null(esm2macro.m)) { - name_trsp=c("fepet","ueLDVt","fedie","ueHDVt","feelt","ueelTt") - name_trsp=name_trsp[name_trsp%in%getNames(esm2macro.m)] + name_trsp <- c("fepet", "ueLDVt", "fedie", "ueHDVt", "feelt", "ueelTt") + name_trsp <- name_trsp[name_trsp %in% getNames(esm2macro.m)] } ##################################### ####### pm_pvp = EPS results in fantasy carbon prices - for(yr in getYears(pm_pvp)){ - if(pm_pvp[,yr,"good"] == 0){ - pm_pvp[,yr,"good"] = 0.0001; + for (yr in getYears(pm_pvp)) { + if (pm_pvp[, yr, "good"] == 0) { + pm_pvp[, yr, "good"] <- 0.0001 } } @@ -131,10 +116,10 @@ reportPrices <- function(gdx, output=NULL, regionSubsetList=NULL, pm_SEPrice <- readGDX(gdx, "pm_SEPrice") pm_PEPrice <- readGDX(gdx, c("p_PEPrice", "pm_PEPrice"), format = "first_found") - vm_demFeSector <- readGDX(gdx, "vm_demFeSector", field = "l", restore_zeros = FALSE)[,t,] - prodSe <- readGDX(gdx, "vm_prodSe", field = "l", restore_zeros = FALSE)[,t,] - try(seAgg <- readGDX(gdx, name="seAgg", type="set")) - try(seAgg2se <- readGDX(gdx, name="seAgg2se", type="set")) + vm_demFeSector <- readGDX(gdx, "vm_demFeSector", field = "l", restore_zeros = FALSE)[, t, ] + prodSe <- readGDX(gdx, "vm_prodSe", field = "l", restore_zeros = FALSE)[, t, ] + try(seAgg <- readGDX(gdx, name = "seAgg", type = "set")) + try(seAgg2se <- readGDX(gdx, name = "seAgg2se", type = "set")) # subset price parameters to those entries used in the model pe2se <- readGDX(gdx, "pe2se") @@ -143,33 +128,33 @@ reportPrices <- function(gdx, output=NULL, regionSubsetList=NULL, entyFe2Sector <- readGDX(gdx, "entyFe2Sector") - sector <- emi_sectors <- emiMkt <- all_emiMkt <- NULL + sector <- emi_sectors <- all_emiMkt <- NULL fe.entries <- entyFe2Sector %>% - left_join(sector2emiMkt, by = "emi_sectors", relationship = "many-to-many") %>% - rename( sector = emi_sectors, emiMkt = all_emiMkt) %>% - filter( sector != "CDR") + left_join(sector2emiMkt, by = "emi_sectors", relationship = "many-to-many") %>% + rename(sector = emi_sectors, emiMkt = all_emiMkt) %>% + filter(sector != "CDR") - fe.entries.dot <- paste(fe.entries[,1],fe.entries[,2], fe.entries[,3], sep = ".") + fe.entries.dot <- paste(fe.entries[, 1], fe.entries[, 2], fe.entries[, 3], sep = ".") ttot <- readGDX(gdx, "ttot") - YearsFrom2005 <- paste0("y",ttot[ttot >= 2005]) + YearsFrom2005 <- paste0("y", ttot[ttot >= 2005]) - pm_PEPrice <- pm_PEPrice[,YearsFrom2005,unique(pe2se$all_enty)] - pm_SEPrice <- pm_SEPrice[,YearsFrom2005,unique(se2fe$all_enty)] - pm_FEPrice <- pm_FEPrice[,YearsFrom2005,fe.entries.dot] + pm_PEPrice <- pm_PEPrice[, YearsFrom2005, unique(pe2se$all_enty)] + pm_SEPrice <- pm_SEPrice[, YearsFrom2005, unique(se2fe$all_enty)] + pm_FEPrice <- pm_FEPrice[, YearsFrom2005, fe.entries.dot] ## weights for market aggregation of prices: FE share of market - p_weights_FEprice_mkt <- dimSums(vm_demFeSector, dim=3.1, na.rm = T) / dimSums(vm_demFeSector, dim=c(3.1,3.4), na.rm = T) + p_weights_FEprice_mkt <- dimSums(vm_demFeSector, dim = 3.1, na.rm = TRUE) / dimSums(vm_demFeSector, dim = c(3.1, 3.4), na.rm = TRUE) p_weights_FEprice_mkt[is.na(p_weights_FEprice_mkt)] <- 0 ## adjust to pm_FEprice dimensions - p_weights_FEprice_mkt <- p_weights_FEprice_mkt[,getYears(pm_FEPrice),getNames(pm_FEPrice)] + p_weights_FEprice_mkt <- p_weights_FEprice_mkt[, getYears(pm_FEPrice), getNames(pm_FEPrice)] ## weights for fepet/fedie to liquids aggregation of prices: FE fepet/fedie share of aggregated markets - p_weights_FEprice_diepet <- dimSums(mselect(vm_demFeSector, all_enty1=c("fedie","fepet"), emi_sectors="trans"), dim=c(3.1,3.4), na.rm = T) / - dimSums(mselect(vm_demFeSector, all_enty1=c("fedie","fepet"), emi_sectors="trans"), dim=c(3.1,3.2,3.4), na.rm = T) + p_weights_FEprice_diepet <- dimSums(mselect(vm_demFeSector, all_enty1 = c("fedie", "fepet"), emi_sectors = "trans"), dim = c(3.1, 3.4), na.rm = TRUE) / + dimSums(mselect(vm_demFeSector, all_enty1 = c("fedie", "fepet"), emi_sectors = "trans"), dim = c(3.1, 3.2, 3.4), na.rm = TRUE) p_weights_FEprice_diepet[is.na(p_weights_FEprice_diepet)] <- 0 - p_weights_FEprice_diepet <- p_weights_FEprice_diepet[,getYears(pm_FEPrice),] + p_weights_FEprice_diepet <- p_weights_FEprice_diepet[, getYears(pm_FEPrice), ] out <- NULL @@ -177,99 +162,99 @@ reportPrices <- function(gdx, output=NULL, regionSubsetList=NULL, ## FE Transport Prices out <- mbind( out, - setNames( dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1="feelt", emi_sectors = "trans"), dim=3.3, na.rm = T)*tdptwyr2dpgj, + setNames(dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1 = "feelt", emi_sectors = "trans"), dim = 3.3, na.rm = TRUE) * tdptwyr2dpgj, "Price|Final Energy|Transport|Electricity (US$2017/GJ)"), ## in case of transport liquids: calculate weighted average of markets first, then calculate weighted average of fepet/fedie - setNames( dimSums( p_weights_FEprice_diepet * dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1=c("fepet","fedie"), emi_sectors = "trans"), dim=3.3, na.rm = T), dim=3.1, na.rm = T)*tdptwyr2dpgj, + setNames(dimSums(p_weights_FEprice_diepet * dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1 = c("fepet", "fedie"), emi_sectors = "trans"), dim = 3.3, na.rm = TRUE), dim = 3.1, na.rm = TRUE) * tdptwyr2dpgj, "Price|Final Energy|Transport|Liquids (US$2017/GJ)"), - setNames( dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1="feh2t", emi_sectors = "trans"), dim=3.3, na.rm = T)*tdptwyr2dpgj, + setNames(dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1 = "feh2t", emi_sectors = "trans"), dim = 3.3, na.rm = TRUE) * tdptwyr2dpgj, "Price|Final Energy|Transport|Hydrogen (US$2017/GJ)") ) - if (module2realisation["transport",2] == "edge_esm") { + if (module2realisation["transport", 2] == "edge_esm") { out <- mbind( out, - setNames( dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1="fegat", emi_sectors = "trans"), dim=3.3, na.rm = T)*tdptwyr2dpgj, + setNames(dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1 = "fegat", emi_sectors = "trans"), dim = 3.3, na.rm = TRUE) * tdptwyr2dpgj, "Price|Final Energy|Transport|Gases (US$2017/GJ)")) } out <- mbind( out, - setNames( dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1="fedie", emi_sectors = "trans"), dim=3.3, na.rm = T)*tdptwyr2dpgj, + setNames(dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1 = "fedie", emi_sectors = "trans"), dim = 3.3, na.rm = TRUE) * tdptwyr2dpgj, "Price|Final Energy|Transport|Liquids|HDV (US$2017/GJ)"), - setNames( dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1="fepet", emi_sectors = "trans"), dim=3.3, na.rm = T)*tdptwyr2dpgj, + setNames(dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1 = "fepet", emi_sectors = "trans"), dim = 3.3, na.rm = TRUE) * tdptwyr2dpgj, "Price|Final Energy|Transport|Liquids|LDV (US$2017/GJ)")) ## FE Buildings Prices out <- mbind(out, - setNames( dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1="feels", emi_sectors = "build"), dim=3.3, na.rm = T)*tdptwyr2dpgj, + setNames(dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1 = "feels", emi_sectors = "build"), dim = 3.3, na.rm = TRUE) * tdptwyr2dpgj, "Price|Final Energy|Buildings|Electricity (US$2017/GJ)"), - setNames( dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1="fegas", emi_sectors = "build"), dim=3.3, na.rm = T)*tdptwyr2dpgj, + setNames(dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1 = "fegas", emi_sectors = "build"), dim = 3.3, na.rm = TRUE) * tdptwyr2dpgj, "Price|Final Energy|Buildings|Gases (US$2017/GJ)"), - setNames( dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1="feh2s", emi_sectors = "build"), dim=3.3, na.rm = T)*tdptwyr2dpgj, + setNames(dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1 = "feh2s", emi_sectors = "build"), dim = 3.3, na.rm = TRUE) * tdptwyr2dpgj, "Price|Final Energy|Buildings|Hydrogen (US$2017/GJ)"), - setNames( dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1="fehos", emi_sectors = "build"), dim=3.3, na.rm = T)*tdptwyr2dpgj, + setNames(dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1 = "fehos", emi_sectors = "build"), dim = 3.3, na.rm = TRUE) * tdptwyr2dpgj, "Price|Final Energy|Buildings|Liquids (US$2017/GJ)"), - setNames( dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1="fehes", emi_sectors = "build"), dim=3.3, na.rm = T)*tdptwyr2dpgj, + setNames(dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1 = "fehes", emi_sectors = "build"), dim = 3.3, na.rm = TRUE) * tdptwyr2dpgj, "Price|Final Energy|Buildings|Heat (US$2017/GJ)"), - setNames( dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1="fesos", emi_sectors = "build"), dim=3.3, na.rm = T)*tdptwyr2dpgj, + setNames(dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1 = "fesos", emi_sectors = "build"), dim = 3.3, na.rm = TRUE) * tdptwyr2dpgj, "Price|Final Energy|Buildings|Solids (US$2017/GJ)") - ) + ) ## FE Industry Prices out <- mbind(out, - setNames( dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1="feels", emi_sectors = "indst"), dim=3.3, na.rm = T)*tdptwyr2dpgj, + setNames(dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1 = "feels", emi_sectors = "indst"), dim = 3.3, na.rm = TRUE) * tdptwyr2dpgj, "Price|Final Energy|Industry|Electricity (US$2017/GJ)"), - setNames( dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1="fegas", emi_sectors = "indst"), dim=3.3, na.rm = T)*tdptwyr2dpgj, + setNames(dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1 = "fegas", emi_sectors = "indst"), dim = 3.3, na.rm = TRUE) * tdptwyr2dpgj, "Price|Final Energy|Industry|Gases (US$2017/GJ)"), - setNames( dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1="feh2s", emi_sectors = "indst"), dim=3.3, na.rm = T)*tdptwyr2dpgj, + setNames(dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1 = "feh2s", emi_sectors = "indst"), dim = 3.3, na.rm = TRUE) * tdptwyr2dpgj, "Price|Final Energy|Industry|Hydrogen (US$2017/GJ)"), - setNames( dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1="fehos", emi_sectors = "indst"), dim=3.3, na.rm = T)*tdptwyr2dpgj, + setNames(dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1 = "fehos", emi_sectors = "indst"), dim = 3.3, na.rm = TRUE) * tdptwyr2dpgj, "Price|Final Energy|Industry|Liquids (US$2017/GJ)"), - setNames( dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1="fehes", emi_sectors = "indst"), dim=3.3, na.rm = T)*tdptwyr2dpgj, + setNames(dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1 = "fehes", emi_sectors = "indst"), dim = 3.3, na.rm = TRUE) * tdptwyr2dpgj, "Price|Final Energy|Industry|Heat (US$2017/GJ)"), - setNames( dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1="fesos", emi_sectors = "indst"), dim=3.3, na.rm = T)*tdptwyr2dpgj, + setNames(dimSums(mselect(p_weights_FEprice_mkt * pm_FEPrice, all_enty1 = "fesos", emi_sectors = "indst"), dim = 3.3, na.rm = TRUE) * tdptwyr2dpgj, "Price|Final Energy|Industry|Solids (US$2017/GJ)") - ) + ) ## SE Prices out <- mbind(out, - setNames(mselect(pm_SEPrice, all_enty="seliqfos")*tdptwyr2dpgj, + setNames(mselect(pm_SEPrice, all_enty = "seliqfos") * tdptwyr2dpgj, "Price|Secondary Energy|Liquids|Fossil (US$2017/GJ)"), - setNames(mselect(pm_SEPrice, all_enty="seliqbio")*tdptwyr2dpgj, + setNames(mselect(pm_SEPrice, all_enty = "seliqbio") * tdptwyr2dpgj, "Price|Secondary Energy|Liquids|Biomass (US$2017/GJ)"), - setNames(mselect(pm_SEPrice, all_enty="seliqsyn")*tdptwyr2dpgj, + setNames(mselect(pm_SEPrice, all_enty = "seliqsyn") * tdptwyr2dpgj, "Price|Secondary Energy|Liquids|Hydrogen (US$2017/GJ)"), - setNames(mselect(pm_SEPrice, all_enty="sesobio")*tdptwyr2dpgj, + setNames(mselect(pm_SEPrice, all_enty = "sesobio") * tdptwyr2dpgj, "Price|Secondary Energy|Solids|Biomass (US$2017/GJ)"), - setNames(mselect(pm_SEPrice, all_enty="sesofos")*tdptwyr2dpgj, + setNames(mselect(pm_SEPrice, all_enty = "sesofos") * tdptwyr2dpgj, "Price|Secondary Energy|Solids|Fossil (US$2017/GJ)"), - setNames(mselect(pm_SEPrice, all_enty="seel")*tdptwyr2dpgj, + setNames(mselect(pm_SEPrice, all_enty = "seel") * tdptwyr2dpgj, "Price|Secondary Energy|Electricity (US$2017/GJ)"), - setNames(mselect(pm_SEPrice, all_enty="seh2")*tdptwyr2dpgj, + setNames(mselect(pm_SEPrice, all_enty = "seh2") * tdptwyr2dpgj, "Price|Secondary Energy|Hydrogen (US$2017/GJ)"), - setNames(mselect(pm_SEPrice, all_enty="segabio")*tdptwyr2dpgj, + setNames(mselect(pm_SEPrice, all_enty = "segabio") * tdptwyr2dpgj, "Price|Secondary Energy|Gases|Biomass (US$2017/GJ)"), - setNames(mselect(pm_SEPrice, all_enty="segafos")*tdptwyr2dpgj, + setNames(mselect(pm_SEPrice, all_enty = "segafos") * tdptwyr2dpgj, "Price|Secondary Energy|Gases|Fossil (US$2017/GJ)"), - setNames(mselect(pm_SEPrice, all_enty="segasyn")*tdptwyr2dpgj, + setNames(mselect(pm_SEPrice, all_enty = "segasyn") * tdptwyr2dpgj, "Price|Secondary Energy|Gases|Hydrogen (US$2017/GJ)"), - setNames(mselect(pm_SEPrice, all_enty="sehe")*tdptwyr2dpgj, + setNames(mselect(pm_SEPrice, all_enty = "sehe") * tdptwyr2dpgj, "Price|Secondary Energy|Heat (US$2017/GJ)") - ) + ) weightsSe <- NULL if (exists("seAgg")) { for (i in seAgg) { # calculate weights based on SE production shares, for all possible secondary energy aggregations (set seAgg in REMIND code) - weightsSe <- mbind(weightsSe,dimSums(mselect(prodSe, all_enty1=unique(filter(seAgg2se ,.data$all_enty==i)[,"all_enty1"])),dim=c(3.1,3.3),na.rm=T)/ - dimSums(mselect(prodSe, all_enty1=filter(seAgg2se ,.data$all_enty==i)[,"all_enty1"]) ,na.rm=T)) + weightsSe <- mbind(weightsSe, dimSums(mselect(prodSe, all_enty1 = unique(filter(seAgg2se, .data$all_enty == i)[, "all_enty1"])), dim = c(3.1, 3.3), na.rm = TRUE) / + dimSums(mselect(prodSe, all_enty1 = filter(seAgg2se, .data$all_enty == i)[, "all_enty1"]), na.rm = TRUE)) } - weightsSe <- weightsSe[,intersect(getYears(weightsSe),getYears(pm_SEPrice)),intersect(getItems(weightsSe,dim=3),getItems(pm_SEPrice,dim=3))] + weightsSe <- weightsSe[, intersect(getYears(weightsSe), getYears(pm_SEPrice)), intersect(getItems(weightsSe, dim = 3), getItems(pm_SEPrice, dim = 3))] out <- mbind(out, - setNames(dimSums(mselect(weightsSe*pm_SEPrice[,,getItems(weightsSe,dim=3)], all_enty1=unique(filter(seAgg2se ,.data$all_enty=="all_seliq")[,"all_enty1"])))*tdptwyr2dpgj, + setNames(dimSums(mselect(weightsSe * pm_SEPrice[, , getItems(weightsSe, dim = 3)], all_enty1 = unique(filter(seAgg2se, .data$all_enty == "all_seliq")[, "all_enty1"]))) * tdptwyr2dpgj, "Price|Secondary Energy|Liquids (US$2017/GJ)"), - setNames(dimSums(mselect(weightsSe*pm_SEPrice[,,getItems(weightsSe,dim=3)], all_enty1=unique(filter(seAgg2se ,.data$all_enty=="all_seso")[,"all_enty1"])))*tdptwyr2dpgj, + setNames(dimSums(mselect(weightsSe * pm_SEPrice[, , getItems(weightsSe, dim = 3)], all_enty1 = unique(filter(seAgg2se, .data$all_enty == "all_seso")[, "all_enty1"]))) * tdptwyr2dpgj, "Price|Secondary Energy|Solids (US$2017/GJ)"), - setNames(dimSums(mselect(weightsSe*pm_SEPrice[,,getItems(weightsSe,dim=3)], all_enty1=unique(filter(seAgg2se ,.data$all_enty=="all_sega")[,"all_enty1"])))*tdptwyr2dpgj, + setNames(dimSums(mselect(weightsSe * pm_SEPrice[, , getItems(weightsSe, dim = 3)], all_enty1 = unique(filter(seAgg2se, .data$all_enty == "all_sega")[, "all_enty1"]))) * tdptwyr2dpgj, "Price|Secondary Energy|Gases (US$2017/GJ)") ) } @@ -277,52 +262,52 @@ reportPrices <- function(gdx, output=NULL, regionSubsetList=NULL, ## PE Prices out <- mbind(out, - setNames(mselect(pm_PEPrice, all_enty="peoil")*tdptwyr2dpgj, + setNames(mselect(pm_PEPrice, all_enty = "peoil") * tdptwyr2dpgj, "Price|Primary Energy|Oil (US$2017/GJ)"), - setNames(mselect(pm_PEPrice, all_enty="pegas")*tdptwyr2dpgj, + setNames(mselect(pm_PEPrice, all_enty = "pegas") * tdptwyr2dpgj, "Price|Primary Energy|Gas (US$2017/GJ)"), - setNames(mselect(pm_PEPrice, all_enty="pecoal")*tdptwyr2dpgj, + setNames(mselect(pm_PEPrice, all_enty = "pecoal") * tdptwyr2dpgj, "Price|Primary Energy|Coal (US$2017/GJ)"), - setNames(mselect(pm_PEPrice, all_enty="peur")*tdptwyr2dpgj, + setNames(mselect(pm_PEPrice, all_enty = "peur") * tdptwyr2dpgj, "Price|Primary Energy|Nuclear (US$2017/GJ)"), ## only modern (ligno-cellulosic) biomass reported - setNames(mselect(pm_PEPrice, all_enty="pebiolc")*tdptwyr2dpgj, + setNames(mselect(pm_PEPrice, all_enty = "pebiolc") * tdptwyr2dpgj, "Price|Primary Energy|Biomass|Modern (US$2017/GJ)"), - setNames(mselect(pm_PEPrice, all_enty="pebios")*tdptwyr2dpgj, + setNames(mselect(pm_PEPrice, all_enty = "pebios") * tdptwyr2dpgj, "Price|Primary Energy|Biomass|1st Generation|Sugar and Starch (US$2017/GJ)"), - setNames(mselect(pm_PEPrice, all_enty="pebios")*tdptwyr2dpgj, + setNames(mselect(pm_PEPrice, all_enty = "pebios") * tdptwyr2dpgj, "Price|Primary Energy|Biomass|1st Generation|Oil-based (US$2017/GJ)") - ) + ) # FE marginal price ---- ## loading values from the model - pm_FEPrice_by_SE_Sector_EmiMkt <- readGDX(gdx, c("pm_FEPrice_by_SE_Sector_EmiMkt","p_FEPrice_by_SE_Sector_EmiMkt"), format="first_found", react = "silent") - pm_FEPrice_by_Sector_EmiMkt <- readGDX(gdx, c("pm_FEPrice_by_Sector_EmiMkt","p_FEPrice_by_Sector_EmiMkt"), format="first_found", react = "silent") - pm_FEPrice_by_SE_Sector <- readGDX(gdx, c("pm_FEPrice_by_SE_Sector","p_FEPrice_by_SE_Sector"), format="first_found", react = "silent") - pm_FEPrice_by_SE_EmiMkt <- readGDX(gdx, c("pm_FEPrice_by_SE_EmiMkt","p_FEPrice_by_SE_EmiMkt"), format="first_found", react = "silent") - pm_FEPrice_by_SE <- readGDX(gdx, c("pm_FEPrice_by_SE","p_FEPrice_by_SE"), format="first_found", react = "silent") - pm_FEPrice_by_Sector <- readGDX(gdx, c("pm_FEPrice_by_Sector","p_FEPrice_by_Sector"), format="first_found", react = "silent") - pm_FEPrice_by_EmiMkt <- readGDX(gdx, c("pm_FEPrice_by_EmiMkt","p_FEPrice_by_EmiMkt"), format="first_found", react = "silent") - pm_FEPrice_by_FE <- readGDX(gdx, c("pm_FEPrice_by_FE","p_FEPrice_by_FE"), format="first_found", react = "silent") - - if(length(pm_FEPrice_by_FE) > 0) { + pm_FEPrice_by_SE_Sector_EmiMkt <- readGDX(gdx, c("pm_FEPrice_by_SE_Sector_EmiMkt", "p_FEPrice_by_SE_Sector_EmiMkt"), format = "first_found", react = "silent") + pm_FEPrice_by_Sector_EmiMkt <- readGDX(gdx, c("pm_FEPrice_by_Sector_EmiMkt", "p_FEPrice_by_Sector_EmiMkt"), format = "first_found", react = "silent") + pm_FEPrice_by_SE_Sector <- readGDX(gdx, c("pm_FEPrice_by_SE_Sector", "p_FEPrice_by_SE_Sector"), format = "first_found", react = "silent") + pm_FEPrice_by_SE_EmiMkt <- readGDX(gdx, c("pm_FEPrice_by_SE_EmiMkt", "p_FEPrice_by_SE_EmiMkt"), format = "first_found", react = "silent") + pm_FEPrice_by_SE <- readGDX(gdx, c("pm_FEPrice_by_SE", "p_FEPrice_by_SE"), format = "first_found", react = "silent") + pm_FEPrice_by_Sector <- readGDX(gdx, c("pm_FEPrice_by_Sector", "p_FEPrice_by_Sector"), format = "first_found", react = "silent") + pm_FEPrice_by_EmiMkt <- readGDX(gdx, c("pm_FEPrice_by_EmiMkt", "p_FEPrice_by_EmiMkt"), format = "first_found", react = "silent") + pm_FEPrice_by_FE <- readGDX(gdx, c("pm_FEPrice_by_FE", "p_FEPrice_by_FE"), format = "first_found", react = "silent") + + if (length(pm_FEPrice_by_FE) > 0) { ## subsetting to those entries used in the model all_enty <- all_enty1 <- NULL - se.fe.sector.emiMkt <- se2fe[,-3] %>% #remove te dimension + se.fe.sector.emiMkt <- se2fe[, -3] %>% # remove te dimension rename(se = all_enty, fe = all_enty1) %>% # rename dimensions left_join(entyFe2Sector %>% rename(fe = all_enty, sector = emi_sectors), by = "fe", relationship = "many-to-many") %>% # adding sectors column left_join(sector2emiMkt %>% rename(emiMkt = all_emiMkt, sector = emi_sectors), by = "sector", relationship = "many-to-many") # adding emiMkt column - pm_FEPrice_by_SE_Sector_EmiMkt <- pm_FEPrice_by_SE_Sector_EmiMkt[,YearsFrom2005,unique(paste(se.fe.sector.emiMkt$se,se.fe.sector.emiMkt$fe,se.fe.sector.emiMkt$sector,se.fe.sector.emiMkt$emiMkt,sep="."))]*tdptwyr2dpgj - pm_FEPrice_by_Sector_EmiMkt <- pm_FEPrice_by_Sector_EmiMkt[,YearsFrom2005,unique(paste( se.fe.sector.emiMkt$fe,se.fe.sector.emiMkt$sector,se.fe.sector.emiMkt$emiMkt,sep="."))]*tdptwyr2dpgj - pm_FEPrice_by_SE_Sector <- pm_FEPrice_by_SE_Sector[,YearsFrom2005,unique(paste(se.fe.sector.emiMkt$se,se.fe.sector.emiMkt$fe,se.fe.sector.emiMkt$sector ,sep="."))]*tdptwyr2dpgj - pm_FEPrice_by_SE_EmiMkt <- pm_FEPrice_by_SE_EmiMkt[,YearsFrom2005,unique(paste(se.fe.sector.emiMkt$se,se.fe.sector.emiMkt$fe ,se.fe.sector.emiMkt$emiMkt,sep="."))]*tdptwyr2dpgj - pm_FEPrice_by_SE <- pm_FEPrice_by_SE[,YearsFrom2005,unique(paste(se.fe.sector.emiMkt$se,se.fe.sector.emiMkt$fe ,sep="."))]*tdptwyr2dpgj - pm_FEPrice_by_Sector <- pm_FEPrice_by_Sector[,YearsFrom2005,unique(paste( se.fe.sector.emiMkt$fe,se.fe.sector.emiMkt$sector ,sep="."))]*tdptwyr2dpgj - pm_FEPrice_by_EmiMkt <- pm_FEPrice_by_EmiMkt[,YearsFrom2005,unique(paste( se.fe.sector.emiMkt$fe ,se.fe.sector.emiMkt$emiMkt,sep="."))]*tdptwyr2dpgj - pm_FEPrice_by_FE <- pm_FEPrice_by_FE[,YearsFrom2005,unique(paste( se.fe.sector.emiMkt$fe ,sep="."))]*tdptwyr2dpgj + pm_FEPrice_by_SE_Sector_EmiMkt <- pm_FEPrice_by_SE_Sector_EmiMkt[, YearsFrom2005, unique(paste(se.fe.sector.emiMkt$se, se.fe.sector.emiMkt$fe, se.fe.sector.emiMkt$sector, se.fe.sector.emiMkt$emiMkt, sep = "."))] * tdptwyr2dpgj + pm_FEPrice_by_Sector_EmiMkt <- pm_FEPrice_by_Sector_EmiMkt[, YearsFrom2005, unique(paste(se.fe.sector.emiMkt$fe, se.fe.sector.emiMkt$sector, se.fe.sector.emiMkt$emiMkt, sep = "."))] * tdptwyr2dpgj + pm_FEPrice_by_SE_Sector <- pm_FEPrice_by_SE_Sector[, YearsFrom2005, unique(paste(se.fe.sector.emiMkt$se, se.fe.sector.emiMkt$fe, se.fe.sector.emiMkt$sector, sep = "."))] * tdptwyr2dpgj + pm_FEPrice_by_SE_EmiMkt <- pm_FEPrice_by_SE_EmiMkt[, YearsFrom2005, unique(paste(se.fe.sector.emiMkt$se, se.fe.sector.emiMkt$fe, se.fe.sector.emiMkt$emiMkt, sep = "."))] * tdptwyr2dpgj + pm_FEPrice_by_SE <- pm_FEPrice_by_SE[, YearsFrom2005, unique(paste(se.fe.sector.emiMkt$se, se.fe.sector.emiMkt$fe, sep = "."))] * tdptwyr2dpgj + pm_FEPrice_by_Sector <- pm_FEPrice_by_Sector[, YearsFrom2005, unique(paste(se.fe.sector.emiMkt$fe, se.fe.sector.emiMkt$sector, sep = "."))] * tdptwyr2dpgj + pm_FEPrice_by_EmiMkt <- pm_FEPrice_by_EmiMkt[, YearsFrom2005, unique(paste(se.fe.sector.emiMkt$fe, se.fe.sector.emiMkt$emiMkt, sep = "."))] * tdptwyr2dpgj + pm_FEPrice_by_FE <- pm_FEPrice_by_FE[, YearsFrom2005, unique(paste(se.fe.sector.emiMkt$fe, sep = "."))] * tdptwyr2dpgj ## reporting naming convention varName <- list( @@ -379,23 +364,23 @@ reportPrices <- function(gdx, output=NULL, regionSubsetList=NULL, ) ## add rawdata price variables, calculated from marginals, to the reporting - addVar <- function(input,var,namevector,fe,se,sector,emiMkt) { # function to add only variables if they were not saved already + addVar <- function(input, var, namevector, fe, se, sector, emiMkt) { # function to add only variables if they were not saved already name <- paste0("Price|Final Energy|", paste(namevector, collapse = "|"), " (US$2017/GJ)") name <- gsub("||", "|", name, fixed = TRUE) name <- gsub("| (", " (", name, fixed = TRUE) if (any(is.na(namevector))) warning("In reportPrices, addVar called with a NA value: ", name) - if (name %in% getItems(input, 3)){ + if (name %in% getItems(input, 3)) { return(NULL) } else { - return(setNames(var[, , paste(c(se,fe,sector,emiMkt),collapse = ".")] , name)) + return(setNames(var[, , paste(c(se, fe, sector, emiMkt), collapse = ".")], name)) } } for (i in 1:nrow(se.fe.sector.emiMkt)) { - curr_fe = se.fe.sector.emiMkt[i,]$fe - curr_se = se.fe.sector.emiMkt[i,]$se - curr_sector = se.fe.sector.emiMkt[i,]$sector - curr_emiMKt = se.fe.sector.emiMkt[i,]$emiMkt + curr_fe <- se.fe.sector.emiMkt[i, ]$fe + curr_se <- se.fe.sector.emiMkt[i, ]$se + curr_sector <- se.fe.sector.emiMkt[i, ]$sector + curr_emiMKt <- se.fe.sector.emiMkt[i, ]$emiMkt out <- mbind(out, addVar(input = out, @@ -412,10 +397,10 @@ reportPrices <- function(gdx, output=NULL, regionSubsetList=NULL, ) ) out <- mbind(out, - addVar(input = out, - var = pm_FEPrice_by_SE_Sector, - namevector = c(varName$sector[curr_sector], varName$fe[curr_fe], varName$se[curr_se]), - fe = curr_fe, se = curr_se, sector = curr_sector, emiMkt = NULL + addVar(input = out, + var = pm_FEPrice_by_SE_Sector, + namevector = c(varName$sector[curr_sector], varName$fe[curr_fe], varName$se[curr_se]), + fe = curr_fe, se = curr_se, sector = curr_sector, emiMkt = NULL ) ) @@ -466,10 +451,10 @@ reportPrices <- function(gdx, output=NULL, regionSubsetList=NULL, "tdh2t", "tdelt", "tdels", "tdhes", "tdbiosos", "tdbiogas", "tdh2s", "tdfoshos", "tdsyngat" ) - vm_costTeCapital <- readGDX(gdx, "vm_costTeCapital", field = "l", restore_zeros = F)[, YearsFrom2005, tech] # [tr USD2017/TWh] - p_teAnnuity <- readGDX(gdx, c("p_teAnnuity","pm_teAnnuity"), restore_zeros = F)[, , tech] - vm_capFac <- readGDX(gdx, "vm_capFac", field = "l", restore_zeros = F)[, YearsFrom2005, tech] * 8760 - pm_data_omf <- readGDX(gdx, "pm_data", restore_zeros = F)[, , "omf"][, , tech] + vm_costTeCapital <- readGDX(gdx, "vm_costTeCapital", field = "l", restore_zeros = FALSE)[, YearsFrom2005, tech] # [tr USD2017/TWh] + p_teAnnuity <- readGDX(gdx, c("p_teAnnuity", "pm_teAnnuity"), restore_zeros = FALSE)[, , tech] + vm_capFac <- readGDX(gdx, "vm_capFac", field = "l", restore_zeros = FALSE)[, YearsFrom2005, tech] * 8760 + pm_data_omf <- readGDX(gdx, "pm_data", restore_zeros = FALSE)[, , "omf"][, , tech] price.investment <- vm_costTeCapital * p_teAnnuity / vm_capFac price.omf <- pm_data_omf * vm_costTeCapital / vm_capFac @@ -519,11 +504,11 @@ reportPrices <- function(gdx, output=NULL, regionSubsetList=NULL, ### Carbon Price Component ---- tech.fossil <- c("tdfospet", "tdfosdie", "tdfosgas", "tdfoshos", "tdfossos") - pm_emifac <- readGDX(gdx, "pm_emifac", field = "l", restore_zeros = F)[, YearsFrom2005, "co2"][, , tech.fossil] # [GtC CO2/TWa] + pm_emifac <- readGDX(gdx, "pm_emifac", field = "l", restore_zeros = FALSE)[, YearsFrom2005, "co2"][, , tech.fossil] # [GtC CO2/TWa] pm_emifac <- pm_emifac * 1e9 / s_twa2mwh / 3.6 # [GtC CO2/TWa] -> [tC CO2/GJ] - p_priceCO2 <- readGDX(gdx,name=c("p_priceCO2","pm_priceCO2"),format="first_found", restore_zeros = F) # [USD2017/tC CO2] - if(length(p_priceCO2) > 0) { + p_priceCO2 <- readGDX(gdx, name = c("p_priceCO2", "pm_priceCO2"), format = "first_found", restore_zeros = FALSE) # [USD2017/tC CO2] + if (length(p_priceCO2) > 0) { p_priceCO2 <- add_columns(p_priceCO2, addnm = setdiff(YearsFrom2005, getYears(p_priceCO2)), dim = 2, fill = NA) p_priceCO2 <- add_columns(p_priceCO2, addnm = setdiff(getRegions(pm_emifac), getRegions(p_priceCO2)), dim = 1, fill = NA) @@ -551,8 +536,8 @@ reportPrices <- function(gdx, output=NULL, regionSubsetList=NULL, "indst.fehos", "indst.fesos", "indst.feels", "indst.feh2s", "indst.fegas", "build.fepet", "indst.fepet" ) - pm_tau_fe_tax <- readGDX(gdx, c("p21_tau_fe_tax","pm_tau_fe_tax"), format="first_found")[, YearsFrom2005, entyFe2Sector] # [tr USD2017/TWa] - pm_tau_fe_sub <- readGDX(gdx, c("p21_tau_fe_sub","pm_tau_fe_sub"), format="first_found")[, YearsFrom2005, entyFe2Sector] # [tr USD2017/TWa] + pm_tau_fe_tax <- readGDX(gdx, c("p21_tau_fe_tax", "pm_tau_fe_tax"), format = "first_found")[, YearsFrom2005, entyFe2Sector] # [tr USD2017/TWa] + pm_tau_fe_sub <- readGDX(gdx, c("p21_tau_fe_sub", "pm_tau_fe_sub"), format = "first_found")[, YearsFrom2005, entyFe2Sector] # [tr USD2017/TWa] price.tax <- (pm_tau_fe_tax + pm_tau_fe_sub) / s_twa2mwh / 3.6 * 1e12 # [USD2017/GJ] out <- mbind( @@ -596,49 +581,49 @@ reportPrices <- function(gdx, output=NULL, regionSubsetList=NULL, se <- c("seliqfos", "segafos", "seliqbio", "segabio", "seliqsyn", "segasyn", "seh2", "seel", "sesofos", "sehe", "sesobio") pm_SEPrice <- readGDX(gdx, "pm_SEPrice")[, YearsFrom2005, se] / s_twa2mwh / 3.6 * 1e12 # [tr USD2017/TWa] -> [USD2017/GJ] - pm_eta_conv <- readGDX(gdx, "pm_eta_conv", restore_zeros = F)[, YearsFrom2005, tech] + pm_eta_conv <- readGDX(gdx, "pm_eta_conv", restore_zeros = FALSE)[, YearsFrom2005, tech] price.fuel <- pm_SEPrice / pm_eta_conv out <- mbind( out, - setNames(price.fuel[,,"seliqfos.tdfospet"], "Price|Final Energy|Transport|Liquids|Petroleum|Fuel Cost (US$2017/GJ)"), - setNames(price.fuel[,,"seliqfos.tdfosdie"], "Price|Final Energy|Transport|Liquids|Diesel|Fuel Cost (US$2017/GJ)"), - setNames(price.fuel[,,"segafos.tdfosgat"], "Price|Final Energy|Transport|Gases|Natural Gas|Fuel Cost (US$2017/GJ)"), - setNames(price.fuel[,,"seliqbio.tdbiopet"], "Price|Final Energy|Transport|Liquids|Biomass|Fuel Cost (US$2017/GJ)"), - setNames(price.fuel[,,"seliqsyn.tdsynpet"], "Price|Final Energy|Transport|Liquids|Efuel|Fuel Cost (US$2017/GJ)"), - setNames(price.fuel[,,"segasyn.tdsyngat"], "Price|Final Energy|Transport|Gases|Efuel|Fuel Cost (US$2017/GJ)"), - setNames(price.fuel[,,"seh2.tdh2t"], "Price|Final Energy|Transport|Hydrogen|Fuel Cost (US$2017/GJ)"), - setNames(price.fuel[,,"seel.tdelt"], "Price|Final Energy|Transport|Electricity|Fuel Cost (US$2017/GJ)"), - - setNames(price.fuel[,,"seliqfos.tdfoshos"], "Price|Final Energy|Buildings|Liquids|Oil|Fuel Cost (US$2017/GJ)"), - setNames(price.fuel[,,"segafos.tdfosgas"], "Price|Final Energy|Buildings|Gases|Natural Gas|Fuel Cost (US$2017/GJ)"), - setNames(price.fuel[,,"sehe.tdhes"], "Price|Final Energy|Buildings|Heat|Fuel Cost (US$2017/GJ)"), - setNames(price.fuel[,,"seliqbio.tdbiopet"], "Price|Final Energy|Buildings|Liquids|Biomass|Fuel Cost (US$2017/GJ)"), - setNames(price.fuel[,,"sesobio.tdbiosos"], "Price|Final Energy|Buildings|Solids|Biomass|Fuel Cost (US$2017/GJ)"), - setNames(price.fuel[,,"segabio.tdbiogas"], "Price|Final Energy|Buildings|Gases|Biomass|Fuel Cost (US$2017/GJ)"), - setNames(price.fuel[,,"seliqsyn.tdsynpet"], "Price|Final Energy|Buildings|Liquids|Efuel|Fuel Cost (US$2017/GJ)"), - setNames(price.fuel[,,"segasyn.tdsyngas"], "Price|Final Energy|Buildings|Gases|Efuel|Fuel Cost (US$2017/GJ)"), - setNames(price.fuel[,,"seh2.tdh2s"], "Price|Final Energy|Buildings|Hydrogen|Fuel Cost (US$2017/GJ)"), - setNames(price.fuel[,,"seel.tdels"], "Price|Final Energy|Buildings|Electricity|Fuel Cost (US$2017/GJ)"), - - setNames(price.fuel[,,"seliqfos.tdfoshos"], "Price|Final Energy|Industry|Liquids|Oil|Fuel Cost (US$2017/GJ)"), - setNames(price.fuel[,,"sesofos.tdfossos"], "Price|Final Energy|Industry|Solids|Coal|Fuel Cost (US$2017/GJ)"), - setNames(price.fuel[,,"segafos.tdfosgas"], "Price|Final Energy|Industry|Gases|Natural Gas|Fuel Cost (US$2017/GJ)"), - setNames(price.fuel[,,"seliqbio.tdbiopet"], "Price|Final Energy|Industry|Liquids|Biomass|Fuel Cost (US$2017/GJ)"), - setNames(price.fuel[,,"sesobio.tdbiosos"], "Price|Final Energy|Industry|Solids|Biomass|Fuel Cost (US$2017/GJ)"), - setNames(price.fuel[,,"segabio.tdbiogas"], "Price|Final Energy|Industry|Gases|Biomass|Fuel Cost (US$2017/GJ)"), - setNames(price.fuel[,,"seliqsyn.tdsynpet"], "Price|Final Energy|Industry|Liquids|Efuel|Fuel Cost (US$2017/GJ)"), - setNames(price.fuel[,,"segasyn.tdsyngas"], "Price|Final Energy|Industry|Gases|Efuel|Fuel Cost (US$2017/GJ)"), - - setNames(price.fuel[,,"seel.tdels"], "Price|Final Energy|Industry|Electricity|Fuel Cost (US$2017/GJ)"), - setNames(price.fuel[,,"seh2.tdh2s"], "Price|Final Energy|Industry|Hydrogen|Fuel Cost (US$2017/GJ)"), - setNames(price.fuel[,,"sehe.tdhes"], "Price|Final Energy|Industry|Heat|Fuel Cost (US$2017/GJ)") + setNames(price.fuel[, , "seliqfos.tdfospet"], "Price|Final Energy|Transport|Liquids|Petroleum|Fuel Cost (US$2017/GJ)"), + setNames(price.fuel[, , "seliqfos.tdfosdie"], "Price|Final Energy|Transport|Liquids|Diesel|Fuel Cost (US$2017/GJ)"), + setNames(price.fuel[, , "segafos.tdfosgat"], "Price|Final Energy|Transport|Gases|Natural Gas|Fuel Cost (US$2017/GJ)"), + setNames(price.fuel[, , "seliqbio.tdbiopet"], "Price|Final Energy|Transport|Liquids|Biomass|Fuel Cost (US$2017/GJ)"), + setNames(price.fuel[, , "seliqsyn.tdsynpet"], "Price|Final Energy|Transport|Liquids|Efuel|Fuel Cost (US$2017/GJ)"), + setNames(price.fuel[, , "segasyn.tdsyngat"], "Price|Final Energy|Transport|Gases|Efuel|Fuel Cost (US$2017/GJ)"), + setNames(price.fuel[, , "seh2.tdh2t"], "Price|Final Energy|Transport|Hydrogen|Fuel Cost (US$2017/GJ)"), + setNames(price.fuel[, , "seel.tdelt"], "Price|Final Energy|Transport|Electricity|Fuel Cost (US$2017/GJ)"), + + setNames(price.fuel[, , "seliqfos.tdfoshos"], "Price|Final Energy|Buildings|Liquids|Oil|Fuel Cost (US$2017/GJ)"), + setNames(price.fuel[, , "segafos.tdfosgas"], "Price|Final Energy|Buildings|Gases|Natural Gas|Fuel Cost (US$2017/GJ)"), + setNames(price.fuel[, , "sehe.tdhes"], "Price|Final Energy|Buildings|Heat|Fuel Cost (US$2017/GJ)"), + setNames(price.fuel[, , "seliqbio.tdbiopet"], "Price|Final Energy|Buildings|Liquids|Biomass|Fuel Cost (US$2017/GJ)"), + setNames(price.fuel[, , "sesobio.tdbiosos"], "Price|Final Energy|Buildings|Solids|Biomass|Fuel Cost (US$2017/GJ)"), + setNames(price.fuel[, , "segabio.tdbiogas"], "Price|Final Energy|Buildings|Gases|Biomass|Fuel Cost (US$2017/GJ)"), + setNames(price.fuel[, , "seliqsyn.tdsynpet"], "Price|Final Energy|Buildings|Liquids|Efuel|Fuel Cost (US$2017/GJ)"), + setNames(price.fuel[, , "segasyn.tdsyngas"], "Price|Final Energy|Buildings|Gases|Efuel|Fuel Cost (US$2017/GJ)"), + setNames(price.fuel[, , "seh2.tdh2s"], "Price|Final Energy|Buildings|Hydrogen|Fuel Cost (US$2017/GJ)"), + setNames(price.fuel[, , "seel.tdels"], "Price|Final Energy|Buildings|Electricity|Fuel Cost (US$2017/GJ)"), + + setNames(price.fuel[, , "seliqfos.tdfoshos"], "Price|Final Energy|Industry|Liquids|Oil|Fuel Cost (US$2017/GJ)"), + setNames(price.fuel[, , "sesofos.tdfossos"], "Price|Final Energy|Industry|Solids|Coal|Fuel Cost (US$2017/GJ)"), + setNames(price.fuel[, , "segafos.tdfosgas"], "Price|Final Energy|Industry|Gases|Natural Gas|Fuel Cost (US$2017/GJ)"), + setNames(price.fuel[, , "seliqbio.tdbiopet"], "Price|Final Energy|Industry|Liquids|Biomass|Fuel Cost (US$2017/GJ)"), + setNames(price.fuel[, , "sesobio.tdbiosos"], "Price|Final Energy|Industry|Solids|Biomass|Fuel Cost (US$2017/GJ)"), + setNames(price.fuel[, , "segabio.tdbiogas"], "Price|Final Energy|Industry|Gases|Biomass|Fuel Cost (US$2017/GJ)"), + setNames(price.fuel[, , "seliqsyn.tdsynpet"], "Price|Final Energy|Industry|Liquids|Efuel|Fuel Cost (US$2017/GJ)"), + setNames(price.fuel[, , "segasyn.tdsyngas"], "Price|Final Energy|Industry|Gases|Efuel|Fuel Cost (US$2017/GJ)"), + + setNames(price.fuel[, , "seel.tdels"], "Price|Final Energy|Industry|Electricity|Fuel Cost (US$2017/GJ)"), + setNames(price.fuel[, , "seh2.tdh2s"], "Price|Final Energy|Industry|Hydrogen|Fuel Cost (US$2017/GJ)"), + setNames(price.fuel[, , "sehe.tdhes"], "Price|Final Energy|Industry|Heat|Fuel Cost (US$2017/GJ)") ) ### Total LCOE ---- .calcLCOE <- function(out, var) { - return(setNames(dimSums(out[, , paste0(var,"|"), pmatch = T], dim = 3, na.rm = T), paste0(var, "|Total LCOE (US$2017/GJ)"))) + return(setNames(dimSums(out[, , paste0(var, "|"), pmatch = TRUE], dim = 3, na.rm = TRUE), paste0(var, "|Total LCOE (US$2017/GJ)"))) } out <- mbind( @@ -686,22 +671,25 @@ reportPrices <- function(gdx, output=NULL, regionSubsetList=NULL, # for cm_startyear and non-SSP2, replace price by average of period before and after # this is a workaround to avoid spikes caused by https://github.com/remindmodel/remind/issues/1068 - if (! grepl("gdp_SSP2", readGDX(gdx, "cm_GDPscen", format = "simplest")) - && cm_startyear > min(getYears(out, as.integer = TRUE))) { + if (!grepl("gdp_SSP2", readGDX(gdx, "cm_GDPscen", format = "simplest")) && + cm_startyear > min(getYears(out, as.integer = TRUE))) { out.reporting[, cm_startyear, ] <- 0.5 * (out[, cm_startyear - 5, ] + out[, cm_startyear + 5, ]) } out.reporting <- lowpass(out.reporting) + # reset values for years smaller than cm_startyear to avoid inconsistencies in cm_startyear - 5 - if (! is.null(gdx_ref)) { + fixedyears <- getYears(out)[getYears(out, as.integer = TRUE) < cm_startyear] + if (!is.null(gdx_ref) && length(fixedyears) > 0) { message("reportPrices loads price for < cm_startyear from gdx_ref.") priceRef <- try(reportPrices(gdx_ref, output = NULL, regionSubsetList = regionSubsetList, t = t)) - fixedyears <- getYears(out)[getYears(out, as.integer = TRUE) < cm_startyear] - if (! inherits(priceRef, "try-error") && length(fixedyears) > 0) { + if (!inherits(priceRef, "try-error")) { joinedNamesRep <- intersect(getNames(out), getNames(priceRef)) - joinedRegions <- intersect(getRegions(priceRef), getRegions(out)) + joinedRegions <- intersect(getItems(priceRef, dim = 1), getItems(out, dim = 1)) out.reporting[joinedRegions, fixedyears, joinedNamesRep] <- priceRef[joinedRegions, fixedyears, joinedNamesRep] joinedNamesRaw <- intersect(getNames(out.rawdata), getNames(priceRef)) out.rawdata[joinedRegions, fixedyears, joinedNamesRaw] <- priceRef[joinedRegions, fixedyears, joinedNamesRaw] + } else { + message("failed to run reportPrices on gdx_ref") } } @@ -715,7 +703,7 @@ reportPrices <- function(gdx, output=NULL, regionSubsetList=NULL, ## add years before cm_startyear (temporary, can be adapted once prices only calculated after cm_startyear in REMIND code) out2 <- new.magpie(getRegions(out), getYears(vm_demFeSector), getNames(out), fill = NA) - out2[,getYears(out), ] <- out + out2[, getYears(out), ] <- out out <- out2 @@ -743,14 +731,14 @@ reportPrices <- function(gdx, output=NULL, regionSubsetList=NULL, if (!is.null(esm2macro.m)) { out <- mbind( out, - setNames(abs(esm2macro.m[,,name_trsp[2]]/(budget.m+1e-10)) * tdptwyr2dpgj , "Price|Energy Service|Transport nonLDV (US$2017/GJ)"), - setNames(abs(esm2macro.m[,,name_trsp[1]]/(budget.m+1e-10)) * tdptwyr2dpgj , "Price|Energy Service|Transport LDV (US$2017/GJ)")) + setNames(abs(esm2macro.m[, , name_trsp[2]] / (budget.m + 1e-10)) * tdptwyr2dpgj, "Price|Energy Service|Transport nonLDV (US$2017/GJ)"), + setNames(abs(esm2macro.m[, , name_trsp[1]] / (budget.m + 1e-10)) * tdptwyr2dpgj, "Price|Energy Service|Transport LDV (US$2017/GJ)")) } # report GHG taxes, differentiated by sector # WARNING: the below sector markup code calculation does not consider regipol sector and emission market specific CO2eq prices. If you use both markups and the model 47 formulations, you will have wrongly reported CO2 sectoral and regional prices. if (all(p21_CO2TaxSectorMarkup == 0)) { # then all are identical - out <- mbind(out, setNames(pm_taxCO2eqSum * 1000 * 12/44, "Price|Carbon (US$2017/t CO2)")) + out <- mbind(out, setNames(pm_taxCO2eqSum * 1000 * 12 / 44, "Price|Carbon (US$2017/t CO2)")) for (pcname in c("Price|Carbon|Demand|Buildings (US$2017/t CO2)", "Price|Carbon|Demand|Transport (US$2017/t CO2)", "Price|Carbon|Demand|Industry (US$2017/t CO2)", "Price|Carbon|Supply (US$2017/t CO2)", "Price|Carbon|AggregatedByGrossCO2 (US$2017/t CO2)")) { @@ -758,103 +746,103 @@ reportPrices <- function(gdx, output=NULL, regionSubsetList=NULL, } } else { # currently p21_CO2TaxSectorMarkup is only implemented for build and trans in REMIND out <- mbind(out, - setNames((1 + p21_CO2TaxSectorMarkup[, , "build"]) * pm_taxCO2eqSum * 1000 * 12/44, - "Price|Carbon|Demand|Buildings (US$2017/t CO2)"), - setNames((1 + p21_CO2TaxSectorMarkup[, , "trans"]) * pm_taxCO2eqSum * 1000 * 12/44, - "Price|Carbon|Demand|Transport (US$2017/t CO2)"), - setNames(pm_taxCO2eqSum * 1000 * 12/44, - "Price|Carbon|Demand|Industry (US$2017/t CO2)"), - setNames(pm_taxCO2eqSum * 1000 * 12/44, - "Price|Carbon|Supply (US$2017/t CO2)") - ) - pm_taxCO2eq_FE <- collapseNames( pm_taxCO2eqSum * (1 + - ( - p21_CO2TaxSectorMarkup[, , "build"] * output_wo_GLO[getRegions(p21_CO2TaxSectorMarkup), , "FE|Buildings (EJ/yr)"] - + p21_CO2TaxSectorMarkup[, , "trans"] * output_wo_GLO[getRegions(p21_CO2TaxSectorMarkup), , "FE|Transport (EJ/yr)"] - ) / output_wo_GLO[getRegions(p21_CO2TaxSectorMarkup), , "FE (EJ/yr)"] - ) ) - pm_taxCO2eq_Emi <- collapseNames( pm_taxCO2eqSum * (1 + - ( - p21_CO2TaxSectorMarkup[, , "build"] * output_wo_GLO[getRegions(p21_CO2TaxSectorMarkup), , "Emi|GHG|Gross|Energy|Demand|Buildings (Mt CO2eq/yr)"] - + p21_CO2TaxSectorMarkup[, , "trans"] * output_wo_GLO[getRegions(p21_CO2TaxSectorMarkup), , "Emi|GHG|Gross|Energy|Demand|Transport (Mt CO2eq/yr)"] - ) / output_wo_GLO[getRegions(p21_CO2TaxSectorMarkup), , "Emi|GHG|Gross|Energy (Mt CO2eq/yr)"] - ) ) - - out <- mbind(out, setNames(pm_taxCO2eq_FE * 1000 * 12/44, - "Price|Carbon (US$2017/t CO2)")) # AggregatedbyFE - out <- mbind(out, setNames(pm_taxCO2eq_Emi * 1000 * 12/44, - "Price|Carbon|AggregatedByGrossCO2 (US$2017/t CO2)")) # AggregatedByEmiGHGGross + setNames((1 + p21_CO2TaxSectorMarkup[, , "build"]) * pm_taxCO2eqSum * 1000 * 12 / 44, + "Price|Carbon|Demand|Buildings (US$2017/t CO2)"), + setNames((1 + p21_CO2TaxSectorMarkup[, , "trans"]) * pm_taxCO2eqSum * 1000 * 12 / 44, + "Price|Carbon|Demand|Transport (US$2017/t CO2)"), + setNames(pm_taxCO2eqSum * 1000 * 12 / 44, + "Price|Carbon|Demand|Industry (US$2017/t CO2)"), + setNames(pm_taxCO2eqSum * 1000 * 12 / 44, + "Price|Carbon|Supply (US$2017/t CO2)") + ) + pm_taxCO2eq_FE <- collapseNames(pm_taxCO2eqSum * (1 + + ( + p21_CO2TaxSectorMarkup[, , "build"] * output_wo_GLO[getRegions(p21_CO2TaxSectorMarkup), , "FE|Buildings (EJ/yr)"] + + p21_CO2TaxSectorMarkup[, , "trans"] * output_wo_GLO[getRegions(p21_CO2TaxSectorMarkup), , "FE|Transport (EJ/yr)"] + ) / output_wo_GLO[getRegions(p21_CO2TaxSectorMarkup), , "FE (EJ/yr)"] + )) + pm_taxCO2eq_Emi <- collapseNames(pm_taxCO2eqSum * (1 + + ( + p21_CO2TaxSectorMarkup[, , "build"] * output_wo_GLO[getRegions(p21_CO2TaxSectorMarkup), , "Emi|GHG|Gross|Energy|Demand|Buildings (Mt CO2eq/yr)"] + + p21_CO2TaxSectorMarkup[, , "trans"] * output_wo_GLO[getRegions(p21_CO2TaxSectorMarkup), , "Emi|GHG|Gross|Energy|Demand|Transport (Mt CO2eq/yr)"] + ) / output_wo_GLO[getRegions(p21_CO2TaxSectorMarkup), , "Emi|GHG|Gross|Energy (Mt CO2eq/yr)"] + )) + + out <- mbind(out, setNames(pm_taxCO2eq_FE * 1000 * 12 / 44, + "Price|Carbon (US$2017/t CO2)")) # AggregatedbyFE + out <- mbind(out, setNames(pm_taxCO2eq_Emi * 1000 * 12 / 44, + "Price|Carbon|AggregatedByGrossCO2 (US$2017/t CO2)")) # AggregatedByEmiGHGGross } peFos <- readGDX(gdx, "peFos") # fossil PE carriers p21_tau_Import <- readGDX(gdx, name = "p21_tau_Import", react = "silent")[, t, peFos] tax_import_type_21 <- readGDX(gdx, name = "tax_import_type_21", react = "silent") - if (! is.null(p21_tau_Import) && ! is.null(tax_import_type_21)) { + if (!is.null(p21_tau_Import) && !is.null(tax_import_type_21)) { pm_taxCO2eqMport <- 0 * pm_taxCO2eqSum if ("CO2taxmarkup" %in% tax_import_type_21) { - pm_taxCO2eqMport <- pm_taxCO2eqMport + dimSums(p21_tau_Import[,,"CO2taxmarkup"], dim = 3.2) * pm_taxCO2eqSum + pm_taxCO2eqMport <- pm_taxCO2eqMport + dimSums(p21_tau_Import[, , "CO2taxmarkup"], dim = 3.2) * pm_taxCO2eqSum } if ("avCO2taxmarkup" %in% tax_import_type_21) { - pm_taxCO2eqMport <- pm_taxCO2eqMport + dimSums(p21_tau_Import[,, "avCO2taxmarkup"], dim = 3.2) * pmax(pm_taxCO2eqSum, magpie_expand(colMeans(pm_taxCO2eqSum), pm_taxCO2eqSum)) + pm_taxCO2eqMport <- pm_taxCO2eqMport + dimSums(p21_tau_Import[, , "avCO2taxmarkup"], dim = 3.2) * pmax(pm_taxCO2eqSum, magpie_expand(colMeans(pm_taxCO2eqSum), pm_taxCO2eqSum)) } - pm_taxCO2eqMport <- pm_taxCO2eqMport * 1000 * 12/44 + pm_taxCO2eqMport <- pm_taxCO2eqMport * 1000 * 12 / 44 # use unweighted average, because weighing according to import volumes might lead to big jumps - out <- mbind(out, setNames(dimSums(pm_taxCO2eqMport, dim = 3.1)/length(peFos), "Price|Carbon|Imported (US$2017/t CO2)")) + out <- mbind(out, setNames(dimSums(pm_taxCO2eqMport, dim = 3.1) / length(peFos), "Price|Carbon|Imported (US$2017/t CO2)")) } # CaptureBal_tmp <- new.magpie(getRegions(out), getYears(out), fill = NA) - CaptureBal_tmp[,getYears(balcapture.m),] <- balcapture.m - out <- mbind(out, setNames(CaptureBal_tmp / (budget.m+1e-10) / 3.66 * 1e3, - "Price|Carbon|Captured (US$2017/t CO2)")) + CaptureBal_tmp[, getYears(balcapture.m), ] <- balcapture.m + out <- mbind(out, setNames(CaptureBal_tmp / (budget.m + 1e-10) / 3.66 * 1e3, + "Price|Carbon|Captured (US$2017/t CO2)")) if (is.null(regionSubsetList$EUR)) { - out <- mbind(out,setNames(pm_taxCO2eq * 1000 * 12/44, "Price|Carbon|EU-wide Regulation For All Sectors (US$2017/t CO2)")) + out <- mbind(out, setNames(pm_taxCO2eq * 1000 * 12 / 44, "Price|Carbon|EU-wide Regulation For All Sectors (US$2017/t CO2)")) } else { - co2EUprice <- pm_taxCO2eq * 1000 * 12/44 - co2EUprice[getRegions(pm_taxCO2eq)[!getRegions(pm_taxCO2eq) %in% regionSubsetList$EUR],,] <- 0 - priceReg <- regionSubsetList$EUR[!regionSubsetList$EUR %in% c("DEU","FRA","EUC","EUW")][1] #select region to define EU price (excluding Germany and France) - for (r in regionSubsetList$EUR){ - co2EUprice[r,,] <- as.vector(co2EUprice[priceReg,,]) + co2EUprice <- pm_taxCO2eq * 1000 * 12 / 44 + co2EUprice[getRegions(pm_taxCO2eq)[!getRegions(pm_taxCO2eq) %in% regionSubsetList$EUR], , ] <- 0 + priceReg <- regionSubsetList$EUR[!regionSubsetList$EUR %in% c("DEU", "FRA", "EUC", "EUW")][1] # select region to define EU price (excluding Germany and France) + for (r in regionSubsetList$EUR) { + co2EUprice[r, , ] <- as.vector(co2EUprice[priceReg, , ]) } - out <- mbind(out,setNames(co2EUprice, "Price|Carbon|EU-wide Regulation For All Sectors (US$2017/t CO2)")) + out <- mbind(out, setNames(co2EUprice, "Price|Carbon|EU-wide Regulation For All Sectors (US$2017/t CO2)")) } # reporting carbon prices considering sectoral and emission market differentiated taxes (it does not consider sectoral CO2 markup formulations) if (!is.null(pm_taxemiMkt)) { - out <- mbind(out,setNames(pm_taxemiMkt[,,"ETS"] * 1000 * 12/44, "Price|Carbon|ETS (US$2017/t CO2)")) - out <- mbind(out,setNames(pm_taxemiMkt[,,"ES"] * 1000 * 12/44, "Price|Carbon|ESR (US$2017/t CO2)")) - if(!is.null(p47_taxCO2eq_AggFE)) { # recalculating carbon prices to take into account emi Mkt taxes - out <- out[,,setdiff(getNames(out), c("Price|Carbon|Demand|Buildings (US$2017/t CO2)","Price|Carbon|Demand|Industry (US$2017/t CO2)", - "Price|Carbon|Demand|Transport (US$2017/t CO2)","Price|Carbon|Supply (US$2017/t CO2)","Price|Carbon (US$2017/t CO2)"))] - out <- mbind(out,setNames(p47_taxCO2eq_SectorAggFE[,,"build"] * 1000 * 12/44, "Price|Carbon|Demand|Buildings (US$2017/t CO2)")) - out <- mbind(out,setNames(p47_taxCO2eq_SectorAggFE[,,"indst"] * 1000 * 12/44, "Price|Carbon|Demand|Industry (US$2017/t CO2)")) - out <- mbind(out,setNames(p47_taxCO2eq_SectorAggFE[,,"trans"] * 1000 * 12/44, "Price|Carbon|Demand|Transport (US$2017/t CO2)")) - out <- mbind(out,setNames(p47_taxCO2eq_AggFE * 1000 * 12/44, "Price|Carbon (US$2017/t CO2)")) - out <- mbind(out,setNames(p47_taxCO2eq_AggFE * 1000 * 12/44, "Price|Carbon|Supply (US$2017/t CO2)")) + out <- mbind(out, setNames(pm_taxemiMkt[, , "ETS"] * 1000 * 12 / 44, "Price|Carbon|ETS (US$2017/t CO2)")) + out <- mbind(out, setNames(pm_taxemiMkt[, , "ES"] * 1000 * 12 / 44, "Price|Carbon|ESR (US$2017/t CO2)")) + if (!is.null(p47_taxCO2eq_AggFE)) { # recalculating carbon prices to take into account emi Mkt taxes + out <- out[, , setdiff(getNames(out), c("Price|Carbon|Demand|Buildings (US$2017/t CO2)", "Price|Carbon|Demand|Industry (US$2017/t CO2)", + "Price|Carbon|Demand|Transport (US$2017/t CO2)", "Price|Carbon|Supply (US$2017/t CO2)", "Price|Carbon (US$2017/t CO2)"))] + out <- mbind(out, setNames(p47_taxCO2eq_SectorAggFE[, , "build"] * 1000 * 12 / 44, "Price|Carbon|Demand|Buildings (US$2017/t CO2)")) + out <- mbind(out, setNames(p47_taxCO2eq_SectorAggFE[, , "indst"] * 1000 * 12 / 44, "Price|Carbon|Demand|Industry (US$2017/t CO2)")) + out <- mbind(out, setNames(p47_taxCO2eq_SectorAggFE[, , "trans"] * 1000 * 12 / 44, "Price|Carbon|Demand|Transport (US$2017/t CO2)")) + out <- mbind(out, setNames(p47_taxCO2eq_AggFE * 1000 * 12 / 44, "Price|Carbon (US$2017/t CO2)")) + out <- mbind(out, setNames(p47_taxCO2eq_AggFE * 1000 * 12 / 44, "Price|Carbon|Supply (US$2017/t CO2)")) } } if (!is.null(pm_taxCO2eqSCC)) { - out <- mbind(out,setNames(abs(pm_taxCO2eqSCC) * 1000 * 12/44, "Price|Carbon|SCC (US$2017/t CO2)")) - out <- mbind(out,setNames(out[, , "Price|Carbon (US$2017/t CO2)"] - out[, , "Price|Carbon|SCC (US$2017/t CO2)"], "Price|Carbon|Guardrail (US$2017/t CO2)")) + out <- mbind(out, setNames(abs(pm_taxCO2eqSCC) * 1000 * 12 / 44, "Price|Carbon|SCC (US$2017/t CO2)")) + out <- mbind(out, setNames(out[, , "Price|Carbon (US$2017/t CO2)"] - out[, , "Price|Carbon|SCC (US$2017/t CO2)"], "Price|Carbon|Guardrail (US$2017/t CO2)")) } else { - out <- mbind(out,setNames(out[, , "Price|Carbon (US$2017/t CO2)"], "Price|Carbon|Guardrail (US$2017/t CO2)")) + out <- mbind(out, setNames(out[, , "Price|Carbon (US$2017/t CO2)"], "Price|Carbon|Guardrail (US$2017/t CO2)")) } - out <- mbind(out,setNames(out[,,"Price|Carbon (US$2017/t CO2)"] * s_GWP_N2O, "Price|N2O (US$2017/t N2O)")) - out <- mbind(out,setNames(out[,,"Price|Carbon (US$2017/t CO2)"] * s_GWP_CH4, "Price|CH4 (US$2017/t CH4)")) + out <- mbind(out, setNames(out[, , "Price|Carbon (US$2017/t CO2)"] * s_GWP_N2O, "Price|N2O (US$2017/t N2O)")) + out <- mbind(out, setNames(out[, , "Price|Carbon (US$2017/t CO2)"] * s_GWP_CH4, "Price|CH4 (US$2017/t CH4)")) # adding extra variables for alternative carbon price aggregation weights - out <- mbind(out,setNames(out[,,"Price|Carbon|Captured (US$2017/t CO2)"], - "Price|Carbon|Captured|AggregatedByGrossCO2 (US$2017/t CO2)")) - out <- mbind(out,setNames(out[,,"Price|Carbon|EU-wide Regulation For All Sectors (US$2017/t CO2)"], - "Price|Carbon|EU-wide Regulation For All Sectors|AggregatedByGrossCO2 (US$2017/t CO2)")) - out <- mbind(out,setNames(out[,,"Price|Carbon|Guardrail (US$2017/t CO2)"], - "Price|Carbon|Guardrail|AggregatedByGrossCO2 (US$2017/t CO2)")) - out <- mbind(out,setNames(out[,,"Price|Carbon|SCC (US$2017/t CO2)"], - "Price|Carbon|SCC|AggregatedByGrossCO2 (US$2017/t CO2)")) + out <- mbind(out, setNames(out[, , "Price|Carbon|Captured (US$2017/t CO2)"], + "Price|Carbon|Captured|AggregatedByGrossCO2 (US$2017/t CO2)")) + out <- mbind(out, setNames(out[, , "Price|Carbon|EU-wide Regulation For All Sectors (US$2017/t CO2)"], + "Price|Carbon|EU-wide Regulation For All Sectors|AggregatedByGrossCO2 (US$2017/t CO2)")) + out <- mbind(out, setNames(out[, , "Price|Carbon|Guardrail (US$2017/t CO2)"], + "Price|Carbon|Guardrail|AggregatedByGrossCO2 (US$2017/t CO2)")) + out <- mbind(out, setNames(out[, , "Price|Carbon|SCC (US$2017/t CO2)"], + "Price|Carbon|SCC|AggregatedByGrossCO2 (US$2017/t CO2)")) # ---- mapping of weights for the variables for global aggregation ---- @@ -909,7 +897,7 @@ reportPrices <- function(gdx, output=NULL, regionSubsetList=NULL, "Price|Carbon|EU-wide Regulation For All Sectors|AggregatedByGrossCO2 (US$2017/t CO2)" = "Emi|GHG|Gross|Energy (Mt CO2eq/yr)", "Price|Carbon|Guardrail|AggregatedByGrossCO2 (US$2017/t CO2)" = "Emi|GHG|Gross|Energy (Mt CO2eq/yr)", "Price|Carbon|SCC|AggregatedByGrossCO2 (US$2017/t CO2)" = "Emi|GHG|Gross|Energy (Mt CO2eq/yr)" - ) + ) if (!is.null(esm2macro.m)) { int2ext <- c(int2ext, @@ -941,14 +929,14 @@ reportPrices <- function(gdx, output=NULL, regionSubsetList=NULL, "Price|Final Energy|Industry|Hydrogen (US$2017/GJ)" = "FE|Industry|Hydrogen (EJ/yr)", "Price|Final Energy|Industry|Heat (US$2017/GJ)" = "FE|Industry|Heat (EJ/yr)", "Price|Final Energy|Industry|Solids (US$2017/GJ)" = "FE|Industry|Solids (EJ/yr)" - ) + ) # transport-specific mappings depending on realization - if (module2realisation["transport",2] == "complex") { + if (module2realisation["transport", 2] == "complex") { int2ext <- c(int2ext, "Price|Final Energy|Transport|Liquids|HDV (US$2017/GJ)" = "FE|Transport|non-LDV|Liquids (EJ/yr)", "Price|Final Energy|Transport|Liquids|LDV (US$2017/GJ)" = "FE|Transport|LDV|Liquids (EJ/yr)") - } else if (module2realisation["transport",2] == "edge_esm") { + } else if (module2realisation["transport", 2] == "edge_esm") { int2ext <- c(int2ext, "Price|Final Energy|Transport|Liquids|HDV (US$2017/GJ)" = "FE|Transport|Diesel Liquids (EJ/yr)", "Price|Final Energy|Transport|Liquids|LDV (US$2017/GJ)" = "FE|Transport|Pass|Liquids (EJ/yr)") @@ -956,16 +944,16 @@ reportPrices <- function(gdx, output=NULL, regionSubsetList=NULL, ## add weights definition for region aggregation for FE prices that were added automatically if (length(pm_FEPrice_by_FE) > 0) { - margPriceVars <- grep("Price|Final Energy|", getItems(out,3), fixed = TRUE, value = TRUE) + margPriceVars <- grep("Price|Final Energy|", getItems(out, 3), fixed = TRUE, value = TRUE) margPriceVars <- setdiff(margPriceVars, names(int2ext)) - vars <- gsub("US\\$2017/GJ", "EJ/yr", gsub("Price\\|Final Energy\\|","FE|",margPriceVars)) + vars <- gsub("US\\$2017/GJ", "EJ/yr", gsub("Price\\|Final Energy\\|", "FE|", margPriceVars)) names(vars) <- margPriceVars - vars <- gsub("Efuel","Hydrogen",vars) ###warning FE variable should be renamed and this line should be removed in the future - vars <- vars[vars %in% getItems(output,3)] + vars <- gsub("Efuel", "Hydrogen", vars) ### warning FE variable should be renamed and this line should be removed in the future + vars <- vars[vars %in% getItems(output, 3)] int2ext <- c(int2ext, vars) } - int2ext <- int2ext[! is.na(int2ext)] + int2ext <- int2ext[!is.na(int2ext)] ## add subvariables, plus moving averages and rawdata. .addSubvariable is defined at the bottom of this file int2ext <- .addSubvariable(int2ext, c("|Fuel Cost", "|Other Taxes", "|Total LCOE", @@ -977,159 +965,120 @@ reportPrices <- function(gdx, output=NULL, regionSubsetList=NULL, # # CES prices and CES markup cost - o01_CESderivatives <- readGDX(gdx, "o01_CESderivatives", restore_zeros = T, react = "silent") + o01_CESderivatives <- readGDX(gdx, "o01_CESderivatives", restore_zeros = TRUE, react = "silent") if (!is.null(o01_CESderivatives)) { - o01_CESderivatives <- o01_CESderivatives[,YearsFrom2005,] + o01_CESderivatives <- o01_CESderivatives[, YearsFrom2005, ] - # CES price is derivative of GDP with respect to production factor (CES node) - # temporay: only report CES prices of primary production factors (bottom-most CES nodes) for buildings and industry for now to not oversize the reporting - # TODO: generate MIF without internal variables as standard such that more internal variables can be added here + # CES price is derivative of GDP with respect to production factor (CES node) + # temporay: only report CES prices of primary production factors (bottom-most CES nodes) for buildings and industry for now to not oversize the reporting + # TODO: generate MIF without internal variables as standard such that more internal variables can be added here - ppfen_industry_dyn37 <- readGDX(gdx, "ppfen_industry_dyn37") - ppfen_buildings_dyn36 <- readGDX(gdx, "ppfen_buildings_dyn36") - ppfen <- c(ppfen_industry_dyn37, ppfen_buildings_dyn36) + ppfen_industry_dyn37 <- readGDX(gdx, "ppfen_industry_dyn37") + ppfen_buildings_dyn36 <- readGDX(gdx, "ppfen_buildings_dyn36") + ppfen <- c(ppfen_industry_dyn37, ppfen_buildings_dyn36) - # CES Prices + # CES Prices - # choose derivative of GDP (inco) with respect to input - ces_price <- collapseDim(mselect(o01_CESderivatives, all_in = "inco", all_in1 = ppfen)) - # variable names - ces_price <- setNames(ces_price, paste0("Internal|Price|CES|",getNames(ces_price)," (tr US$2017/input unit)")) + # choose derivative of GDP (inco) with respect to input + ces_price <- collapseDim(mselect(o01_CESderivatives, all_in = "inco", all_in1 = ppfen)) + # variable names + ces_price <- setNames(ces_price, paste0("Internal|Price|CES|", getNames(ces_price), " (tr US$2017/input unit)")) - # CES Markup Cost - p37_CESMkup <- readGDX(gdx, "p37_CESMkup") # markup in industry - p36_CESMkup <- readGDX(gdx, "p36_CESMkup") # markup in buildings + # CES Markup Cost + p37_CESMkup <- readGDX(gdx, "p37_CESMkup") # markup in industry + p36_CESMkup <- readGDX(gdx, "p36_CESMkup") # markup in buildings - CESMkup <- mbind( - mselect(p36_CESMkup[,YearsFrom2005,], all_in = ppfen_buildings_dyn36), - mselect(p37_CESMkup[,YearsFrom2005,], all_in = ppfen_industry_dyn37)) + CESMkup <- mbind( + mselect(p36_CESMkup[, YearsFrom2005, ], all_in = ppfen_buildings_dyn36), + mselect(p37_CESMkup[, YearsFrom2005, ], all_in = ppfen_industry_dyn37)) - CESMkup <- setNames( CESMkup, paste0("Internal|CES Markup Cost|",getNames(CESMkup)," (tr US$2017/input unit)")) + CESMkup <- setNames(CESMkup, paste0("Internal|CES Markup Cost|", getNames(CESMkup), " (tr US$2017/input unit)")) - out <- mbind(out, ces_price, CESMkup) + out <- mbind(out, ces_price, CESMkup) } # add global prices - map <- data.frame(region=getRegions(out),world="GLO",stringsAsFactors=FALSE) + map <- data.frame(region = getRegions(out), world = "GLO", stringsAsFactors = FALSE) tmp_GLO <- new.magpie("GLO", getYears(out), getNames(out), fill = NA) int2ext <- int2ext[intersect(names(int2ext), getNames(out))] # select only data that exists for (i2e in names(int2ext)) { - tmp_GLO["GLO",,i2e] <- toolAggregate(out[,,i2e], rel = map, weight = output[map$region,,int2ext[i2e]], zeroWeight = "setNA") + tmp_GLO["GLO", , i2e] <- toolAggregate(out[, , i2e], rel = map, weight = output[map$region, , int2ext[i2e]], zeroWeight = "setNA") } out <- mbind(out, tmp_GLO) # add other region aggregations - if (!is.null(regionSubsetList)){ - tmp_RegAgg <- new.magpie(names(regionSubsetList),getYears(out),getNames(out),fill=0) - for(region in names(regionSubsetList)){ - tmp_RegAgg_ie2 <- do.call("mbind",lapply(names(int2ext), function(i2e) { - map <- data.frame(region=regionSubsetList[[region]],parentRegion=region,stringsAsFactors=FALSE) - result <- toolAggregate(out[regionSubsetList[[region]],,i2e], rel = map, weight = output[regionSubsetList[[region]],,int2ext[i2e]], zeroWeight = "allow") + if (!is.null(regionSubsetList)) { + tmp_RegAgg <- new.magpie(names(regionSubsetList), getYears(out), getNames(out), fill = 0) + for (region in names(regionSubsetList)) { + tmp_RegAgg_ie2 <- do.call("mbind", lapply(names(int2ext), function(i2e) { + map <- data.frame(region = regionSubsetList[[region]], parentRegion = region, stringsAsFactors = FALSE) + result <- toolAggregate(out[regionSubsetList[[region]], , i2e], rel = map, weight = output[regionSubsetList[[region]], , int2ext[i2e]], zeroWeight = "allow") getItems(result, dim = 1) <- region - for(t in getYears(out)){ - if(all(output[regionSubsetList[[region]],t,int2ext[i2e]]==0)){ - result[region,t,i2e] <- NA + for (t in getYears(out)) { + if (all(output[regionSubsetList[[region]], t, int2ext[i2e]] == 0)) { + result[region, t, i2e] <- NA } } return(result) })) - tmp_RegAgg[region,,names(int2ext)] <- tmp_RegAgg_ie2[region,,names(int2ext)] + tmp_RegAgg[region, , names(int2ext)] <- tmp_RegAgg_ie2[region, , names(int2ext)] } out <- mbind(out, tmp_RegAgg) } # prices that are same for all regions, including GLO - glob_price <- new.magpie(getRegions(out),getYears(out),fill=0) - for(i in getRegions(out)) glob_price[i,,] <- pm_pvp[,,"peur"]*1000 - out <- mbind(out,setNames(glob_price, "PVP1|Uranium (billionDpktU)")) - for(i in getRegions(out)) glob_price[i,,] <- pm_pvp[,,"peoil"]*1000 - out <- mbind(out,setNames(glob_price, "PVP1|Oil (billionDpTWyr)")) - for(i in getRegions(out)) glob_price[i,,] <- pm_pvp[,,"pegas"]*1000 - out <- mbind(out,setNames(glob_price, "PVP1|Gas (billionDpTWyr)")) - for(i in getRegions(out)) glob_price[i,,] <- pm_pvp[,,"pecoal"]*1000 - out <- mbind(out,setNames(glob_price, "PVP1|Coal (billionDpTWyr)")) - for(i in getRegions(out)) glob_price[i,,] <- pm_pvp[,,"pebiolc"]*1000 - out <- mbind(out,setNames(glob_price, "PVP1|Biomass (billionDpTWyr)")) - for(i in getRegions(out)) glob_price[i,,] <- pm_pvp[,,"good"]*1000 - out <- mbind(out,setNames(glob_price, "PVP2|Good ()")) - for(i in getRegions(out)) glob_price[i,,] <- pm_pvp[,,"perm"] - out <- mbind(out,setNames(glob_price, "PVP3|Permit ()")) - - for(i in getRegions(out)) glob_price[i,,] <- pm_pvp[,,"peur"]/pm_pvp[,,"good"] * tdptwyr2dpgj - out <- mbind(out,setNames(glob_price, "Price|Uranium|World Market (US$2017/GJ)")) - for(i in getRegions(out)) glob_price[i,,] <- pm_pvp[,,"peoil"]/pm_pvp[,,"good"] * tdptwyr2dpgj - out <- mbind(out,setNames(glob_price, "Price|Oil|World Market (US$2017/GJ)")) - for(i in getRegions(out)) glob_price[i,,] <- pm_pvp[,,"pegas"]/pm_pvp[,,"good"] * tdptwyr2dpgj - out <- mbind(out,setNames(glob_price, "Price|Gas|World Market (US$2017/GJ)")) - for(i in getRegions(out)) glob_price[i,,] <- pm_pvp[,,"pecoal"]/pm_pvp[,,"good"] * tdptwyr2dpgj - out <- mbind(out,setNames(glob_price, "Price|Coal|World Market (US$2017/GJ)")) - for(i in getRegions(out)) glob_price[i,,] <- pm_pvp[,,"pebiolc"]/pm_pvp[,,"good"] * tdptwyr2dpgj - out <- mbind(out,setNames(glob_price, "Price|Biomass|World Market (US$2017/GJ)")) + glob_price <- new.magpie(getRegions(out), getYears(out), fill = 0) + for (i in getRegions(out)) glob_price[i, , ] <- pm_pvp[, , "peur"] * 1000 + out <- mbind(out, setNames(glob_price, "PVP1|Uranium (billionDpktU)")) + for (i in getRegions(out)) glob_price[i, , ] <- pm_pvp[, , "peoil"] * 1000 + out <- mbind(out, setNames(glob_price, "PVP1|Oil (billionDpTWyr)")) + for (i in getRegions(out)) glob_price[i, , ] <- pm_pvp[, , "pegas"] * 1000 + out <- mbind(out, setNames(glob_price, "PVP1|Gas (billionDpTWyr)")) + for (i in getRegions(out)) glob_price[i, , ] <- pm_pvp[, , "pecoal"] * 1000 + out <- mbind(out, setNames(glob_price, "PVP1|Coal (billionDpTWyr)")) + for (i in getRegions(out)) glob_price[i, , ] <- pm_pvp[, , "pebiolc"] * 1000 + out <- mbind(out, setNames(glob_price, "PVP1|Biomass (billionDpTWyr)")) + for (i in getRegions(out)) glob_price[i, , ] <- pm_pvp[, , "good"] * 1000 + out <- mbind(out, setNames(glob_price, "PVP2|Good ()")) + for (i in getRegions(out)) glob_price[i, , ] <- pm_pvp[, , "perm"] + out <- mbind(out, setNames(glob_price, "PVP3|Permit ()")) + + for (i in getRegions(out)) glob_price[i, , ] <- pm_pvp[, , "peur"] / pm_pvp[, , "good"] * tdptwyr2dpgj + out <- mbind(out, setNames(glob_price, "Price|Uranium|World Market (US$2017/GJ)")) + for (i in getRegions(out)) glob_price[i, , ] <- pm_pvp[, , "peoil"] / pm_pvp[, , "good"] * tdptwyr2dpgj + out <- mbind(out, setNames(glob_price, "Price|Oil|World Market (US$2017/GJ)")) + for (i in getRegions(out)) glob_price[i, , ] <- pm_pvp[, , "pegas"] / pm_pvp[, , "good"] * tdptwyr2dpgj + out <- mbind(out, setNames(glob_price, "Price|Gas|World Market (US$2017/GJ)")) + for (i in getRegions(out)) glob_price[i, , ] <- pm_pvp[, , "pecoal"] / pm_pvp[, , "good"] * tdptwyr2dpgj + out <- mbind(out, setNames(glob_price, "Price|Coal|World Market (US$2017/GJ)")) + for (i in getRegions(out)) glob_price[i, , ] <- pm_pvp[, , "pebiolc"] / pm_pvp[, , "good"] * tdptwyr2dpgj + out <- mbind(out, setNames(glob_price, "Price|Biomass|World Market (US$2017/GJ)")) ## special global prices ## not meaningful global prices set to NA - out["GLO",,"Internal|Price|Biomass|Shiftfactor ()"] <- NA + out["GLO", , "Internal|Price|Biomass|Shiftfactor ()"] <- NA if (!is.null(regionSubsetList$EUR)) - out["EUR",,"Price|Carbon|EU-wide Regulation For All Sectors (US$2017/t CO2)"] <- as.vector(co2EUprice[priceReg,,]) + out["EUR", , "Price|Carbon|EU-wide Regulation For All Sectors (US$2017/t CO2)"] <- as.vector(co2EUprice[priceReg, , ]) ## not meaningful region aggregation prices set to NA - if (!is.null(regionSubsetList)){ - for(region in names(regionSubsetList)){ - out[region,,"Internal|Price|Biomass|Shiftfactor ()"] <- NA + if (!is.null(regionSubsetList)) { + for (region in names(regionSubsetList)) { + out[region, , "Internal|Price|Biomass|Shiftfactor ()"] <- NA } } - # comment out this section for now as errors, debug this section if needed - # # ---- debug information for industry/subsectors ---- - # if ('subsectors' == indu_mod & !is.null(q37_limit_secondary_steel_share.m)) { - # - # t <- getYears(budget.m) - # - # .x <- q37_limit_secondary_steel_share.m[, t,] / budget.m - # .x <- mbind(.x, calc_regionSubset_sums(.x, regionSubsetList)) - # - # - # tmp2 <- mbind( - # # fake some GLO data - # setNames( - # mbind(.x, dimSums(.x * NA, dim = 1)), - # 'Debug|Industry|Secondary Steel Premium (US$2017)'), - # - # mbind( - # lapply( - # list( - # c('ue_industry', '', 'arbitrary unit', 1), - # c('ue_cement', '|Cement', 't cement', 1e3), - # c('ue_chemicals', '|Chemicals', 'arbitrary unit', 1), - # c('ue_steel', '|Steel', 't Steel', 1e3), - # c('ue_steel_primary', '|Steel|Primary', 't Steel', 1e3), - # c('ue_steel_secondary', '|Steel|Secondary', 't Steel', 1e3), - # c('ue_otherInd', '|other', 'arbitrary unit', 1)), - # function(x) { - # setNames( - # ( out[,,paste0('Price|CES_input|', x[1], ' ('), pmatch = 'left'] - # * as.numeric(x[4]) - # ), - # paste0('Debug|Price|Industry', x[2], ' (US$2017/', x[3], ')')) - # }) - # ) - # ) - # - # out <- mbind(out, tmp2) - # } - getSets(out)[3] <- "variable" return(out) } diff --git a/README.md b/README.md index 624ae878..324dd9bf 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # The REMIND R package (2nd generation) -R package **remind2**, version **1.158.2** +R package **remind2**, version **1.159.0** [![CRAN status](https://www.r-pkg.org/badges/version/remind2)](https://cran.r-project.org/package=remind2) [![R build status](https://github.com/pik-piam/remind2/workflows/check/badge.svg)](https://github.com/pik-piam/remind2/actions) [![codecov](https://codecov.io/gh/pik-piam/remind2/branch/master/graph/badge.svg)](https://app.codecov.io/gh/pik-piam/remind2) [![r-universe](https://pik-piam.r-universe.dev/badges/remind2)](https://pik-piam.r-universe.dev/builds) @@ -49,7 +49,7 @@ In case of questions / problems please contact Renato Rodrigues . +Rodrigues R, Baumstark L, Benke F, Dietrich J, Dirnaichner A, Duerrwaechter J, Führlich P, Giannousakis A, Hasse R, Hilaire J, Klein D, Koch J, Kowalczyk K, Levesque A, Malik A, Merfort A, Merfort L, Morena-Leiva S, Pehl M, Pietzcker R, Rauner S, Richters O, Rottoli M, Schötz C, Schreyer F, Siala K, Sörgel B, Spahr M, Strefler J, Verpoort P, Weigmann P, Rüter T (2024). _remind2: The REMIND R package (2nd generation)_. R package version 1.159.0, . A BibTeX entry for LaTeX users is @@ -58,7 +58,7 @@ A BibTeX entry for LaTeX users is title = {remind2: The REMIND R package (2nd generation)}, author = {Renato Rodrigues and Lavinia Baumstark and Falk Benke and Jan Philipp Dietrich and Alois Dirnaichner and Jakob Duerrwaechter and Pascal Führlich and Anastasis Giannousakis and Robin Hasse and Jérome Hilaire and David Klein and Johannes Koch and Katarzyna Kowalczyk and Antoine Levesque and Aman Malik and Anne Merfort and Leon Merfort and Simón Morena-Leiva and Michaja Pehl and Robert Pietzcker and Sebastian Rauner and Oliver Richters and Marianna Rottoli and Christof Schötz and Felix Schreyer and Kais Siala and Björn Sörgel and Mike Spahr and Jessica Strefler and Philipp Verpoort and Pascal Weigmann and Tonn Rüter}, year = {2024}, - note = {R package version 1.158.2}, + note = {R package version 1.159.0}, url = {https://github.com/pik-piam/remind2}, } ``` diff --git a/man/convGDX2MIF.Rd b/man/convGDX2MIF.Rd index 7a863962..1dc7bbbe 100644 --- a/man/convGDX2MIF.Rd +++ b/man/convGDX2MIF.Rd @@ -37,8 +37,9 @@ Read in all information from GDX file and create the *.mif reporting } \examples{ - -\dontrun{convGDX2MIF(gdx,gdx_refpolicycost,file="REMIND_generic_default.csv",scenario="default")} +\dontrun{ +convGDX2MIF(gdx, gdx_refpolicycost, file = "REMIND_generic_default.csv", scenario = "default") +} } \author{ diff --git a/man/convGDX2MIF_LCOE.Rd b/man/convGDX2MIF_LCOE.Rd index 83c3216c..090db5d7 100644 --- a/man/convGDX2MIF_LCOE.Rd +++ b/man/convGDX2MIF_LCOE.Rd @@ -6,7 +6,6 @@ \usage{ convGDX2MIF_LCOE( gdx, - gdx_ref, file = NULL, scenario = "default", t = c(seq(2005, 2060, 5), seq(2070, 2110, 10), 2130, 2150) @@ -15,8 +14,6 @@ convGDX2MIF_LCOE( \arguments{ \item{gdx}{a GDX as created by readGDX, or the file name of a gdx} -\item{gdx_ref}{a GDX as created by readGDX of the reference run} - \item{file}{name of the mif file which will be written, if no name is provided a magpie object containing all the reporting information is returned} @@ -32,7 +29,7 @@ the LCOE .mif reporting } \examples{ \dontrun{ -convGDX2MIF(gdx, gdx_ref, file = "REMIND_generic_LCOE.csv", scenario = "default") +convGDX2MIF_LCOE(gdx, file = "REMIND_generic_LCOE.csv", scenario = "default") } } diff --git a/man/modifyInvestmentVariables.Rd b/man/modifyInvestmentVariables.Rd new file mode 100644 index 00000000..c72a8a38 --- /dev/null +++ b/man/modifyInvestmentVariables.Rd @@ -0,0 +1,30 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/modifyInvestmentVariables.R +\name{modifyInvestmentVariables} +\alias{modifyInvestmentVariables} +\title{Modify Investment Variables} +\usage{ +modifyInvestmentVariables(x, ref = NULL, startYear = NULL) +} +\arguments{ +\item{x}{a magclass object to be manipulated, must have timesteps in 'ttot'} + +\item{ref}{an optional magclass object to be used for fixing values before 'startYear'} + +\item{startYear}{years before will be overwritten with values from 'ref'} +} +\description{ +This function transforms investment variables to the normal reporting convention. +} +\details{ +Years represented by investment variables in the energy system +('vm_deltaCap', 'vm_costInvTeDir' and 'vm_costInvTeAdj') are different from the normal +reporting convention. In the current REMIND version, vm_deltacap(t) represents +the average of the years t-4..t, while in the general reporting +convention it represents the average of t-2.5..t+2.5 (for 5 year time steps). + +See also: https://github.com/remindmodel/remind/pull/1238. +} +\author{ +Falk Benke +} diff --git a/man/reportCapacity.Rd b/man/reportCapacity.Rd index e9003306..4ff4ca80 100644 --- a/man/reportCapacity.Rd +++ b/man/reportCapacity.Rd @@ -7,7 +7,8 @@ reportCapacity( gdx, regionSubsetList = NULL, - t = c(seq(2005, 2060, 5), seq(2070, 2110, 10), 2130, 2150) + t = c(seq(2005, 2060, 5), seq(2070, 2110, 10), 2130, 2150), + gdx_ref = gdx_ref ) } \arguments{ @@ -19,6 +20,10 @@ be created.} \item{t}{temporal resolution of the reporting, default: t=c(seq(2005,2060,5),seq(2070,2110,10),2130,2150)} + +\item{gdx_ref}{a GDX object as created by readGDX, or the path to a gdx of the reference run. +It is used to guarantee consistency before 'cm_startyear' for capacity variables +using time averaging.} } \value{ MAgPIE object - contains the capacity variables diff --git a/man/reportCapitalStock.Rd b/man/reportCapitalStock.Rd index e74c52a7..dedaa1ce 100644 --- a/man/reportCapitalStock.Rd +++ b/man/reportCapitalStock.Rd @@ -7,7 +7,8 @@ reportCapitalStock( gdx, regionSubsetList = NULL, - t = c(seq(2005, 2060, 5), seq(2070, 2110, 10), 2130, 2150) + t = c(seq(2005, 2060, 5), seq(2070, 2110, 10), 2130, 2150), + gdx_ref = gdx_ref ) } \arguments{ @@ -19,6 +20,10 @@ be created.} \item{t}{temporal resolution of the reporting, default: t=c(seq(2005,2060,5),seq(2070,2110,10),2130,2150)} + +\item{gdx_ref}{a GDX object as created by readGDX, or the path to a gdx of the reference run. +It is used to guarantee consistency before 'cm_startyear' for capacity variables +using time averaging.} } \value{ MAgPIE object - contains the capital stock variables @@ -28,8 +33,9 @@ Read in capital stock information from GDX file, information used in convGDX2MIF for the reporting } \examples{ - -\dontrun{reportCapitalStock(gdx)} +\dontrun{ +reportCapitalStock(gdx) +} } \seealso{ diff --git a/man/reportEnergyInvestment.Rd b/man/reportEnergyInvestment.Rd index dea8d654..638966d9 100644 --- a/man/reportEnergyInvestment.Rd +++ b/man/reportEnergyInvestment.Rd @@ -7,7 +7,8 @@ reportEnergyInvestment( gdx, regionSubsetList = NULL, - t = c(seq(2005, 2060, 5), seq(2070, 2110, 10), 2130, 2150) + t = c(seq(2005, 2060, 5), seq(2070, 2110, 10), 2130, 2150), + gdx_ref = NULL ) } \arguments{ @@ -19,6 +20,10 @@ be created.} \item{t}{temporal resolution of the reporting, default: t=c(seq(2005,2060,5),seq(2070,2110,10),2130,2150)} + +\item{gdx_ref}{a GDX object as created by readGDX, or the path to a gdx of the reference run. +It is used to guarantee consistency before 'cm_startyear' for investment variables +using time averaging.} } \value{ MAgPIE object - contains the price variables diff --git a/man/reportPrices.Rd b/man/reportPrices.Rd index 8afa76b5..1b981f0c 100644 --- a/man/reportPrices.Rd +++ b/man/reportPrices.Rd @@ -35,8 +35,9 @@ Read in price information from GDX file, information used in convGDX2MIF.R for the reporting } \examples{ - -\dontrun{reportPrices(gdx)} +\dontrun{ +reportPrices(gdx) +} } \seealso{