added most recent version of the script that adds combined and divisi…

code/cobine_so2_so4.R

@@ -0,0 +1,113 @@
+# ------------------------------------------------------------------------------
+# Program Name: combine_so2_so4.R
+# Authors: Harrison Suchyta
+# Date Last Modified: September 15, 2022
+# Program Purpose: this script combines the high SO2 and SO4 at height values in order to compare
+# that value with the high-SO4-at-height experiment
+# Input Files: ~Emissions-MIP/input/
+# Output Files: ~Emissions-MIP/output/
+# TODO: ADD GISS functionality
+# ------------------------------------------------------------------------------
+
+# Load required libraries
+library(dplyr)
+library(tidyr)
+library(purrr)
+library(ggplot2)
+library(gridExtra)
+library(grid)
+
+
+
+
+#set emissions directory
+emi_dir <- paste0('C:/Users/such559/Documents/Emissions-MIP_Data')
+#-------------------------------------------------------------------------------
+#Read in variable, region, and experiment names and sort them into their own lists
+var_master_list <- read.csv(file = paste0(emi_dir, '/input/var_master_list.csv'), fileEncoding="UTF-8-BOM", stringsAsFactors = FALSE)
+master_vars <- var_master_list$vars
+master_com_var <- var_master_list$com_vars
+master_region <- var_master_list$reg
+master_exp <- var_master_list$exper
+#remove blank strings
+master_vars <- master_vars[master_vars != ""]
+master_com_var<- master_com_var[master_com_var != ""]
+master_region <- master_region[master_region != ""]
+master_exp <- master_exp[master_exp != ""]
+#-------------------------------------------------------------------------------
+#create a function to accumulate data from csv files
+data_accumulation <- function(emi_dir, reg_name, exper){
+
+    setwd(paste0(emi_dir,'/input/', reg_name,'/', exper, '/diff'))
+
+    # Read in csv files and bind into single data frame
+    target_filename <- list.files(getwd(), "*.csv")
+    regional_data <- rbind(map(target_filename, read.csv))
+    regional_data <- lapply(regional_data, function(x) {x["unit"] <- NULL; x})
+    regional_data <- bind_rows(regional_data)
+
+    # Extract model from file names (fifth segment) and bind to experiment data frame
+    models <- sapply(strsplit(target_filename, "[-.]+"),function(x) x[5])
+    rep_models <- rep(models, each = 5) # four years
+    regional_data$model <- rep_models
+
+    # Convert SO2 volume mixing ratio to mass mixing ratio by multiplying by molar
+    # mass of SO2 and dividing by molar mass of air, invert sign of forcing variables
+    # to be consistent with convention (i.e. positive value denotes a heating effect),
+    # then take the average over all years for each variable and calculate std dev
+    regional_data_summary <- regional_data %>%
+        within(value <- ifelse(variable == "so2", 64.066 / 28.96, 1) * value) %>%
+        within(value <- ifelse(variable %in% c("rlut", "rsut", "rlutcs", "rsutcs"), -1, 1) * value) %>%
+        within(value <- ifelse(variable %in% c("wetbc", "wetso2", "wetso4") & model == "CESM2", -1, 1) * value) %>%
+        within(value <- ifelse(variable == "dms", 62.13 / 28.96, 1) * value) %>%
+        dplyr::group_by(variable, model) %>%
+        dplyr::summarise(regional_data = mean(value), regional_data_sd = sd(value))
+
+    return(regional_data_summary)
+}
+#-------------------------------------------------------------------------------
+#Make a function to combine experiment outputs
+combine_expers <- function(region_name,curr_names,summary_list){
+
+    #select both dataframes associated with the current region
+    curr_reg <- curr_names[grepl(region_name,curr_names) ==TRUE]
+
+    curr_df1 <- summary_list[[curr_reg[1]]]
+    curr_df2 <- summary_list[[curr_reg[2]]]
+
+    #add the values and standard deviations
+    sum_df <- curr_df1 %>%
+        left_join(curr_df2,by = c('variable','model')) %>%
+        mutate(regional_data = regional_data.x + regional_data.y,
+               regional_data_sd = sqrt(regional_data_sd.x^2 + regional_data_sd.y^2),
+               experiment = 'so2_so4_sum') %>%
+        subset(select = -c(regional_data.x,regional_data.y,regional_data_sd.x,regional_data_sd.y))
+
+    return(sum_df)
+
+}
+#-------------------------------------------------------------------------------
+#read in high-So2 and SO4 at height values
+
+experiments <- c('high-so4','so2-at-height')
+
+summary_data_list <- list()
+
+#loop through each region and extract values. Filter to just E3SM
+for (exper in experiments){
+    for (reg_name in master_region){
+        summary_data_list[[paste0(exper,'_summary_',reg_name)]] <- assign(paste0(exper,'_summary_',reg_name),data_accumulation(emi_dir,reg_name,exper)) %>%
+            filter(model %in% c('E3SM','GISS'))
+    }
+}
+
+#create a list of the names in summary_data_list
+data_names <- names(summary_data_list)
+
+summed_so2_so4 <- list()
+#loop through each region and add the results together
+for (reg_name in master_region){
+
+    summed_so2_so4[[paste0('so4_so2_sum_',reg_name)]] <- assign(paste0('so4_so2_sum_',reg_name),combine_expers(reg_name,data_names,summary_data_list))
+
+}

