R/diagnostic_plots.R

In JGCRI/osiris: Process climate impacts on agricultural yields for GCAM

#' diagnostic_plots
#'
#' Function to generate diagnostic plots to compare the agprodchange and agyield
#' for the different climate projections, for the specified crops and irrigation.
#'
#' @param hot_ssp3_apc Default = NULL. Path to hot_ssp3 agprodchange csv file.
#' @param hot_ssp5_apc Default = NULL. Path to hot_ssp3 agprodchange csv file.
#' @param cold_ssp3_apc Default = NULL. Path to hot_ssp3 agprodchange csv file.
#' @param cold_ssp5_apc Default = NULL. Path to hot_ssp3 agprodchange csv file.
#' @param reference_apc Default = NULL. Path to hot_ssp3 agprodchange csv file.
#' @param agyield Default = NULL. Path to agyield gcamdata file L101.ag_Yield_kgm2_R_C_Y_GLU.csv.
#' The 2015 ag yield value is extracted from this file.
#' @param iso_GCAM_basin_mapping Default = NULL. Path to mapping file of GCAM basin id, GLU, ISO.
#' @param write_dir Default = "diagnostic_plots". Output folder for plots.
#' @keywords test
#' @return number
#' @importFrom rlang :=
#' @importFrom magrittr %>%
#' @export
#' @examples
#' \dontrun{
#' library(osiris)
#' osiris::diagnostic_plots()
#' }

