First version of new TIRF script

scripts/input/climate_assessment_temperatureImpulseResponse.R

@@ -5,68 +5,285 @@
 # |  REMIND License Exception, version 1.0 (see LICENSE file).
 # |  Contact: [email protected]
 # Extract temperature impulse response function (TIRF) from MAGICC. 
-# needs ./magicc/addEmissionPulse.awk
-# AJS 2016
 print("Acquiring temperature impulse response from MAGICC.. ")
 require(dplyr)
+require(tidyr)
+require(purrr)
 require(gdxrrw)
+require(quitte)
+require(lucode2)
+require(yaml)
 igdx(system("dirname $( which gams )", intern = TRUE))
 
-#find MAGICC SCEN file:
+# find base scenario file generated by REMIND's OpenSCM running between iterations
+# it contains emissions already harmonized and infilled
 load(file.path('config.Rdata'),envir = e <- new.env())
-file_scen =  paste0('./magicc/','REMIND_',as.character(e$cfg$title),'.SCEN')
+cfg <- e$cfg
+file_base_scen <- paste0("./climate-assessment-data/ar6_climate_assessment_",as.character(cfg$title),"_harmonized_infilled.csv")
+
+
+# Get output directory to build full paths, some commands require full paths
+outputDir <- getwd()
+
+# Create an extra log file for TIRF calculation
+logFile <- file.path(outputDir, paste0("log_climate_tirf.txt"))
+if (!file.exists(logFile)) {
+    file.create(logFile)
+    createdLogFile <- TRUE
+} else {
+    createdLogFile <- FALSE
+}
+
+# Use a separate temporary climate-assessment folder for the TIRF calculation
+# Not doing so would interfere with the regular non-TIRF MAGICC files, and 
+# will make it tricky to manipulate emissions between iterations
+climateTempDir <- file.path(outputDir, "climate-assessment-data-tirf")
+if (!dir.exists(climateTempDir)) {
+    dir.create(climateTempDir, showWarnings = FALSE)
+    createdClimateTempDir <- TRUE
+} else {
+    createdClimateTempDir <- FALSE
+}
+cat(climateTempDir)
+
+# Set up climate-assessment environment
+# TONN TODO: Most of what's down here is exactly the same as in climate_assessment_run.R
+# It probably makes sense to create a function or a sourceable script in a separate file
+# that can be sourced from both scripts
+
+# START REPEATED PART ========================================================================
+# Get the scenario name
+scenarioName <- getScenNames(outputDir)
+
+# Set default values for the climate assessment config data in case they are not available for backward compatibility
+if (is.null(cfg$climate_assessment_root)) cfg$climate_assessment_root <- "/p/projects/rd3mod/python/climate-assessment/src/"
+if (is.null(cfg$climate_assessment_infiller_db)) cfg$climate_assessment_infiller_db <- "/p/projects/rd3mod/climate-assessment-files/1652361598937-ar6_emissions_vetted_infillerdatabase_10.5281-zenodo.6390768.csv"
+if (is.null(cfg$climate_assessment_magicc_bin)) cfg$climate_assessment_magicc_bin <- "/p/projects/rd3mod/climate-assessment-files/magicc-v7.5.3/bin/magicc"
+if (is.null(cfg$climate_assessment_magicc_prob_file_iteration)) cfg$climate_assessment_magicc_prob_file_iteration <- "/p/projects/rd3mod/climate-assessment-files/parsets/RCP20_50.json"
+
+# Set up climate-assessment related configuration and output files
+climateAssessmentEmi <- file.path(climateTempDir, paste0("ar6_climate_assessment_", scenarioName, ".csv"))
+if (!file.exists(climateAssessmentEmi)) {
+    file.create(climateAssessmentEmi)
+    createdOutputCsv <- TRUE
+} else {
+    createdOutputCsv <- FALSE
+}
+
+# The base name, that climate-assessment uses to derive it's output names
+baseFn <- sub("\\.csv$", "", basename(climateAssessmentEmi))
+
+# Auxiliary input data for climate-assessment and MAGICC7
+infillingDatabaseFile <- normalizePath(cfg$climate_assessment_infiller_db, mustWork = TRUE)
+probabilisticFile <- normalizePath(cfg$climate_assessment_magicc_prob_file_iteration, mustWork = TRUE)
+
+# Extract the location of the climate-assessment scripts and the MAGICC binary from cfg.txt
+scriptsDir <- normalizePath(file.path(cfg$climate_assessment_root, "scripts"))
+magiccBinFile <- normalizePath(file.path(cfg$climate_assessment_magicc_bin))
+magiccWorkersDir <- file.path(normalizePath(climateTempDir), "workers")
+gamsRDir <- Sys.getenv("GAMSROOT")
+if (nchar(gamsRDir) <= 0) {
+    warning("Empty GAMSROOT environment variable")
+}
+
+# Read parameter sets file to ascertain how many parsets there are
+allparsets <- read_yaml(probabilisticFile)
+nparsets <- length(allparsets$configurations)
+
+# Write parameter file with required modifications
+allparsets$configurations[[1]]$nml_allcfgs$endyear <- 2300
+probabilisticFileModified <- normalizePath("probmod.json")
+jsonlite::write_json(allparsets,probabilisticFileModified, pretty=T, auto_unbox=T)
+
+
+
+
+logMsg <- paste0(date(), " =================== SET UP climate-assessment scripts environment ===================\n")
+cat(logMsg)
+capture.output(cat(logMsg), file = logFile, append = TRUE)
+
+# Create working folder for climate-assessment files
+dir.create(magiccWorkersDir, recursive = TRUE, showWarnings = FALSE)
+
+#
+# SET UP MAGICC ENVIRONMENT VARIABLES
+#
+
+# Character vector of all required MAGICC7 environment variables
+magiccEnvs <- c(
+    "MAGICC_EXECUTABLE_7"    = magiccBinFile,    # Specifies the path to the MAGICC executable
+    "MAGICC_WORKER_ROOT_DIR" = magiccWorkersDir, # Directory of magicc workers
+    "MAGICC_WORKER_NUMBER"   = 1                 # TODO: Get this from slurm or nproc
+)
+
+gamsEnvs <- c(
+    "R_GAMS_SYSDIR" = gamsRDir
+)
+
+environmentVariables <- c(magiccEnvs, gamsEnvs)
+
+# Check if all necessary environment variables are set
+alreadySet <- lapply(Sys.getenv(names(environmentVariables)), nchar) > 0
+# alreadySet[]=F
+# Only set those environment variables that are not already set
+if (any(!alreadySet)) do.call(Sys.setenv, as.list(environmentVariables[!alreadySet]))
+
+# END REPEATED PART ========================================================================
+
 
 # define years and pulse size
 scan_vals_pulse = c(0,1) # pulse size in GtC. Keep 0 in here as the baseline against which the pulse run is compared to
