remind: The REMIND R Package

Documented in reportTechnology

#' Read in GDX and calculate technology information, used in convGDX2MIF.R for the reporting
#'
#' Read in technology 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 output a magpie object containing all needed variables generated by other report*.R functions
#' @return MAgPIE object - contains the technology variables
#' @author Michaja Pehl, Lavinia Baumstark
#' @seealso \code{\link{convGDX2MIF}}
#' @examples
#'
#' \dontrun{reportTechnology(gdx)}
#'
#' @export
#' @importFrom gdx readGDX
#' @importFrom magclass getYears getSets collapseNames new.magpie getRegions getSets<- mbind setNames getNames
#' @importFrom luscale speed_aggregate
reportTechnology <- function(gdx,output=NULL,regionSubsetList=NULL) {
  if(is.null(output)){
    output <- mbind(
      reportSE(gdx, regionSubsetList),
      reportFE(gdx, regionSubsetList),
      reportEmi(gdx, regionSubsetList)
    )
  }

  ## Check realisations
  module2realisation <- readGDX(gdx, "module2realisation", react = "silent")
  tran_mod = module2realisation[module2realisation$modules == "transport", 2]
  if ("CCU" %in% module2realisation[,1]) {
    CCU_mod = module2realisation[module2realisation$modules == "CCU", 2]
  } else {
    CCU_mod <- "off"
  }
  
  CDR_mod = module2realisation[module2realisation$modules == "CDR", 2]
  
  sety       <- readGDX(gdx,c("entySe","sety"),format="first_found")
  # calculate maximal temporal resolution
  p_dataeta    <- readGDX(gdx,name=c("pm_dataeta","p_dataeta"),format = "first_found")
  p_eta_conv   <- readGDX(gdx,name=c("pm_eta_conv","p_eta_conv"),format = "first_found")
  pm_inco0_t   <- readGDX(gdx,name=c("pm_inco0_t","p_inco0_t"),format = "first_found")
  v_investcost <- readGDX(gdx,name=c("vm_costTeCapital","v_costTeCapital","v_investcost"),field = "l",format = "first_found")

  y <- Reduce(intersect,list(getYears(p_dataeta),getYears(p_eta_conv),getYears(v_investcost)))
  p_dataeta    <- p_dataeta[,y,]
  p_eta_conv   <- p_eta_conv[,y,]
  pm_inco0_t   <- pm_inco0_t[,y,]
  v_investcost <- v_investcost[,y,]

  pm_data      <- collapseNames(readGDX(gdx,name=c("pm_data"),format = "first_found")[,,"inco0"])
  inco0 <- new.magpie(getRegions(pm_data),y,magclass::getNames(pm_data))
  for(i in y){
      inco0[,i,] <- pm_data
  }
  # use pm_inco0_t for v_investcost for REMIND 1.6 where it was 0 for all non-learning technologies
  v_investcost[v_investcost == 0] <- pm_inco0_t[v_investcost == 0]

  pm_data      <- collapseNames(readGDX(gdx,name=c("pm_data"),format = "first_found")[,,"omf"])
  omf <- new.magpie(getRegions(pm_data),y,magclass::getNames(pm_data))
  for(i in y){
      omf[,i,] <- pm_data
  }

  pm_data      <- collapseNames(readGDX(gdx,name=c("pm_data"),format = "first_found")[,,"omv"])
  omv <- new.magpie(getRegions(pm_data),y,magclass::getNames(pm_data))
  for(i in y){
    omv[,i,] <- pm_data
  }

  pm_data      <- collapseNames(readGDX(gdx,name=c("pm_data"),format = "first_found")[,,"lifetime"])
  lifetime <- new.magpie(getRegions(pm_data),y,magclass::getNames(pm_data))
  for(i in y){
    lifetime[,i,] <- pm_data
  }
  getSets(v_investcost)[1] <- getSets(inco0)[1]
  getSets(v_investcost)[2] <- getSets(inco0)[2]
  getSets(v_investcost)[3] <- getSets(inco0)[3]
  getSets(p_dataeta)[1]    <- getSets(inco0)[1]
  getSets(p_dataeta)[2]    <- getSets(inco0)[2]
  getSets(p_dataeta)[3]    <- getSets(inco0)[3]
  getSets(p_eta_conv)[1]   <- getSets(inco0)[1]
  getSets(p_eta_conv)[2]   <- getSets(inco0)[2]
  getSets(p_eta_conv)[3]   <- getSets(inco0)[3]

  ############ build reporting #####################

  techmap <- c(
      "bioigccc"="Electricity|Biomass|IGCCC|w/ CCS",
      "bioigcc"="Electricity|Biomass|IGCC|w/o CCS",
      "biochp"="Electricity|Biomass|CHP|w/o CCS",
      "igccc"="Electricity|Coal|IGCCC|w/ CCS",
      "igcc"="Electricity|Coal|IGCC|w/o CCS",
      "pcc"="Electricity|Coal|PCC|w/ CCS",
      "pc"="Electricity|Coal|PC|w/o CCS",
      "coalchp"="Electricity|Coal|CHP|w/o CCS",
      "ngccc"="Electricity|Gas|CCC|w/ CCS",
      "ngcc"="Electricity|Gas|CC|w/o CCS",
      "gaschp"="Electricity|Gas|CHP|w/o CCS",
      "ngt"="Electricity|Gas|GT",
      "dot"="Electricity|Oil|DOT",
      "geohdr"="Electricity|Geothermal",
      "hydro"="Electricity|Hydro",
      "tnrs"="Electricity|Nuclear",
      "spv"="Electricity|Solar|PV",
      "csp"="Electricity|Solar|CSP",
      "wind"="Electricity|Wind",
      "storspv"="Electricity|Storage|Battery|For PV",
      "storcsp"="Electricity|Storage|Battery|For CSP",
      "storwind"="Electricity|Storage|Battery|For Wind",
      "biogas"="Gases|Biomass|w/o CCS",
      "coalgas"="Gases|Coal|w/o CCS",
      "bioh2c"="Hydrogen|Biomass|w/ CCS",
      "bioh2"="Hydrogen|Biomass|w/o CCS",
      "coalh2c"="Hydrogen|Coal|w/ CCS",
      "coalh2"="Hydrogen|Coal|w/o CCS",
      "elh2"="Hydrogen|Electricity",
      "gash2c"="Hydrogen|Gas|w/ CCS",
      "gash2"="Hydrogen|Gas|w/o CCS",
      "bioftcrec"="Liquids|Biomass|Biofuel|BioFTRC|w/ CCS",
      "bioftrec"="Liquids|Biomass|Biofuel|BioFTR|w/o CCS",
      "bioethl"="Liquids|Biomass|Biofuel|Ethanol|Cellulosic|w/o CCS",
      "bioeths"="Liquids|Biomass|Biofuel|Ethanol|Conventional|w/o CCS",
      "biodiesel"="Liquids|Biomass|Biofuel|Biodiesel|w/o CCS",
      "coalftcrec"="Liquids|Coal|w/ CCS",
      "coalftrec"="Liquids|Coal|w/o CCS")
  
  if(tran_mod == "complex"){
      carmap <- c(
          "apCarPeT"="Transport|Pass|Road|LDV|ICE",
          "apCarElT"="Transport|Pass|Road|LDV|EV",
          "apCarH2T"="Transport|Pass|Road|LDV|H2")
  }else{
      carmap <- c()
  }
  
  if (CCU_mod == "on") {
    techmap <- append(techmap, c("MeOH" = "Liquids|Hydrogen",
                                 "h22ch4" = "Gases|Hydrogen"))
  }
  
  if (CDR_mod != "off") {
    cdrmap <- c("dac" = "DAC",
                "ccsinje" = "CO2 Storage")
  } else {
    cdrmap <- c()
  }
  
  if (("seliq" %in% sety) || ("seliqbio" %in% sety)) {
      techmap[["refliq"]] <- "Liquids|Oil"
  }else{
      techmap[["refdip"]] <- "Liquids|Oil"
  }

  bar_and <- function(str){
      ## prepend pipe if not empty
      ifelse(str == "", str, paste0("|", str))
  }

  report_str <- function(tech, category="", unit="", predicate="Tech"){
      ## Construct a reporting string of the form predicate|tech|category (unit)
      if(unit != "")
          unit <- paste0(" (", unit, ")")
      paste0(predicate, bar_and(tech), bar_and(category), unit)
  }

  bind_element <- function(data, parameters, key, category, unit, factor=1., map=techmap){
      ## Add an element, described by *key*, from a *parameters* (magpie) dataset to *data*
      ## The new key in *data* is given by the report string.
      rep_str <- report_str(map[[key]], category, unit)
      return(mbind(data, setNames(parameters[,,key] * factor, rep_str)))
  }

  bind_category <- function(data, parameters, category, unit, factor=1., map=techmap){
      ## Like bind_element, for all keys in *map*
      for(key in names(map)){
          data <- bind_element(data, parameters, key, category, unit, factor, map)
      }
      return(data)
  }

  get_global_mapping <- function(category, unit, map=techmap){
      ## When calculating global averages, the regions have to be weighted by their respective
      ## contributions. Conventionally, the SE value of the respective energy technology is used.
      ## For cars (technically a Energy Service), we use FEs.
      int2ext <- c()
      if(all(map %in% techmap)){
          for(label in techmap){
              int2ext[[report_str(label, category, unit)]] <- report_str(label, unit="EJ/yr", predicate="SE")
          }
          ## storage needs a special mapping
          int2ext[[report_str("Electricity|Storage|Battery|For PV", category, unit)]] <- report_str("Electricity|Solar|PV", unit="EJ/yr", predicate="SE")
          int2ext[[report_str("Electricity|Storage|Battery|For CSP", category, unit)]] <- report_str("Electricity|Solar|CSP", unit="EJ/yr", predicate="SE")
          int2ext[[report_str("Electricity|Storage|Battery|For Wind", category, unit)]] <- report_str("Electricity|Wind", unit="EJ/yr", predicate="SE")

      }else if(all(map %in% carmap)){
          ## cars need a special mapping, too
          ## for global avgs we use FEs as weights
          int2ext[[report_str("Transport|Pass|Road|LDV|ICE", category, unit)]] <- report_str("Transport|Pass|Road|LDV|Liquids", unit="EJ/yr", predicate="FE")
          int2ext[[report_str("Transport|Pass|Road|LDV|EV", category, unit)]] <- report_str("Transport|Pass|Road|LDV|Electricity", unit="EJ/yr", predicate="FE")
          int2ext[[report_str("Transport|Pass|Road|LDV|H2", category, unit)]] <- report_str("Transport|Pass|Road|LDV|Hydrogen", unit="EJ/yr", predicate="FE")
      } else if(all(map %in% cdrmap)){
        # CDR technologies need special mapping
        # for global avgs we use CO2 flows as weights
        int2ext[[report_str("DAC", category, unit)]] <- report_str("DAC", unit="Mt CO2/yr", predicate="Emi|CO2")
        int2ext[[report_str("CO2 Storage", category, unit)]] <- report_str("Carbon Capture and Storage", unit="Mt CO2/yr", predicate="Emi|CO2")
      }
      return(int2ext)
  }

  tmp <- NULL

  ### capital costs ###

  category <- "Capital Costs"
  unit <- "US$2005/kW"
  factor <- 1000.

  tmp <- bind_category(tmp, v_investcost, category, unit, factor, techmap)
  int2ext <- get_global_mapping(category, unit, techmap)

  if(tran_mod == "complex"){
      unit <- "US$2005/veh"
      tmp <- bind_category(tmp, v_investcost, category, unit, factor, carmap)
      int2ext <- c(int2ext, get_global_mapping(category, unit, carmap))
  }
  
  if(CDR_mod != "off"){
    unit <- "US$2005/tCO2 yr"
    factor <- 1000/3.6
    tmp <- bind_category(tmp, v_investcost, category, unit, factor, cdrmap)
    int2ext <- c(int2ext, get_global_mapping(category, unit, cdrmap))
  }

  ### efficiency ###
  ## efficiency variables can be found in both p_dataeta and p_eta_conv, so do it one-by-one
  ## for cars, efficiencies are not given

  category <- "Efficiency"
  unit <- "%"
  factor <- 100.

  in_dataeta <- c("bioigccc", "bioigcc", "igccc", "igcc", "pc", "ngccc", "ngcc", "ngt")

  for(key in names(techmap)){
      if(key %in% in_dataeta){
          tmp <- bind_element(tmp, p_dataeta, key, category, unit, factor)
      }else{
          tmp <- bind_element(tmp, p_eta_conv, key, category, unit, factor)
      }
  }

  int2ext <- c(int2ext, get_global_mapping(category, unit, techmap))


  ### lifetime ###

  category <- "Lifetime"
  unit <- "years"

  tmp <- bind_category(tmp, lifetime, category, unit, 1., techmap)
  int2ext <- c(int2ext, get_global_mapping(category, unit, techmap))
  if(tran_mod == "complex"){
      tmp <- bind_category(tmp, lifetime, category, unit, 1., carmap)
      int2ext <- c(int2ext, get_global_mapping(category, unit, carmap))
  }
  
  if(CDR_mod != "off"){
    tmp <- bind_category(tmp, lifetime, category, unit, 1., cdrmap)
    int2ext <- c(int2ext, get_global_mapping(category, unit, cdrmap))
  }
  

  
  ### o&m fix costs ###
  category <- "OM Cost|fixed"
  unit <- "US$2005/kW/yr"
  tmp <- bind_category(tmp, omf * v_investcost, category, unit, 1000.)
  int2ext <- c(int2ext, get_global_mapping(category, unit, techmap))

  if(tran_mod == "complex"){
      ## op costs for cars ###
      category <- "Op Costs"
      unit <- "US$2005/veh/yr"
      tmp <- bind_category(tmp, omf * v_investcost, category, unit, 1000., carmap)
      int2ext <- c(int2ext, get_global_mapping(category, unit, carmap))
  }
  

  if(CDR_mod != "off"){
    ## op costs for CDR technologies ###
    category <- "OM Cost|fixed"
    unit <- "US$2005/tCO2 yr"
    tmp <- bind_category(tmp, omf * v_investcost, category, unit, 1000/3.66, cdrmap)
    int2ext <- c(int2ext, get_global_mapping(category, unit, cdrmap))
  }

  ### o&m variable costs ###
  category <- "OM Cost|variable"
  unit <- "US$2005/GJ"
  tmp <- bind_category(tmp, omv, category, unit, 1000./31.7098)
  int2ext <- c(int2ext, get_global_mapping(category, unit, techmap))


  ### write to output ###
  output[is.na(output)] <- 0  # substitute na by 0

  # add global values
  map <- data.frame(region=getRegions(tmp),world="GLO",stringsAsFactors=FALSE)
  y <- Reduce(intersect,list(getYears(tmp),getYears(output)))
  tmp <- tmp[,y,]
  tmp_GLO <- new.magpie("GLO",getYears(tmp),magclass::getNames(tmp),fill=0)

  for (i2e in names(int2ext)){
    tmp_GLO["GLO",,i2e] <- speed_aggregate(tmp[,,i2e],map,weight=output[map$region,,int2ext[[i2e]]])
  }
  tmp <- mbind(tmp,tmp_GLO)

  # add other region aggregations
  if (!is.null(regionSubsetList)){
    tmp_RegAgg <- new.magpie(names(regionSubsetList),getYears(tmp),magclass::getNames(tmp),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 <- speed_aggregate(tmp[regionSubsetList[[region]],,i2e],map,weight=output[regionSubsetList[[region]],,int2ext[i2e]])
        getRegions(result) <- region
        for(t in getYears(tmp)){
          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 <- mbind(tmp,tmp_RegAgg)
  }

  tmp[is.na(tmp)] <- 0  # tmp is NA if weight is zero for all regions within the GLO or the specific region aggregation. Therefore, we replace all NAs with zeros.

  return(tmp)
}