Add capability to use different SSP yield growth rate from gcamdata b…

…ased on the climate_scenario varaible
JGCRI · Oct 1, 2024 · 982cff3 · 982cff3
1 parent a71a2f9
commit 982cff3
Show file tree

Hide file tree

Showing 17 changed files with 192 additions and 71 deletions.
diff --git a/R/gcam_agprodchange.R b/R/gcam_agprodchange.R
@@ -15,10 +15,12 @@
 #' Mapping from mirca crop to GCAM commodities
 #'
 #' @param gcam_version Default = 'gcam7'. string for the GCAM version. Only support gcam6 and gcam7
+#' @param climate_scenario Default = NULL. string for climate scenario (e.g., 'ssp245')
 #' @keywords internal
 #' @export
 
-mirca_to_gcam <- function(gcam_version = 'gcam7')
+mirca_to_gcam <- function(gcam_version = 'gcam7',
+                          climate_scenario= NULL)
 {
 
   subRegionMap <- country_name <- AgSupplySector <- crop_type <- GCAM_commod <-
@@ -32,7 +34,8 @@ mirca_to_gcam <- function(gcam_version = 'gcam7')
   # ----------------------------------------------------------------------------
   # input reference ag productivity change
   # note that different GCAM versions will have different agprodchange structure
-  agprodchange_ni <- agprodchange_ref(gcam_version = gcam_version)
+  agprodchange_ni <- agprodchange_ref(gcam_version = gcam_version,
+                                      climate_scenario = climate_scenario)
 
 
   # ----------------------------------------------------------------------------
@@ -188,17 +191,20 @@ get_mirca_cropland <- function(raster_brick = NULL,
 #' cropland area within intersected region-glu
 #'
 #' @param gcam_version Default = 'gcam7'. string for the GCAM version. Only support gcam6 and gcam7
+#' @param climate_scenario Default = NULL. string for climate scenario (e.g., 'ssp245')
 #' @keywords internal
 #' @export
 
-get_cropland_weight <- function(gcam_version = 'gcam7')
+get_cropland_weight <- function(gcam_version = 'gcam7',
+                                climate_scenario = NULL)
 {
 
   region_name <- basin_name <- glu <- region_id <- glu_id <- glu_name <- crop <-
     croparea_to <- croparea_from <- country_name <- irr <- irrtype <- iso <- NULL
 
   # get mapping
-  mp <- mirca_to_gcam(gcam_version = gcam_version)
+  mp <- mirca_to_gcam(gcam_version = gcam_version,
+                      climate_scenario = climate_scenario)
 
   mp_iso <- mp$mp_iso
   mp_rmap <- mp$mp_rmap
@@ -278,6 +284,7 @@ get_cropland_weight <- function(gcam_version = 'gcam7')
 #'
 #' @param data Default = NULL. output data from function yield_shock_projection, or similar format of data
 #' @param gcam_version Default = 'gcam7'. string for the GCAM version. Only support gcam6 and gcam7
+#' @param climate_scenario Default = NULL. string for climate scenario (e.g., 'ssp245')
 #' @param diagnostics Default = TRUE. Logical for performing diagnostic plot
 #' @param output_dir Default = file.path(getwd(), 'output'). String for output directory
 #'
@@ -286,6 +293,7 @@ get_cropland_weight <- function(gcam_version = 'gcam7')
 
 get_weighted_yield_impact <- function(data = NULL,
                                       gcam_version = 'gcam7',
+                                      climate_scenario = NULL,
                                       diagnostics = TRUE,
                                       output_dir = file.path(getwd(), 'output'))
 {
@@ -296,7 +304,8 @@ get_weighted_yield_impact <- function(data = NULL,
 
   # get weight of cropland area within the intersected region-glu-country to
   # cropland area within intersected region-glu
-  out_list <- get_cropland_weight(gcam_version = gcam_version)
+  out_list <- get_cropland_weight(gcam_version = gcam_version,
+                                  climate_scenario = climate_scenario)
 
   weight <- out_list$weight
   cropland_glu_region <- out_list$cropland_glu_region
@@ -474,12 +483,14 @@ gcam_agprodchange <- function(data = NULL,
 
   yield_impact_clean <- get_weighted_yield_impact(data = data,
                                                   gcam_version = gcam_version,
+                                                  climate_scenario = climate_scenario,
                                                   diagnostics = diagnostics,
                                                   output_dir = output_dir)
 
   # input reference ag productivity change without climate impact
   # note that different GCAM versions will have different agprodchange structure
-  agprodchange_ni <- agprodchange_ref(gcam_version = gcam_version)
+  agprodchange_ni <- agprodchange_ref(gcam_version = gcam_version,
+                                      climate_scenario = climate_scenario)
 
   # linear interpolate at 5 year interval
   yield_impact <- dplyr::bind_rows(

diff --git a/R/pkg_utils.R b/R/pkg_utils.R
@@ -235,23 +235,37 @@ merge_data <- function( d1, d2, x1, x2 )
 #' Get the reference agricultural productivity change based GCAM version
 #'
 #' @param gcam_version Default = 'gcam7'. string for the GCAM version. Only support gcam6 and gcam7
+#' @param climate_scenario Default = NULL. string for climate scenario (e.g., 'ssp245')
 #' @param gcamdata_dir Default = NULL. string for directory to the gcamdata folder within the specific GCAM version. The gcamdata need to be run with drake to have the CSV outputs beforehand.
 #' @keywords internal
 #' @export
 
 agprodchange_ref <- function(gcam_version = 'gcam7',
+                             climate_scenario = NULL,
                              gcamdata_dir = NULL)
 {
 
   year <- AgProdChange <- AgProdChange_ni <- NULL
 
   if(!is.null(gcamdata_dir)){
 
+    # If user provide their own gcamdata dirctory, then use user provided data
     gaea::path_check(gcamdata_dir)
 
+    if(grepl('ssp1|ssp5', climate_scenario)){
+      agprodchange_ag <- data.table::fread(file.path(gcamdata_dir, 'outputs', 'L2052.AgProdChange_irr_high.csv'))
+    }else if(grepl('ssp3', climate_scenario)){
+      agprodchange_ag <- data.table::fread(file.path(gcamdata_dir, 'outputs', 'L2052.AgProdChange_irr_low.csv'))
+    }else if(grepl('ssp4', climate_scenario)){
+      agprodchange_ag <- data.table::fread(file.path(gcamdata_dir, 'outputs', 'L2052.AgProdChange_irr_ssp4.csv'))
+    }else{
+      agprodchange_ag <- data.table::fread(file.path(gcamdata_dir, 'outputs', 'L2052.AgProdChange_ag_irr_ref.csv'))
+    }
+
     agprodchange_ni <- dplyr::bind_rows(
       data.table::fread(file.path(gcamdata_dir, 'outputs', 'L2052.AgProdChange_bio_irr_ref.csv')),
-      data.table::fread(file.path(gcamdata_dir, 'outputs', 'L2052.AgProdChange_ag_irr_ref.csv'))) %>%
+      agprodchange_ag
+      ) %>%
       dplyr::mutate(year = paste0('X', year)) %>%
       dplyr::rename(AgProdChange_ni = AgProdChange)
 
@@ -263,8 +277,10 @@ agprodchange_ref <- function(gcam_version = 'gcam7',
         dplyr::mutate(year = 'X2015',
                       AgProdChange_ni = 0)
     )
+
   } else {
 
+    # if user doesn't provide any gcamdata, then use default
     if(gcam_version == 'gcam6'){
       agprodchange_ni <- gaea::agprodchange_ni_gcam6
     }
@@ -273,6 +289,28 @@ agprodchange_ref <- function(gcam_version = 'gcam7',
       agprodchange_ni <- gaea::agprodchange_ni_gcam7
     }
 
+    # filter based on the scenario
+    if(grepl('ssp1|ssp5', climate_scenario)){
+      # for ssp1 and ssp5
+      agprodchange_ni <- agprodchange_ni %>%
+        dplyr::select(region, AgSupplySector, AgSupplySubsector, AgProductionTechnology, year, AgProdChange_ni = high)
+
+    }else if(grepl('ssp3', climate_scenario)){
+      # for ssp3
+      agprodchange_ni <- agprodchange_ni %>%
+        dplyr::select(region, AgSupplySector, AgSupplySubsector, AgProductionTechnology, year, AgProdChange_ni = low)
+
+    }else if(grepl('ssp4', climate_scenario)){
+      # for ssp4
+      agprodchange_ni <-agprodchange_ni %>%
+        dplyr::select(region, AgSupplySector, AgSupplySubsector, AgProductionTechnology, year, AgProdChange_ni = ssp4)
+
+    }else{
+      # for both ssp2 and other non-ssp related climate scenarios (e.g., rcp45)
+      agprodchange_ni <- agprodchange_ni %>%
+        dplyr::select(region, AgSupplySector, AgSupplySubsector, AgProductionTechnology, year, AgProdChange_ni = ref)
+    }
+
   }
 
   return(agprodchange_ni)

diff --git a/R/yield_impacts_functions.R b/R/yield_impacts_functions.R
@@ -239,7 +239,7 @@ data_merge <- function(data = NULL,
 #'
 #' @param data Default = NULL. data table for crop data created by crop_month
 #' @param climate_model Default = NULL. string for climate model
-#' @param climate_scenario Default = NULL. string for climate scenario
+#' @param climate_scenario Default = NULL. string for climate scenario (e.g., 'ssp245')
 #' @param crop_name Default = NULL. string for crop name
 #' @param output_dir Default = NULL. string for path to the output folder
 #' @param co2_proj Default = NULL. data table for future CO2 concentration [year, co2_conc]
@@ -459,7 +459,7 @@ plot_fit <- function(data = NULL,
 
   d <- data
 
-  p <- ggplot2::ggplot( d, ggplot2::aes( "yield", fit_name, size = "area_harvest", color = "GCAM_region_name" ) ) +
+  p <- ggplot2::ggplot( d, ggplot2::aes_string( x = 'yield', y = fit_name, size = 'area_harvest', color = 'GCAM_region_name' ) ) +
     ggplot2::geom_point( shape = 21, stroke = 0.5 ) +
     ggplot2::scale_size_area( max_size = 20 ) +
     ggplot2::guides( color = ggplot2::guide_legend( ncol = 1 ) ) +
@@ -502,7 +502,7 @@ plot_fit <- function(data = NULL,
 #'
 #' @param use_default_coeff Default = FALSE. binary for using default regression coefficients. Set to TRUE will use the default coefficients instead of calculating coefficients from the historical climate data.
 #' @param climate_model Default = NULL. string for climate model name
-#' @param climate_scenario Default = NULL. string for RCP
+#' @param climate_scenario Default = NULL. string for climate scenario (e.g., 'ssp245')
 #' @param crop_name Default = NULL. string for crop name
 #' @param output_dir Default = file.path(getwd(), 'output'). String for output directory
 #' @keywords internal
@@ -615,7 +615,7 @@ z_estimate <- function(use_default_coeff = FALSE,
 #'
 #' @param use_default_coeff Default = FALSE. binary for using default regression coefficients. Set to TRUE will use the default coefficients instead of calculating coefficients from the historical climate data.
 #' @param climate_model Default = NULL. string for climate model name
-#' @param climate_scenario Default = NULL. string for RCP
+#' @param climate_scenario Default = NULL. string for climate scenario (e.g., 'ssp245')
 #' @param crop_name Default = NULL. string for crop name
 #' @param base_year Default = NULL. integer for the base year (for GCAM)
 #' @param start_year Default = NULL. integer for the  start year of the data
@@ -695,7 +695,7 @@ climate_impact <- function(use_default_coeff = FALSE,
 #' E.g., 2020_avg = mean(2011:2030)
 #'
 #' @param climate_model Default = NULL. string for climate model name
-#' @param climate_scenario Default = NULL. string for RCP
+#' @param climate_scenario Default = NULL. string for climate scenario (e.g., 'ssp245')
 #' @param crop_name Default = NULL. string for crop name
 #' @param base_year Default = NULL. integer for the base year (for GCAM)
 #' @param start_year Default = NULL. integer for the  start year of the data
@@ -803,7 +803,7 @@ smooth_impacts <- function(data = NULL,
 #'
 #' Format smoothed yield impact projection to wider
 #'
-#' @param data Default = NULL. data frame from yield_projections with columns [crop, model, iso, years]
+#' @param data Default = NULL. data frame from yield_shock_projection with columns [crop, model, iso, years]
 #' @param base_year Default = NULL. integer for the base year (for GCAM)
 #' @param select_years Default = NULL. vector of integers of selected years to format
 #' @param output_dir Default = file.path(getwd(), 'output'). String for output directory
@@ -867,7 +867,7 @@ format_projection <- function(data = NULL,
 #'
 #' @param data Default = NULL. data frame for the data to plot
 #' @param climate_model Default = NULL. string for climate model name
-#' @param climate_scenario Default = NULL. string for RCP
+#' @param climate_scenario Default = NULL. string for climate scenario (e.g., 'ssp245')
 #' @param crop_name Default = NULL. string for crop name
 #' @param base_year Default = NULL. integer for the base year (for GCAM)
 #' @param output_dir Default = file.path(getwd(), 'output'). String for output directory
@@ -916,7 +916,7 @@ plot_projection <- function(data = NULL,
 #'
 #' @param data Default = NULL. data frame for the data to plot
 #' @param climate_model Default = NULL. string for climate model name
-#' @param climate_scenario Default = NULL. string for RCP
+#' @param climate_scenario Default = NULL. string for climate scenario (e.g., 'ssp245')
 #' @param crop_name Default = NULL. string for crop name
 #' @param base_year Default = NULL. integer for the base year (for GCAM)
 #' @param output_dir Default = file.path(getwd(), 'output'). String for output directory
@@ -1082,7 +1082,14 @@ plot_yield_impact <- function(data = NULL,
     GCAM_commod <- year <- glu <- irrtype <- yield_multiplier <- region_name <-
     AgProductionTechnology <- NULL
 
-  print(paste0('Plotting ', commodity, crop_type))
+  save_path <- file.path(output_dir, 'figures_yield_impacts')
+  if(!dir.exists(save_path)){
+    dir.create(save_path, recursive = TRUE)
+  }
+
+
+  print(paste0('Plotting Yield Shock for ', commodity, crop_type, ' to: ',
+               file.path(save_path, paste0(commodity, crop_type, '.png'))))
 
   df_plot <- data %>%
     dplyr::filter(GCAM_commod == commodity, crop_type %in% crop_type) %>%
@@ -1116,14 +1123,10 @@ plot_yield_impact <- function(data = NULL,
                                            'RFD' = '#5773CC')) +
     ggplot2::theme_bw()
 
-  save_path <- file.path(output_dir, 'figures_yield_impacts')
-  if(!dir.exists(save_path)){
-    dir.create(save_path, recursive = TRUE)
-  }
 
   ggplot2::ggsave(p,
                   filename = file.path(save_path, paste0(commodity, crop_type, '.png')),
-                  height = 10, width = 10, dpi = 300)
+                  height = 10, width = 10, dpi = 150)
 
 }
 
@@ -1147,8 +1150,12 @@ plot_agprodchange <- function(data = NULL,
   AgProductionTechnology <- crop <- mgmt <- year <- AgProdChange <- region <-
     irrtype <- NULL
 
+  save_path <- file.path(output_dir, 'figures_agprodchange')
+  if(!dir.exists(save_path)){
+    dir.create(save_path, recursive = TRUE)
+  }
 
-  print(paste0('Plotting ', commodity))
+  print(paste0('Plotting Productivity Change for ', commodity, ' to: ', file.path(save_path, paste0(commodity, '.png'))))
 
   df_plot <- data %>%
     tidyr::separate(col = AgProductionTechnology, sep = '_',
@@ -1162,7 +1169,8 @@ plot_agprodchange <- function(data = NULL,
 
 
   if(nrow(df_plot > 0)){
-    ggplot2::ggplot(data = df_plot,
+
+    p <- ggplot2::ggplot(data = df_plot,
                     ggplot2::aes(x = year, y = AgProdChange,
                                  group = interaction(region, AgProductionTechnology, irrtype))) +
       ggplot2::geom_line(ggplot2::aes(color = irrtype), show.legend = T) +
@@ -1174,13 +1182,10 @@ plot_agprodchange <- function(data = NULL,
                                              'RFD' = '#5773CC')) +
       ggplot2::theme_bw()
 
-    save_path <- file.path(output_dir, 'figures_agprodchange')
-    if(!dir.exists(save_path)){
-      dir.create(save_path, recursive = TRUE)
-    }
 
-    ggplot2::ggsave(filename = file.path(save_path, paste0(commodity, '.png')),
-                    height = 10, width = 10, dpi = 300)
+    ggplot2::ggsave(p,
+                    filename = file.path(save_path, paste0(commodity, '.png')),
+                    height = 10, width = 10, dpi = 150)
   } else {
     print(paste0('No data for ', commodity))
   }