Creation of FSDP_process & revision of FSDP_collect

CHANGELOG.md

@@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased]
 
 ### changed
+- **scripts** FSDP_collect handles the health impacts data provided by Marco Springmann, distributed it into the scenario's various reports. It performs a similar operation for global nutrient surplus (which must be calculated on the grid-level and then aggregated).
 - **13_tc** relaxed vm_tau upper limit
 - **scripts** updated FSEC start and output scripts
 - **scripts** update of rds_report to allow gridded intermediate outputs
@@ -21,9 +22,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - **scripts** log files are now written in a subfolder "logs"
 
 ### added
-- **scripts** added output script creating a merged .csv for dietaryIndicators and caloricSupply outputs
+- **scripts** added output script, FSDP_process creating a merged .csv and .gdx for dietaryIndicators and caloricSupply outputs
 - **scripts** added output script creating a set of outputs for Simon Dietz in the FSEC context
-- **scripts** added output script distributing Marco Springmann's health impacts into the respective output scenario report.mif, report.rds, and report_iso.rds
 - **scripts** added output script running MAGICC7 on a MAgPIE scenario
 - **scripts** added output script for gridded crop diversity indices
 - **scripts** added output scripts for FSEC FSDP runs

scripts/output/projects/FSDP_collect.R

@@ -11,14 +11,16 @@
 # ---------------------------------------------------------------
 
 # Version 1.0, Florian Humpenoeder
-#
+
 library(lucode2)
 library(magclass)
 library(gms)
 library(magpiesets)
 library(data.table)
 library(gdx)
 library(quitte)