code/excluded_data.csv

@@ -1 +1,6 @@
 Model,Scenario,Variable
+,bc-no-season,
+,high-so4,
+,no-so4,
+,so2-no-season,
+,so2-at-height,

code/summary_plots_diff.R

@@ -1,7 +1,7 @@
 # ------------------------------------------------------------------------------
 # Program Name: summary_plots_diff.R
 # Authors: Hamza Ahsan, Harrison Suchyta
-# Date Last Modified: January 14, 2022
+# Date Last Modified: April 4, 2022
 # Program Purpose: Produces summary plots of the difference between the
 # perturbations and the reference case averaged over all years
 # Input Files: ~Emissions-MIP/input/
@@ -26,8 +26,8 @@ setwd(paste0(emi_dir))
 sort_by <- 'region'
 
 #determines which region or experiment is sorted by
-region <- 'arctic'
-exper <- 'bc-no-season'
+region <- 'global'
+exper <- 'high-SO4-at-height'
 
 # Specify what you are sorting by and either the region (i.e., global, land, sea, arctic, NH-land, NH-sea, SH-land, SH-sea) or experiment (i.e., bc-no-season, high-so4, no-so4, reference, so2-at-height, so2-no-season)
 #The command line would look like: rscript <rscript>.r <"experiment" or "region"> <specific experiment or region you are sorting by>
@@ -95,14 +95,15 @@ cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#920000",
                "#F0E442", "#0072B2", "#D55E00", "#CC79A7", "#490092",
                "#117733")
 
-model_colors <- c(CESM1 = cbPalette[1], E3SM = cbPalette[2], GISS = cbPalette[3],
-                  CESM2 = cbPalette[4], MIROC = cbPalette[5], NorESM2 = cbPalette[6],
-                  GFDL = cbPalette[7], OsloCTM3 = cbPalette[8], UKESM = cbPalette[9],
-                  GEOS = cbPalette[10], CAM5 = cbPalette[11])
+model_colors <- c(CESM1 = cbPalette[1], E3SM = cbPalette[2], GISS = cbPalette[3])
+                  #,
+                  # CESM2 = cbPalette[4], MIROC = cbPalette[5], NorESM2 = cbPalette[6],
+                  # GFDL = cbPalette[7], OsloCTM3 = cbPalette[8], UKESM = cbPalette[9],
+                  # GEOS = cbPalette[10], CAM5 = cbPalette[11])
 
-model_symbols <- c(CESM1 = 15, E3SM = 15, GISS = 17, CESM2 = 19, MIROC = 15,
-                   NorESM2 = 17, GFDL = 19, OsloCTM3 = 19, UKESM = 15, GEOS = 17,
-                   CAM5 = 17)
+model_symbols <- c(CESM1 = 15, E3SM = 15, GISS = 17) #, CESM2 = 19, MIROC = 15,
+                   #NorESM2 = 17, GFDL = 19, OsloCTM3 = 19, UKESM = 15, GEOS = 17,
+                   #CAM5 = 17)
 
 #-------------------------------------------------------------------------------
 #extracts data for each perturbation experiment from csv files
