adapt NGFS to NPi structure, add integratedDamageCosts script

CHANGELOG.md

@@ -11,6 +11,8 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 - **config** regex tests for many parameters
 - **scripts** add script to check fixing of runs to reference run
     [[#1410](https://github.com/remindmodel/remind/pull/1410)]
+- **scripts** add script for cost decomposition of integrated damage runs
+    [[#1445](https://github.com/remindmodel/remind/pull/1445)]
 - **testthat** test and compile all config files
 
 ### fixed

config/scenario_config_coupled_NGFS_v4.csv

@@ -1,25 +1,24 @@
 title;start;copyConfigFrom;oldrun;path_report;qos;magpie_scen;no_ghgprices_land_until;cm_nash_autoconverge_lastrun;path_gdx;path_gdx_ref;path_gdx_bau
 SSP2-Base;NGFS;;SSP2-Base;;;SSP2EU|NPI|nocc_hist|rcp4p5;y2150;;;;
-h_ndc_preUkraine;NGFS;;h_ndc;;;SSP2EU|NDC|nocc_hist|rcp4p5;y2150;2;;;
-h_ndc;NGFS;;h_ndc;;;SSP2EU|NDC|nocc_hist|rcp4p5;y2150;2;;;
 h_cpol;NGFS;;h_cpol;;;SSP2EU|NPI|nocc_hist|rcp4p5;y2150;2;;;
+h_ndc;NGFS;;h_ndc;;;SSP2EU|NDC|nocc_hist|rcp4p5;y2150;2;;;
 o_1p5c;NGFS;;o_1p5c;;;SSP2EU|NDC|nocc_hist|rcp1p9;y2030;2;;;
 o_2c;NGFS;;o_2c;;;SSP2EU|NDC|nocc_hist|rcp2p6;y2030;2;;;
 o_lowdem;NGFS;;o_lowdem;;;SSP2EU|NDC|nocc_hist|rcp1p9|NGFS_o_lowdem;y2030;2;;;
 d_delfrag;NGFS;;d_delfrag;;;SSP2EU|NDC|nocc_hist|rcp2p6;y2030;2;;;
 d_strain;NGFS;;d_strain;;;SSP2EU|NDC|nocc_hist|rcp4p5;y2030;2;;;
-SSP2-Base_d50;d50;SSP2-Base;SSP2-Base;;;SSP2EU|NPI|cc|rcp6p0;;;;;
-SSP2-Base_d95;d95;SSP2-Base_d50;SSP2-Base;;;;;;;;
-SSP2-Base_d50high;d50high;SSP2-Base_d50;SSP2-Base;;;;;;;;
-SSP2-Base_d95high;d95high;SSP2-Base_d50;SSP2-Base;;;;;;;;
+SSP2-Base_d50;0;SSP2-Base;SSP2-Base;;;SSP2EU|NPI|cc|rcp6p0;;;;;
+SSP2-Base_d95;0;SSP2-Base_d50;SSP2-Base;;;;;;;;
+SSP2-Base_d50high;0;SSP2-Base_d50;SSP2-Base;;;;;;;;
+SSP2-Base_d95high;0;SSP2-Base_d50;SSP2-Base;;;;;;;;
+h_cpol_d50;d50;h_cpol;h_cpol;;;SSP2EU|NPI|cc|rcp4p5;;;;;
+h_cpol_d95;d95;h_cpol_d50;h_cpol;;;;;;;;
+h_cpol_d50high;d50high;h_cpol_d50;h_cpol;;;;;;;;
+h_cpol_d95high;d95high;h_cpol_d50;h_cpol;;;;;;;;
 h_ndc_d50;d50;h_ndc;h_ndc;;;SSP2EU|NDC|cc|rcp4p5;;;;;
 h_ndc_d95;d95;h_ndc_d50;h_ndc;;;;;;;;
 h_ndc_d50high;d50high;h_ndc_d50;h_ndc;;;;;;;;
 h_ndc_d95high;d95high;h_ndc_d50;h_ndc;;;;;;;;
-h_cpol_d50;d50;h_cpol;h_cpol;;;SSP2EU|NPI|cc|rcp6p0;;;;;
-h_cpol_d95;d95;h_cpol_d50;h_cpol;;;;;;;;
-h_cpol_d50high;d50high;h_cpol_d50;h_cpol;;;;;;;;
-h_cpol_d95high;d95high;h_cpol_d50;h_cpol;;;;;;;;
 o_1p5c_d50;d50;o_1p5c;o_1p5c;;;SSP2EU|NDC|cc|rcp1p9;;;;;
 o_1p5c_d95;d95;o_1p5c_d50;o_1p5c;;;;;;;;
 o_1p5c_d50high;d50high;o_1p5c_d50;o_1p5c;;;;;;;;

modules/45_carbonprice/NDC/datainput.gms

@@ -79,7 +79,7 @@ $offdelim
   /             ;
 
 *** parameters for selecting NDC years
-Scalar p45_ignoreNDCbefore          "NDC targets before this years are ignored, for example to exclude 2030 targets" /2020/;
+Scalar p45_ignoreNDCbefore          "NDC targets before this years are ignored, for example to exclude 2030 targets" /2024/;
 p45_ignoreNDCbefore = max(p45_ignoreNDCbefore, cm_startyear)
 Scalar p45_ignoreNDCafter           "NDC targets after  this years are ignored, for example to exclude 2050 net zero targets" /2030/;
 Scalar p45_minRatioOfCoverageToMax  "only targets whose coverage is this times p45_bestNDCcoverage are considered. Use 1 for only best." /1.0/;

output.R

@@ -197,7 +197,7 @@ if (comp %in% c("comparison", "export")) {
   }
 } else { # comp = single
   # define slurm class or direct execution
-  outputInteractive <- c("plotIterations", "fixOnRef")
+  outputInteractive <- c("plotIterations", "fixOnRef", "integratedDamageCosts")
   if (! exists("source_include")) {
     # for selected output scripts, only slurm configurations matching these regex are available
     slurmExceptions <- if ("reporting" %in% output) "--mem=[0-9]*[0-9]{3}" else NULL

scripts/output/single/integratedDamageCosts.R

@@ -0,0 +1,230 @@
+library(dplyr)
+library(tidyr)
+library(quitte)
+
+#' @param i_data quitte object
+#' @param scenBaseNoDamage string with scenario name of baseline (NPi) without damages
+#' @param scenNoDamage string with scenario name without damages
+#' @param scenDamage string with scenario name with damages
+#' @param damageName description of damage type (such as high, medium, ...)
+
+computeCostsScen <- function(i_data, scenBaseNoDamage, scenNoDamage, scenDamage, damageName) {
+  i_data <- i_data %>%
+    mutate(model = paste(model)) %>%
+    mutate(scenario = paste(scenario)) %>%
+    mutate(region = paste(region)) %>%
+    mutate(variable = paste(variable)) %>%
+    mutate(unit = paste(unit))
+
+  message("Compute GDP with damages")
+  tmp0_gdpwdamages <- i_data %>%
+    filter(scenario == scenDamage, variable %in% c("GDP|MER", "GDP|PPP", "Damage factor")) %>%
+    select(-unit) %>%
+    pivot_wider(names_from="variable", values_from="value") %>%
+    mutate(`GDP|MER|w/ Macro-Economic Climate Damage` = `GDP|MER` * `Damage factor`) %>%
+    mutate(`GDP|PPP|w/ Macro-Economic Climate Damage` = `GDP|PPP` * `Damage factor`) %>%
+    mutate(`GDP|PPP|including chronic physical risk damage estimate` = `GDP|PPP|w/ Macro-Economic Climate Damage`) %>%
+    mutate(`GDP|MER|including chronic physical risk damage estimate` = `GDP|MER|w/ Macro-Economic Climate Damage`) %>%
+    mutate(unit = "billion US$2005/yr") %>%
+    select(-`GDP|MER`, -`GDP|PPP`, -`Damage factor`) %>%
+    pivot_longer(names_to="variable", values_to="value", cols=c(-"model", -"scenario", -"region", -"unit", -"period")) %>%
+    select(model, scenario, region, variable, unit, period, value)
+  tmp0_gdpwdamages$variable[tmp0_gdpwdamages$variable == "GDP|PPP|including chronic physical risk damage estimate"] <- paste0("GDP|PPP|including ", damageName, " chronic physical risk damage estimate")
+  tmp0_gdpwdamages$variable[tmp0_gdpwdamages$variable == "GDP|MER|including chronic physical risk damage estimate"] <- paste0("GDP|MER|including ", damageName, " chronic physical risk damage estimate")
+
+  message("Compute climate damage costs")
+  tmp1 <- i_data %>%
+    filter(scenario == scenDamage, variable %in% c("GDP|MER", "Damage factor")) %>% #, "Consumption"
+    select(-unit) %>%
+    pivot_wider(names_from="variable", values_from="value") %>%
+    mutate(value = `GDP|MER` * `Damage factor`) %>%
+    mutate(variable = "GDP|MER") %>%
+    mutate(unit = "billion US$2005/yr") %>%
+    select(-`GDP|MER`, -`Damage factor`) %>%
+    select(model, scenario, region, variable, unit, period, value) %>%
+    left_join(
+      i_data %>%
+        filter(scenario == scenBaseNoDamage) %>%
+        filter(variable %in% c("GDP|MER", "Damage factor")) %>% #, "Consumption"
+        select(-unit) %>%
+        pivot_wider(names_from = "variable", values_from = "value") %>%
+        mutate(value = `GDP|MER` * `Damage factor`) %>%
+        mutate(variable = "GDP|MER") %>%
+        mutate(unit = "billion US$2010/yr") %>%
+        select(-scenario, -`GDP|MER`, -`Damage factor`) %>%
+        select(model, region, variable, unit, period, value) %>%
+        rename(value_ref = value),
+      by = c("model", "region", "variable", "unit", "period")
+    )
+
+  # Absolute differences
+  tmp1_diffabs <- tmp1 %>%
+    mutate(value = value - value_ref) %>%
+    select(-value_ref) %>%
+    mutate(variable = ifelse(variable == "GDP|MER", "Macro-Economic Climate Damage|GDP Change", variable)) #%>%
+    #mutate(variable = ifelse(variable == "Consumption", "Macro-Economic Climate Damage|Consumption Change", variable))
+
+  # Recompute global damage values
+  tmp1_diffabs <- rbind(
+    tmp1_diffabs %>%
+      filter(region != "World"),
+    tmp1_diffabs %>%
+      filter(region != "World") %>%
+      group_by(model, scenario, variable, unit, period) %>%
+      summarise(value = sum(value)) %>%
+      ungroup() %>%
+      mutate(region = "World") %>%
+      select(model, scenario, region, variable, unit, period, value)
+  )
+
+#  # Relative differences
+#  tmp1_diffrel <- tmp1 %>%
+#    mutate(value = (value - value_ref)/value_ref * 100) %>%
+#    select(-value_ref) %>%
+#    mutate(variable = ifelse(variable == "GDP|MER", "Macro-Economic Climate Damage|GDP Change", variable)) %>%
+#    #mutate(variable = ifelse(variable == "Consumption", "Macro-Economic Climate Damage|Consumption Change", variable)) %>%
+#    mutate(unit = "%")
+
+  # Get climate policy costs from scenario without damages and write it to scenario with damages
+  tmp2 <- i_data %>%
+    filter(scenario == scenNoDamage,
+           variable %in% c("Policy Cost|Consumption Loss", "Policy Cost|GDP Loss")) %>%
+           mutate(scenario = scenDamage)
+
+  joint <- as_tibble(rbind(tmp1_diffabs, tmp2))
+
+  message("Combine GDP w/ damage, climate damage and policy costs")
+  out <- rbind(tmp1_diffabs, tmp2) %>%
+    pivot_wider(names_from = "variable", values_from = "value") %>%
+    mutate(`Policy Cost and Macro-Economic Climate Damage|GDP Change` = `Macro-Economic Climate Damage|GDP Change`) %>%
+    mutate(`Macro-Economic Climate Damage|GDP Change` = `Macro-Economic Climate Damage|GDP Change` - `Policy Cost|GDP Loss`) %>%
+    #mutate(`Policy Cost and Macro-Economic Climate Damage|Consumption Change` = `Macro-Economic Climate Damage|Consumption Change`) %>%
+    #mutate(`Macro-Economic Climate Damage|Consumption Change` = `Macro-Economic Climate Damage|Consumption Change` + `Policy Cost|Consumption Loss`) %>%
+    pivot_longer(names_to = "variable", values_to = "value",
+      cols = c("Macro-Economic Climate Damage|GDP Change",
+               "Policy Cost|GDP Loss",
+               "Policy Cost|Consumption Loss",
+               "Policy Cost and Macro-Economic Climate Damage|GDP Change"
+    )) %>%
+    select(model, scenario, region, variable, unit, period, value) %>%
+    rbind(tmp0_gdpwdamages)
+  message("Returning damage costs")
+  return(out)
+}
+
+
+#' @param mifdata quitte object or csv or xlsx file name
+#' @param damagestrings string that differentiate runs with integrated damage
+#'    from those without. names(damagestrings) contains the classification
+#'    of the risk such as medium or high.
+#' @param scenBaseNoDamage string with scenario name of baseline (NPi) without damages
+#' @param keepNoDamages boolean whether runs without damages should be kept in the returned data
+calculateDamages <- function(mifdata, damagestrings, scenBaseNoDamage, keepNoDamages = FALSE) {
+  mifdata <- droplevels(quitte::as.quitte(mifdata))
+  if (all(damagestrings == "")) {
+    message("no damagestringe found, returning data")
+    return(mifdata)
+  }
+  if (is.null(names(damagestrings))) names(damagestrings) <- gsub("^_|_$", "", damagestrings)
+  names(damagestrings) <- ifelse(names(damagestrings) == "", damagestrings, names(damagestrings))
+  if (! scenBaseNoDamage %in% levels(mifdata$scenario)) {
+    stop("Baseline without damages is missing in mifdata: ", scenBaseNoDamage)
+  }
+
+  # construct vectors of damage scenarios and corresponding no-damage scenarios
+  scenDamages <- NULL
+  scenNoDamages <- NULL
+  for (d in seq_along(damagestrings)) {
+    scenD <- grep(damagestrings[d], levels(mifdata$scenario), value = TRUE)
+    names(scenD) <- rep(names(damagestrings)[d], length(scenD))
+    scenNoD <- gsub(damagestrings[d], "", scenD)
+    scenDamages <- c(scenDamages, scenD[scenNoD %in% levels(mifdata$scenario)])
+    scenNoDamages <- c(scenNoDamages, scenNoD[scenNoD %in% levels(mifdata$scenario)])
+  }
+  message("Mappings found: ", paste0(scenDamages, " -> ", scenNoDamages, collapse = ", "))
+  if (length(scenDamages) == 0) {
+    message("No scenario containing damagestring=", damagestrings, " found in mifdata.")
+    return(mifdata)
+  }
+
+  message("Recompute global damage factor")
+  # Recompute global damage factor (weighted average)
+  mifdata <- rbind(
+    mifdata %>%
+      filter(!(variable == "Damage factor" & region == "World")),
+    mifdata %>%
+      filter(variable %in% c("GDP|PPP", "Damage factor"), region != "World") %>%
+      select(-unit) %>%
+      pivot_wider(names_from="variable", values_from="value") %>%
+      left_join(
+        mifdata %>%
+          filter(variable == c("GDP|PPP"), region == "World") %>%
+          select(-region, -variable, -unit) %>%
+          rename(gdp_world = value),
+        by = c("model", "scenario", "period")) %>%
+      group_by(model, scenario, period) %>%
+      summarise(value = sum(`Damage factor` * `GDP|PPP`/gdp_world)) %>%
+      ungroup() %>%
+      mutate(region = "World") %>%
+      mutate(variable = "Damage factor") %>%
+      mutate(unit = "1") %>%
+      select(model, scenario, region, variable, unit, period, value)
+  )
+
+  # Recompute global SCC (SCC is uniform across regions, so we just pick USA here)
+  mifdata <- rbind(
+    mifdata %>%
+      filter(!(variable == "Price|Carbon|SCC" & region == "World")),
+    mifdata %>%
+      filter(variable == "Price|Carbon|SCC" & region == "USA") %>%
+      mutate(region = "World")
+  )
+
+  deleteFromMifdata <- c("Policy Cost|GDP Loss", "Policy Cost and Macro-Economic Climate Damage|GDP Change",
+                         "Macro-Economic Climate Damage|GDP Change", "Policy Cost|Consumption Loss")
+  deleteFromCostdata <- c("GDP|PPP", "GDP|MER", "Damage factor", "Macro-Economic Climate Damage|GDP Change",
+                          "Policy Cost and Macro-Economic Climate Damage|GDP Change")
+  # Compute costs
+  cat("Compute costs...\n")
+  returndata <- NULL
+  for (s in seq_along(scenDamages)) {
+    costdata <- computeCostsScen(mifdata, scenBaseNoDamage, scenNoDamages[s], scenDamages[s], names(scenDamages)[s]) %>%
+        filter(! (variable == "Policy Cost|Consumption Loss" & is.na(value))) %>%
+        filter(! variable %in% deleteFromCostdata)
+    dataDamage <- mifdata %>%
+      filter(scenario %in% scenDamages[s]) %>%
+      filter(! variable %in% deleteFromMifdata) %>%
+      rbind(costdata)
+    returndata <- rbind(returndata, dataDamage)
+  }
+  if (keepNoDamages) {
+    returndata <- rbind(filter(mifdata, ! scenario %in% scenDamages), returndata)
+  }
+  return(droplevels(returndata))
+}
+
+nodamagefolder <- "output"
+nodamagescenarios <- c("d_delfrag", "d_strain", "o_2c", "h_cpol", "o_lowdem", "o_1p5c", "h_ndc")
+remmagiter <- 6
+nodamagemifs <- file.path(nodamagefolder, paste0("C_", nodamagescenarios, "-rem-", remmagiter),
+                          paste0("REMIND_generic_C_", nodamagescenarios, "-rem-", remmagiter, ".mif"))
+nodamage <- quitte::as.quitte(nodamagemifs)
+levels(nodamage$scenario) <- gsub("-rem-[0-9]+", "", levels(nodamage$scenario))
+
+runnames <- lucode2::getScenNames(outputdir)
+
+damagestrings <- c(medium = "_d50", high = "_d95high")
+
+for (r in runnames) {
+  message("\n## Deriving mif for ", r)
+  mif <- file.path(outputdir, paste0("REMIND_generic_", r, ".mif"))
+  file.copy(mif, gsub("REMIND_generic", "REMIND_beforedamagecosts", mif))
+  q <- quitte::as.quitte(mif)
+  levels(q$scenario) <- gsub("-rem-[0-9]+", "", levels(q$scenario))
+  damcosts <- calculateDamages(rbind(q, nodamage), damagestrings, "C_h_cpol", keepNoDamages = FALSE)
+  q <- quitte::as.quitte(damcosts)
+  message("Finished this run (should only be a single name): ", paste(levels(q$scenario), collapse = ", "))
+  levels(q$scenario) <- r
+  quitte::write.mif(q, mif)
+  message("See ", gsub("REMIND_generic", "REMIND_beforedamagecosts", mif), " and ", mif, ".")
+}