+library(m4fsdp)
+library(stringr)
 
 ############################# BASIC CONFIGURATION #############################
 if(!exists("source_include")) {
@@ -27,15 +29,8 @@ if(!exists("source_include")) {
   lucode2::readArgs("outputdir")
 }
 ###############################################################################
-cat("\nStarting output generation\n")
-
-reg <- NULL
-iso <- NULL
-grid <- NULL
-missing <- NULL
-
-saveRDS(outputdir,"outputdir.rds")
 
+##########
 #filter out calibration run
 x <- unlist(lapply(strsplit(basename(outputdir),"_"),function(x) x[2]))
 outputdir <- outputdir[which(x %in% c("FSECa", "FSECb", "FSECc", "FSECd", "FSECe"))]
@@ -44,6 +39,71 @@ outputdir <- outputdir[which(x %in% c("FSECa", "FSECb", "FSECc", "FSECd", "FSECe
 x <- unlist(lapply(strsplit(basename(outputdir),"_"),function(x) x[1]))
 if (length(unique(x)) == 1) rev <- unique(x) else stop("version prefix is not identical. Check your selection of runs")
 
+##########
+# Append health impacts reports
+hi_datasets_path <- "/p/projects/magpie/data/FSEC_healthImpactsDatasets_raw"
+hi_datasets      <- list.files(hi_datasets_path)
+hi_versionToUse  <- grep(rev, hi_datasets, value = TRUE)
+
+if (length(hi_versionToUse) == 0) {
+
+    message("In FSDP_collect.R: No corresponding version ID was found within the FSEC health impacts datasets.
+            Using the highest current version.")
+
+    highestVersionNr <- max(as.numeric(str_extract(hi_datasets, "(?<=v)(.*?)(?=_)")))
+    hi_versionToUse <- grep(paste0("v", highestVersionNr), hi_datasets, value = TRUE)
+
+} else if (length(hi_versionToUse) >= 2) {
+    stop("In FSDP_collect.R: More than one health impacts datasets with this scenario's version ID were found.
+          Only one is expected.")
+}
+
+hi_versionToUse_path <- file.path(hi_datasets_path, hi_versionToUse)
+hi_gdx <- suppressWarnings(readGDX(hi_versionToUse_path))
+
+.appendHealthImpacts <- function(.x) {
+    cfg <- gms::loadConfig(file.path(.x, "config.yml"))
+    title <- cfg$title
+
+    message("Appending health impact report: ", title)
+    tryCatch(
+        expr = {
+            appendReportHealthImpacts(healthImpacts_path = hi_gdx,
+                                      scenario = title,
+                                      dir = .x)
+        }, error = function(e) {
+            message("In FSDP_collect.R: Unable to append health impacts!\n", e)
+        }
+    )
+}
+lapply(X = outputdir, FUN = .appendHealthImpacts)
+
+##########
+# Append nutrient surplus reports
+.appendNutrientSurplus <- function(.x) {
+    cfg <- gms::loadConfig(file.path(.x, "config.yml"))
+    title <- cfg$title
+
+    tryCatch(
+        expr = {
+            appendReportNutrientSurplus(scenario = title, dir = .x)
+        }, error = function(e) {
+            message("In FSDP_collect.R: Unable to append the nutrient surplus dataset!\n", e)
+        }
+    )
+}
+lapply(X = outputdir, FUN = .appendNutrientSurplus)
+
+##########
+# Generate output files
+cat("\nStarting output generation\n")
+
+reg <- NULL
+iso <- NULL
+grid <- NULL
+missing <- NULL
+
+saveRDS(outputdir,"outputdir.rds")
 
 for (i in 1:length(outputdir)) {
   print(paste("Processing",outputdir[i]))
@@ -109,7 +169,7 @@ for (i in 1:length(outputdir)) {
     } else missing <- c(missing,outputdir[i])
 
     ## Nitrogen
-    nc_file <- file.path(outputdir[i], paste(cfg$title,"nutrientSurplus_anthropogenic_unaggregated.mz",sep="-"))
+    nc_file <- file.path(outputdir[i], paste(cfg$title,"nutrientSurplus_intensity.mz",sep="-"))
     if(file.exists(nc_file)) {
       a <- read.magpie(nc_file)[,years,]
       getNames(a) <- "nutrientSurplus (kg N per ha)"
@@ -156,7 +216,6 @@ val <- as.data.table(read.quitte(val))
 saveRDS(val,file = file.path("output",paste(rev,"FSDP_validation.rds",sep="_")), version = 2,compress = "xz")
 
 message("Plotting figures ...")
-library(m4fsdp)
 heatmapFSDP(reg,tableType=1,file=file.path("output",paste(rev,"FSDP_heatmap1.jpg",sep="_")))
 heatmapFSDP(reg,tableType=2,file=file.path("output",paste(rev,"FSDP_heatmap2.jpg",sep="_")))
 spatialMapsFSDP(reg,iso,grid,reg2iso,file = file.path("output",paste(rev,"FSDP_spatialMaps.jpg",sep="_")))

scripts/output/projects/FSDP_process.R

@@ -0,0 +1,112 @@
+# |  (C) 2008-2022 Potsdam Institute for Climate Impact Research (PIK)
+# |  authors, and contributors see CITATION.cff file. This file is part
+# |  of MAgPIE and licensed under AGPL-3.0-or-later. Under Section 7 of
+# |  AGPL-3.0, you are granted additional permissions described in the
+# |  MAgPIE License Exception, version 1.0 (see LICENSE file).
+# |  Contact: magpie@pik-potsdam.de
+
+# -------------------------------------------------------------------------------------------------
+# description: Post-processing of FSEC scenarios, specifically generating dietary data
+# comparison script: TRUE
+# -------------------------------------------------------------------------------------------------
+
+# Version 1.00 - Michael Crawford
+
+library(gms)
+library(gdxrrw)
+library(withr)
+
+message("Starting FSDP_process output runscript")
+
+############################# BASIC CONFIGURATION #######################################
+
+if (!exists("source_include")) {
+
+    outputdir <- file.path("output/", list.dirs("output/", full.names = FALSE, recursive = FALSE))
+    lucode2::readArgs("outputdir")
+
+}
+
+#########################################################################################
+
+# -----------------------------------------------------------------------------------------------------------------
+# Merge dietary outputs from FSEC scenarios into two .csv files
+
+message("Merging dietary outputs")
+
+caloricSupply_mergePath     <- file.path("output", "caloricSupply.csv")
+dietaryIndicators_mergePath <- file.path("output", "dietaryIndicators.csv")
+
+if (file.exists(caloricSupply_mergePath) | file.exists(dietaryIndicators_mergePath)) {
+    message("Merge files are already present for these runs. Removing these old files.")
+    file.remove(caloricSupply_mergePath, dietaryIndicators_mergePath)
+}
+
+message("Creating merge files \"caloricSupply.csv\" and \"dietaryIndicators.csv\"")
+file.create(caloricSupply_mergePath, dietaryIndicators_mergePath)
+
+.writeDietaryIndicators <- function(.dir) {
+
+    cfg <- gms::loadConfig(file.path(.dir, "config.yml"))
+    title <- cfg$title
+
+    caloricSupply     <- read.csv(file.path(.dir, paste0(title, "_caloricSupply.csv")), check.names = FALSE)
+    dietaryIndicators <- read.csv(file.path(.dir, paste0(title, "_dietaryIndicators.csv")), check.names = FALSE)
+
+    if (file.info(caloricSupply_mergePath)$size == 0) {
+        # Include the header
+        write.table(caloricSupply, file = caloricSupply_mergePath,
+                    quote = TRUE, sep = ",", row.names = FALSE, col.names = FALSE)
+    } else {
+        # Otherwise simply append
+        write.table(caloricSupply, file = caloricSupply_mergePath,
+                    quote = TRUE, sep = ",", row.names = FALSE, col.names = FALSE,
+                    append = TRUE)
+    }
+
+    if (file.info(dietaryIndicators_mergePath)$size == 0) {
+        # Include the header
+        write.table(dietaryIndicators, file = dietaryIndicators_mergePath,
+                    quote = TRUE, sep = ",", row.names = FALSE, col.names = FALSE)
+    } else {
+        # Otherwise simply append
+        write.table(dietaryIndicators, file = dietaryIndicators_mergePath,
+                    quote = TRUE, sep = ",", row.names = FALSE, col.names = FALSE,
+                    append = TRUE)
+    }
+}
+
+message("Writing dietary datasets")
+
+# Only merge dietary indicators from selected, dietary-related, scenarios
+dietRelatedScenarios <- c("DietEmptyCals", "DietFish", "DietLegumes", "DietMonogastrics", "DietRuminants", "DietVegFruitsNutsSeeds",
+                          "SoilMonogastric",
+                          "LessFoodWaste",
+                          "NoOverweight", "NoUnderweight",
+                          "BAU",
+                          "FSDP",
+                          "SSP1", "SSP3", "SSP4", "SSP5")
+
+outputdir_diets <- lapply(X = dietRelatedScenarios, FUN = function(.x) grep(x = outputdir, pattern = .x, value = TRUE))
+outputdir_diets <- unlist(outputdir_diets)
+
+lapply(X = outputdir_diets, FUN = .writeDietaryIndicators)
+
+
+# -----------------------------------------------------------------------------------------------------------------
+# Produce .gdx files from two .csv files for Marco Springmann
+
+message("Saving dietaryIndicators.csv and caloricSupply.csv as .gdx files")
+
+with_dir(file.path("output"), {
+    gamsScript <- "csv2gdx_dietaryIndicators.gms"
+    gamsScriptLst <- "csv2gdx_dietaryIndicators.lst"
+    file.create(gamsScript)
+    cat("$call csv2gdx dietaryIndicators.csv output=dietaryIndicators.gdx id=dietaryIndicators index=1..6 values=7..11 useHeader=y\n",
+        file = gamsScript, append = T)
+    cat("$call csv2gdx caloricSupply.csv output=caloricSupply.gdx id=caloricSupply index=1,2,3,4 values=5 useHeader=y",
+        file = gamsScript, append = T)
+    gams(gamsScript)
+    file.remove(gamsScript)
+    file.remove(gamsScriptLst)
+})