@@ -126,11 +127,11 @@ data_accumulation <- function(emi_dir, reg_name, exper){
   # to be consistent with convention (i.e. positive value denotes a heating effect),
   # then take the average over all years for each variable and calculate std dev
   regional_data_summary <- regional_data %>%
-    dplyr::group_by(variable, model) %>%
     within(value <- ifelse(variable == "so2", 64.066 / 28.96, 1) * value) %>%
     within(value <- ifelse(variable %in% c("rlut", "rsut", "rlutcs", "rsutcs"), -1, 1) * value) %>%
     within(value <- ifelse(variable %in% c("wetbc", "wetso2", "wetso4") & model == "CESM2", -1, 1) * value) %>%
     within(value <- ifelse(variable == "dms", 62.13 / 28.96, 1) * value) %>%
+    dplyr::group_by(variable, model) %>%
     dplyr::summarise(regional_data = mean(value), regional_data_sd = sd(value))
 
   return(regional_data_summary)
@@ -394,6 +395,16 @@ if (sort_by == "region"){
   #creates a summary_long with the data being focused on
   summary_long <- create_summary_long(region, emi_dir)
 
+  # source("C:/Users/such559/Documents/Emissions-MIP_Data/code/cobine_so2_so4.R")
+  #
+  # summed_so2_so4_adjusted <- summed_so2_so4[[paste0('so4_so2_sum_',region)]] %>%
+  #   dplyr::rename(value = regional_data) %>%
+  #   dplyr::rename(sd = regional_data_sd)
+  # summed_so2_so4_adjusted <- summed_so2_so4_adjusted[,c('variable','model','experiment','value','sd')]
+  #
+  # summary_long <- summary_long %>%
+  #   rbind(summed_so2_so4_adjusted)
+
   if(fixed_data[1, 'fixed_by'] == 'group'){
     #filter out each variable and determine max and min
 
@@ -493,6 +504,16 @@ if (sort_by == "experiment"){
   #read in data for each region
   summary_long <- create_summary_long(exper, emi_dir)
 
+  #source("C:/Users/such559/Documents/Emissions-MIP_Data/code/cobine_so2_so4.R")
+#
+#   summed_so2_so4_adjusted <- summed_so2_so4[[paste0('so4_so2_sum_',region)]] %>%
+#     dplyr::rename(value = regional_data) %>%
+#     dplyr::rename(sd = regional_data_sd)
+#   summed_so2_so4_adjusted <- summed_so2_so4_adjusted[,c('variable','model','experiment','value','sd')]
+#
+#   summary_long <- summary_long %>%
+#     rbind(summed_so2_so4_adjusted)
+
   if(fixed_data[1, 'fixed_by'] == 'group'){
 
     #Create an empty variable dataframe list
@@ -587,6 +608,13 @@ so2_timescale_to_fix <- total_vars[[so2_timescale_index]]
 so4_lifetime_to_fix <- total_vars[[so4_lifetime_index]]
 total_vars[[so2_timescale_index]] <- mutate(so2_timescale_to_fix, value = value / 86400, sd = sd / 86400)
 total_vars[[so4_lifetime_index]] <- mutate(so4_lifetime_to_fix, value = value / 86400, sd = sd / 86400)
+#Convert the max min data from seconds to days as well
+variables['so2_timescale','Max'] <- variables['so2_timescale','Max'] / 86400
+variables['so2_timescale','Min'] <- variables['so2_timescale','Min'] / 86400
+variables['so4_lifetime','Max'] <- variables['so4_lifetime','Max'] / 86400
+variables['so4_lifetime','Min'] <- variables['so4_lifetime','Min'] / 86400
+
+test <- as.character(variables[2,3])
 
 #Creates a function that creates plots for the data based on each species
 if (sort_by == "region"){

code/summary_plots_per_diff.R

@@ -1,7 +1,7 @@
 # ------------------------------------------------------------------------------
 # Program Name: summary_plots_per-diff.R
 # Authors: Hamza Ahsan, Harrison Suchyta
-# Date Last Modified: January 14, 2022
+# Date Last Modified: October 9, 2022
 # Program Purpose: Produces summary plots of the difference between the
 # perturbations and the reference case averaged over all years
 # Input Files: ~Emissions-MIP/input/
@@ -26,7 +26,7 @@ setwd(paste0(emi_dir))
 sort_by <- 'region'
 
 #determines which region or experiment is sorted by
-region <- 'arctic'
+region <- 'global'
 exper <- 'bc-no-season'
 
 # Specify what you are sorting by and either the region (i.e., global, land, sea, arctic, NH-land, NH-sea, SH-land, SH-sea) or experiment (i.e., bc-no-season, high-so4, no-so4, reference, so2-at-height, so2-no-season)
@@ -238,7 +238,6 @@ if(sort_by == "region"){
     }
 
     # Bind data together
-    #summary_data <- list(bc_no_season_summary, high_so4_summary, no_so4_summary, so2_at_height_summary, so2_no_season_summary) %>% reduce(left_join, by = c("variable", "model"))
     summary_data <- summary_data_list %>% reduce(left_join, by = c("variable", "model"))
 
     # Correct model names for CESM and CESM2
@@ -325,7 +324,6 @@ find_max_min <- function(variable_data, variable, varname){
   return(variable_data)
 }
 #-------------------------------------------------------------------------------
-#Creates a function that filters species
 #creates a function that filters species out of a database
 filter_species <- function(database, species){
   species <- dplyr::filter(database, variable == species)
@@ -745,4 +743,4 @@ grid.newpage()
 grid.draw(deposition_plot)
 grid.newpage()
 grid.draw(column_plot)
-dev.off()
+dev.off()