diagnostic_plots <- function(hot_ssp3_apc = NULL,
                             hot_ssp5_apc = NULL,
                             cold_ssp3_apc = NULL,
                             cold_ssp5_apc = NULL,
                             reference_apc = NULL,
                             agyield = NULL,
                             iso_GCAM_basin_mapping = NULL,
                             write_dir = "diagnostic_plots") {


  #.........................
  # Initialize
  #.........................

  rlang::inform("Starting diagnostic_plots")

  # Check write dir
  if(!dir.exists(write_dir)){dir.create(write_dir)}

  # Initialize values
  NULL -> AgSupplySector -> Irrigation -> region -> year -> AgSupplySubsector ->
    value -> name -> AgSupplySector -> Irrigation -> hot_ssp3 -> hot_ssp5 ->
    cold_ssp3 -> cold_ssp5 -> reference -> AgProductionTechnology -> AgProdChange ->
    GCAM_region_ID -> ISO -> GLU_code -> GLU_name -> hot_ssp3_agyield ->
    hot_ssp5_agyield -> cold_ssp3_agyield -> cold_ssp5_agyield -> reference_agyield -> .


  #.........................
  # Custom functions
  #.........................

  # Plotting function for AgProdChange
  apc_plot <- function(crop, irr) {
    plot_title <- paste("Comparison of AgProdChange outputs using WRF vs reference for", crop, irr)
    apc_us <- dplyr::filter(proj_apc_us, AgSupplySector == crop & Irrigation == irr) %>%
      dplyr::select(-c(region, AgSupplySector, Irrigation))
    apc_subplot <- apc_us %>% tidyr::pivot_longer(-c(year, AgSupplySubsector)) %>%
      ggplot2::ggplot(ggplot2::aes(x=year, y=value, color=name))+
      ggplot2::geom_line()+
      ggplot2::facet_wrap(~AgSupplySubsector) +
      ggplot2::labs(title = plot_title, y = "AgProdChange") +
      ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5),
                     legend.title = ggplot2::element_blank(),
                     axis.title.x = ggplot2::element_blank())

    # Save plot
    ggplot2::ggsave(paste0(write_dir,"/agprodchange_",crop,"_",irr,".png"), width=15, height=10, units="in", device="png")
  }

  # Plotting function for ag yield
  agyield_plot <- function(crop, irr) {
    plot_title <- paste("Comparison of Agyield calculated for WRF vs reference for", crop, irr)
    apc_us <- dplyr::filter(proj_apc_us, AgSupplySector == crop & Irrigation == irr) %>%
      dplyr::select(-c(region, AgSupplySector, Irrigation, hot_ssp3, hot_ssp5, cold_ssp3, cold_ssp5, reference, value))
    apc_subplot <- apc_us %>% tidyr::pivot_longer(-c(year, AgSupplySubsector)) %>%
      ggplot2::ggplot(ggplot2::aes(x=year, y=value, color=name))+
      ggplot2::geom_line()+
      ggplot2::facet_wrap(~AgSupplySubsector) +
      ggplot2::labs(title = plot_title, y = expression(Agyield~(kg/m^2))) +
      ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5),
                     legend.title = ggplot2::element_blank(),
                     axis.title.x = ggplot2::element_blank()) +
      ggplot2::scale_color_hue(labels=c("hot_ssp3_agyield" = "hot_ssp3",
                                        "hot_ssp5_agyield" = "hot_ssp5",
                                        "cold_ssp3_agyield" = "cold_ssp3",
                                        "cold_ssp5_agyield" = "cold_ssp5",
                                        "reference_agyield" = "reference"))

    # Save plot
    ggplot2::ggsave(paste0(write_dir,"/agyield_",crop,"_",irr,".png"), width=15, height=10, units="in", device="png")
  }

  #...............................
  # Read in AgProdChange, agyield, and mapping files
  #...............................

  # The projection files generated by the create_AgProdChange_xml function have
  # four lines above the header which we remove
  hot_ssp3_headers = utils::read.csv(hot_ssp3_apc, skip = 4, header = F, nrows = 1, as.is = T)
  hot_ssp3_apc = utils::read.csv(hot_ssp3_apc, skip = 5, header = F)
  colnames(hot_ssp3_apc) = hot_ssp3_headers

  hot_ssp5_headers = utils::read.csv(hot_ssp5_apc, skip = 4, header = F, nrows = 1, as.is = T)
  hot_ssp5_apc = utils::read.csv(hot_ssp5_apc, skip = 5, header = F)
  colnames(hot_ssp5_apc) = hot_ssp5_headers

  cold_ssp3_headers = utils::read.csv(cold_ssp3_apc, skip = 4, header = F, nrows = 1, as.is = T)
  cold_ssp3_apc = utils::read.csv(cold_ssp3_apc, skip = 5, header = F)
  colnames(cold_ssp3_apc) = cold_ssp3_headers

  cold_ssp5_headers = utils::read.csv(cold_ssp5_apc, skip = 4, header = F, nrows = 1, as.is = T)
  cold_ssp5_apc = utils::read.csv(cold_ssp5_apc, skip = 5, header = F)
  colnames(cold_ssp5_apc) = cold_ssp5_headers

  # Read in reference AgProdChange file
  ref_apc <- utils::read.csv(reference_apc) %>%
    dplyr::select(region, AgSupplySector, AgSupplySubsector, AgProductionTechnology, year, AgProdChange)

  # Read in agyield file
  agyield_headers <- utils::read.csv(agyield, skip = 1, header = F, nrows = 1, as.is = T)
  agyield <- utils::read.csv(agyield, skip = 2, header = F)
  colnames(agyield) <- agyield_headers

  # Extract US data
  agyield <- dplyr::filter(agyield, GCAM_region_ID == 1) %>%
    dplyr::filter(year == 2015)

  # Read in mapping file
  mapping_file <- dplyr::filter(utils::read.csv(iso_GCAM_basin_mapping), ISO == "USA") %>%
    dplyr::select(GLU_code, GLU_name)

  #.........................
  # Process AgProdChange data
  #.........................

  # Extract US data and compare
  hot_ssp3_apc_us <- dplyr::filter(hot_ssp3_apc, region == "USA")
  hot_ssp5_apc_us <- dplyr::filter(hot_ssp5_apc, region == "USA")
  cold_ssp3_apc_us <- dplyr::filter(cold_ssp3_apc, region == "USA")
  cold_ssp5_apc_us <- dplyr::filter(cold_ssp5_apc, region == "USA")
  ref_apc_us <- dplyr::filter(ref_apc, region == "USA")

  proj_apc_us <- list(hot_ssp3_apc_us, hot_ssp5_apc_us, cold_ssp3_apc_us, cold_ssp5_apc_us, ref_apc_us) %>%
    Reduce(function(dtf1,dtf2) dplyr::left_join(dtf1, dtf2, by = c("region", "AgSupplySector", "AgSupplySubsector", "AgProductionTechnology", "year")), .) %>%
    dplyr::rename("hot_ssp3" = "AgProdChange.x") %>%
    dplyr::rename("hot_ssp5" = "AgProdChange.y") %>%
    dplyr::rename("cold_ssp3" = "AgProdChange.x.x") %>%
    dplyr::rename("cold_ssp5" = "AgProdChange.y.y") %>%
    dplyr::rename("reference" = "AgProdChange")

  # Take average of hi and lo irrigation technologies
  proj_apc_us$AgProductionTechnology <- gsub('.{3}$', '', proj_apc_us$AgProductionTechnology) # Removes hi and lo
  proj_apc_us <- tidyr::separate(proj_apc_us, col=AgProductionTechnology, into=c('x', 'y', 'Irrigation'), sep='_') %>%
    dplyr::select(-c('x', 'y')) %>%
    dplyr::group_by(region, AgSupplySector, AgSupplySubsector, Irrigation, year) %>%
    dplyr::summarise_at(.vars = dplyr::vars(hot_ssp3, hot_ssp5, cold_ssp3, cold_ssp5, reference),
                        .funs = c(mean = "mean")) %>%
    dplyr::rename("hot_ssp3" = "hot_ssp3_mean") %>%
    dplyr::rename("hot_ssp5" = "hot_ssp5_mean") %>%
    dplyr::rename("cold_ssp3" = "cold_ssp3_mean") %>%
    dplyr::rename("cold_ssp5" = "cold_ssp5_mean") %>%
    dplyr::rename("reference" = "reference_mean") %>%
    dplyr::ungroup()

  # Generate AgProdChange plots
  apc_plot(crop = "Corn", irr = "IRR")
  apc_plot(crop = "Corn", irr = "RFD")
  apc_plot(crop = "Wheat", irr = "IRR")
  apc_plot(crop = "Wheat", irr = "RFD")
  apc_plot(crop = "Rice", irr = "IRR")
  apc_plot(crop = "Rice", irr = "RFD")
  apc_plot(crop = "OilCrop", irr = "IRR")
  apc_plot(crop = "OilCrop", irr = "RFD")

  #.........................
  # Process Agyield data
  #.........................

  # Map agyield value to each AgSupplySubsector
  mapping_file <- merge(mapping_file, agyield, by.x = 'GLU_code', by.y = 'GLU') %>%
    tidyr::unite(AgSupplySubsector, c("GCAM_commodity", "GLU_name")) %>%
    dplyr::select(AgSupplySubsector, value)

  # Merge mapping file with apc to get agyield into the dataframe
  proj_apc_us <- merge(proj_apc_us, mapping_file, by.x = 'AgSupplySubsector', by.y = 'AgSupplySubsector')

  # Calculate the ag yield for each case using an iterative approach over 5 year timesteps
  # The equation is y_t = y_(t-1)*(1 + A_t)^5
  proj_apc_us[ ,c('hot_ssp3_agyield', 'hot_ssp5_agyield', 'cold_ssp3_agyield', 'cold_ssp5_agyield', 'reference_agyield')] <- NA
  proj_apc_us <- dplyr::group_by(proj_apc_us, region, AgSupplySector, AgSupplySubsector, Irrigation) %>%
    dplyr::do(dplyr::add_row(., year = 2015) %>%
                tidyr::fill(region, AgSupplySector, AgSupplySubsector, Irrigation, value) %>%
                dplyr::mutate(hot_ssp3_agyield = replace(hot_ssp3_agyield, year == 2015, unique(value))) %>%
                dplyr::mutate(hot_ssp5_agyield = replace(hot_ssp5_agyield, year == 2015, unique(value))) %>%
                dplyr::mutate(cold_ssp3_agyield = replace(cold_ssp3_agyield, year == 2015, unique(value))) %>%
                dplyr::mutate(cold_ssp5_agyield = replace(cold_ssp5_agyield, year == 2015, unique(value))) %>%
                dplyr::mutate(reference_agyield = replace(reference_agyield, year == 2015, unique(value))) %>%
                dplyr::arrange(year) %>%
                dplyr::mutate(hot_ssp3_agyield = purrr::accumulate(hot_ssp3[-1], .init = hot_ssp3_agyield[1], ~ .x * (1 + .y)^5)) %>%
                dplyr::mutate(hot_ssp5_agyield = purrr::accumulate(hot_ssp5[-1], .init = hot_ssp5_agyield[1], ~ .x * (1 + .y)^5)) %>%
                dplyr::mutate(cold_ssp3_agyield = purrr::accumulate(cold_ssp3[-1], .init = cold_ssp3_agyield[1], ~ .x * (1 + .y)^5)) %>%
                dplyr::mutate(cold_ssp5_agyield = purrr::accumulate(cold_ssp5[-1], .init = cold_ssp5_agyield[1], ~ .x * (1 + .y)^5)) %>%
                dplyr::mutate(reference_agyield = purrr::accumulate(reference[-1], .init = reference_agyield[1], ~ .x * (1 + .y)^5))) %>%
    dplyr::ungroup() %>%
    tidyr::drop_na()

  # Generate Agyield plots
  agyield_plot(crop = "Corn", irr = "IRR")
  agyield_plot(crop = "Corn", irr = "RFD")
  agyield_plot(crop = "Wheat", irr = "IRR")
  agyield_plot(crop = "Wheat", irr = "RFD")
  agyield_plot(crop = "Rice", irr = "IRR")
  agyield_plot(crop = "Rice", irr = "RFD")
  agyield_plot(crop = "OilCrop", irr = "IRR")
  agyield_plot(crop = "OilCrop", irr = "RFD")

  #.........................
  # Close Out
  #.........................

  rlang::inform("diagnostic_plots completed.")
}