-scan_vals_years = c(2010,2020,2030,2040,2050,2060,2070,2080,2090,2100,2110,2130)
+scan_vals_years = c(2020,2030,2040,2050,2060,2070,2080,2090,2100,2110,2130)
+
+# Read base emissions scenario file
+mif_base_scen <- read.quitte(file_base_scen)
 
-# backup default scenario file
-file.copy(file_scen,paste0(file_scen,'_backup'),overwrite = T)
+# TONN START CHANGE
+# Add emissions data up to 2300, assuming constant emissions of all species after 2100
+mif_base_scen <- bind_rows(mif_base_scen,
+  lapply(2101:2200,function(x){
+    mif_base_scen %>%
+    filter(period == 2100) %>% mutate(period = x)
+  })) %>%
+  arrange(variable,period) 
 
+# The variable which to apply the pulse
+# Keep flexible, as we may want to do this with other gases as well
+# Here we assume the pulse comes from energy emissions. AFOLU emissions
+# imply deforestation in some SCMs and would lead to spurious regrowth
+emis_var_name <- "AR6 climate diagnostics|Infilled|Emissions|CO2|Energy and Industrial Processes"
 
-getTemperatureMagicc = function(file = "./magicc/DAT_SURFACE_TEMP.OUT"){
-  x = read.table(file, skip = 19,header = T)
-  x = x[,c(1,2)]
-  names(x) = c('period','value')
-  x$period = as.integer(as.character(x$period))
-  # Get relevant years
-  years <- x[x$period >= 2005 & x$period <= 2300,1]
-  return(x[x$period %in% years,])
+# Build a quitte with scenarios that scan the range of pulse sizes and years set above
+# This is pretty fast, so opting for a readable loop instead of an apply
+separator_scen <- "---"
+mif_allpulses_scen <- tibble()
+for (val_pulse in scan_vals_pulse) {
+  for (val_year in scan_vals_years) {
+    # val_pulse <- 1
+    # val_year <-2020
+
+    mif_pulse <- mif_base_scen %>%
+      mutate(value = case_when(
+        variable == emis_var_name & period == val_year ~ value + val_pulse*1e3*3.66,#GtC to MtCO2
+        .default = value
+        ),
+        scenario = paste0("P",separator_scen,val_pulse,separator_scen,val_year)) 
+
+    mif_allpulses_scen <- bind_rows(mif_allpulses_scen, mif_pulse)
+
+  }
 }
+# TONN END CHANGE
+
+# Write the emissions file
+file_allpulses_scen <- paste0(normalizePath(climateTempDir),"/allpulses.xlsx")
+write.IAMCxlsx(mif_allpulses_scen,file_allpulses_scen)
+
+# Where we expect the climate output to be
+file_allpulses_climate <- paste0(normalizePath(climateTempDir),"/allpulses_IAMC_climateassessment.xlsx")
+
+# BUILD climate-assessment RUN COMMAND
+runClimateEmulatorCmd <- paste(
+    "python", file.path(scriptsDir, "run_clim.py"),
+    # normalizePath(file.path(climateTempDir, paste0(baseFn, "_harmonized_infilled.csv"))),
+    # normalizePath(file_base_scen),
+    # normalizePath("./climate-assessment-data-tirf/dummy.csv"),
+    file_allpulses_scen,
+    climateTempDir,
+    "--year-filter-last", 2305, # REQUIRES CHANGES FROM https://github.com/gabriel-abrahao/climate-assessment/tree/yearfilter
+    "--num-cfgs", nparsets,
+    "--scenario-batch-size", 1,
+    "--probabilistic-file", probabilisticFileModified
+)
+
+# Get conda environment folder
+condaDir <- file.path(dirname(dirname(outputDir)), "ca_env/")
+
+# Command to activate the conda environment
+condaCmd <- paste0("module load conda/2023.09; source activate ",condaDir,";")
+
+logMsg <- paste0(
+    date(), "  CLIMATE-ASSESSMENT ENVIRONMENT:\n",
+    "  climateTempDir        = '", climateTempDir, "' exists? ", dir.exists(climateTempDir), "\n",
+    "  baseFn                = '", baseFn, "'\n",
+    "  probabilisticFile     = '", probabilisticFile, "' exists? ", file.exists(probabilisticFile), "\n",
+    "  probabilisticFileModified     = '", probabilisticFileModified, "' exists? ", file.exists(probabilisticFileModified), "\n",
+    "  infillingDatabaseFile = '", infillingDatabaseFile, "' exists? ", file.exists(infillingDatabaseFile), "\n",
+    "  scriptsDir            = '", scriptsDir, "' exists? ", dir.exists(scriptsDir), "\n",
+    "  magiccBinFile         = '", magiccBinFile, "' exists? ", file.exists(magiccBinFile), "\n",
+    "  magiccWorkersDir      = '", magiccWorkersDir, "' exists? ", dir.exists(magiccWorkersDir), "\n",
+    "  condaDir      = '", condaDir, "' exists? ", dir.exists(condaDir), "\n\n",    
+    "  ENVIRONMENT VARIABLES:\n",
+    "  MAGICC_EXECUTABLE_7    = ", Sys.getenv("MAGICC_EXECUTABLE_7"), "\n",
+    "  MAGICC_WORKER_ROOT_DIR = ", Sys.getenv("MAGICC_WORKER_ROOT_DIR"), "\n",
+    "  MAGICC_WORKER_NUMBER   = ", Sys.getenv("MAGICC_WORKER_NUMBER"), "\n",
+    "  R_GAMS_SYSDIR          = ", Sys.getenv("R_GAMS_SYSDIR"), "\n",
+    date(), " =================== CONDA ACTIVATION COMMAND ===================\n",    
+    condaCmd, "'\n",
+    date(), " =================== SKIP climate-assessment infilling & harmonization ===================\n",    
+    date(), " =================== RUN climate-assessment model runs ===================\n",
+    runClimateEmulatorCmd, "'\n"
+)
+
+cat(logMsg)
+capture.output(cat(logMsg), file = logFile, append = TRUE)
+
+# Start actual runs ====================================================
+timeStartEmulation <- Sys.time()
+system(paste0(condaCmd,runClimateEmulatorCmd))
+timeStopEmulation <- Sys.time()
+
+# Actual runs done, read the output, already filtering what we need
+temperature_var_name <- "AR6 climate diagnostics|Surface Temperature (GSAT)|MAGICCv7.5.3|50.0th Percentile"
+mif_allpulses_climate <- read.quitte(file_allpulses_climate) %>%
+  filter(variable == temperature_var_name) %>%
+  filter(between(period,2020,2300))
+
+
+tirf <- mif_allpulses_climate %>%
+  select(scenario,period,value) %>%
+  separate(scenario,c("dummy","size_pulse","year_pulse"), sep = separator_scen) %>%
+  # select(-dummy) %>%
+  select(year_pulse, size_pulse, period, value) %>%
+  mutate(
+    year_pulse = as.numeric(year_pulse),
+    size_pulse = as.numeric(size_pulse)
+  ) %>%
+  as.data.frame 
 
-#initialize
-scan_params = list('year_pulse'=2020,"size_pulse"=1)
-
-
-# loop over pulse experiments
-tirf = 
-  do.call(rbind,lapply(scan_vals_years, function(y){
-    scan_params[['year_pulse']] = y  
-    tirfYear = do.call(rbind,lapply(scan_vals_pulse, function(i) {
-      # manipulate scenario file; add pulse
-      scan_params[['size_pulse']] = i
-      opts = paste0(paste0('-v ',paste0(names(scan_params),'=',scan_params)),collapse = ' ')
-      cmd = paste0("awk -f ./magicc/addEmissionPulse.awk ",opts,' ',file_scen,' > ','tmp',' && ','mv tmp ',file_scen)
-      system(cmd)
-      #run MAGICC
-      system('Rscript run_magicc.R')
-      #reset scenario file
-      file.copy(paste0(file_scen,'_backup'),file_scen,overwrite = T)
-      #read results
-      tirfYearPulsesize = getTemperatureMagicc()
-      tirfYearPulsesize = merge(as.data.frame(scan_params),tirfYearPulsesize) 
-      #return
-      tirfYearPulsesize
-    }))
-    tirfYear
-  }))
 
 #calculate difference to baseline to get TIRF, normalize to 1 GtCO2eq emission.
-tirf = tirf %>% 
+tirf <- tirf %>% 
   group_by(period,year_pulse) %>% 
   summarize(tirf = (value[size_pulse != 0] - value[size_pulse == 0])/size_pulse[size_pulse!=0]/(44/12)) %>% 
   ungroup()
 
+# Extend period back to 2010 to all year_pulses, assuming tirf 0
+tirf <- tirf %>%
+  nest(-year_pulse) %>%
+  mutate(data = map(data,function(x){
+    rbind(data.frame(
+      period = 2010:2019, tirf = 0
+    ),x) %>%
+    arrange(period)
+  })) %>%
+  unnest()
+
+# Assume tirf has the same shape in year_pulse 2010 as in 2020
+tirf <- rbind(
+  tirf %>%
+    filter(year_pulse == 2020) %>%
+    mutate(year_pulse = year_pulse - 10, period = period - 10),
+  tirf
+) %>%
+filter(period >= 2010)
+
+
+saveRDS(tirf,"newtirf.rds")
+tirf <- readRDS("newtirf.rds")
+write.table(tirf,"newtirf.csv",sep=";",row.names=F)
+
+
 #FIXME the result for 2150 is just zero, don't know why. work around by assuming the TIRF in 2150 is equal to the one in 2130. From 2150 to 2250, assume the same. 
 tirf = rbind(tirf,
   tirf %>% 
@@ -102,7 +319,6 @@ oupt = oupt %>%
   filter(tall<=2250,tall>=2005)  %>%
   select(tall,tall1,tirf)
 
-
 writeToGdx = function(file="pm_magicc_temperatureImpulseResponse",df){
   df$tall = factor(df$tall)
   df$tall1 = factor(df$tall1)
@@ -117,6 +333,3 @@ writeToGdx = function(file="pm_magicc_temperatureImpulseResponse",df){
 # write to GDX:
 writeToGdx('pm_magicc_temperatureImpulseResponse',oupt)
 print("...done.")
-
-
-