R/zchunk_LA143.HDDCDD-korea.R

In rohmin9122/gcam-korea-release: Build the GCAM-Korea Input Datasets

#' module_gcam.korea_LA143.HDDCDD
#'
#' Calculate share of heating and cooling degree days.
#'
#' @param command API command to execute
#' @param ... other optional parameters, depending on command
#' @return Depends on \code{command}: either a vector of required inputs,
#' a vector of output names, or (if \code{command} is "MAKE") all
#' the generated outputs: \code{L143.share_korea_Pop_HDD}, \code{L143.share_korea_Pop_CDD}, \code{L143.HDDCDD_scen_korea}. The corresponding file in the
#' original data system was \code{LA143.HDDCDD.R} (gcam-korea level1).
#' @details Estimate heating and cooling degree days for gcam-korea.
#' @importFrom assertthat assert_that
#' @importFrom dplyr filter mutate select
#' @importFrom tidyr gather spread
#' @author M. Roh
module_gcam.korea_LA143.HDDCDD <- function(command, ...) {
  if(command == driver.DECLARE_INPUTS) {
    return(c(FILE = "gcam-korea/states_subregions",
             FILE = "gcam-korea/Census_pop_hist",
             FILE = "gcam-korea/GIS/HDD_His",
             FILE = "gcam-korea/GIS/HDD_GFDL_A2",
             FILE = "gcam-korea/GIS/CDD_His",
             FILE = "gcam-korea/GIS/CDD_GFDL_A2"))
  } else if(command == driver.DECLARE_OUTPUTS) {
    return(c("L143.share_korea_Pop_HDD",
             "L143.share_korea_Pop_CDD",
             "L143.HDDCDD_scen_korea"))
  } else if(command == driver.MAKE) {

    all_data <- list(...)[[1]]

    # Load required inputs
    states_subregions <- get_data(all_data, "gcam-korea/states_subregions")
    Census_pop_hist <- get_data(all_data, "gcam-korea/Census_pop_hist")
    CDD_His <- get_data(all_data, "gcam-korea/GIS/CDD_His")
    CDD_GFDL_A2 <- get_data(all_data, "gcam-korea/GIS/CDD_GFDL_A2")
    HDD_His <- get_data(all_data, "gcam-korea/GIS/HDD_His")
    HDD_GFDL_A2 <- get_data(all_data, "gcam-korea/GIS/HDD_GFDL_A2")

    # Silence package
    state <- subregion9 <- subregion13 <- year <- value <- degree_day <- population <-
      value_kor <- variable  <- GCM <- Scen <- historical_value <-
      hist_value <- base_year <- scn_final_historical_year <- NULL

    # ===================================================

    # Part 1: Create tibbles of state share of subregion population's degree days.

    # Add subregion13 (for RECS) and subregion9 (for CBECS) columns to the population by state tibble
    # and transform to long format. This tibble will be used to calculate the share of person heating
    # degree days within census division (subregion) and state.
    Census_pop_hist %>%
      left_join_error_no_match(states_subregions %>%
                                 select(state),
                               by = "state") %>%
      gather_years ->
      Census_pop_hist_subregion

    # Add subregion13 (for RECS) and subregion9 (for CBECS) columns to the historic cooling degree days
    # latter and heating degree days. Latter on this will be used to  to calculate the share of person heating
    # degree days within census division (subregion).
    #
    # Historic cooling degree days.
    CDD_His %>%
      left_join_error_no_match(states_subregions %>%
                                 select(state),
                               by = "state") %>%
      gather_years -> CDD_His_subregion

    # Historic heating degree days.
    HDD_His %>%
      left_join_error_no_match(states_subregions %>%
                                 select(state),
                               by = "state") %>%
      gather_years -> HDD_His_subregion


    # Interpolate to fill out any missing historical years cooling degree days.
    CDD_His_subregion %>%
      complete(nesting(state), year = c(year, HISTORICAL_YEARS)) %>%
      arrange(state, year) %>%
      group_by(state) %>%
      mutate(value = approx_fun(year, value, rule = 2)) %>%
      filter(year %in% HISTORICAL_YEARS) %>%
      ungroup() -> L143.CDD_state

    # Interpolate to fill out any missing historical years heating degree days.
    HDD_His_subregion %>%
      complete(nesting(state), year = c(year, HISTORICAL_YEARS)) %>%
      arrange(state, year) %>%
      group_by(state) %>%
      mutate(value = approx_fun(year, value, rule = 2)) %>%
      filter(year %in% HISTORICAL_YEARS) %>%
      ungroup() -> L143.HDD_state


    # Multiply the historic cooling degree days tibble by the census population to calculate
    # the population cooling degree days.
    L143.CDD_state %>%
      rename(degree_day = value) %>%
      left_join_error_no_match(Census_pop_hist_subregion %>% rename(population = value),
                               by = c("year", "state")) %>%
      mutate(value = degree_day * population) %>%
      select(state, year, value) ->
      L143.Pop_CDD_state

    # Multiply the historic heating degree days tibble by the census population to calculate
    # the population heating degree days.
    L143.HDD_state %>%
      rename(degree_day = value) %>%
      left_join_error_no_match(Census_pop_hist_subregion %>% rename(population = value),
                               by = c("year", "state")) %>%
      mutate(value = degree_day * population) %>%
      select(state, year, value) ->
      L143.Pop_HDD_state

    # Aggregate the total population for each census subregion. The total populatin will be used
    # to calcualte the state's share of the subregion's population heating and cooling degree days.
    #
    # Population cooling degree days in subregion 9.
    #L143.Pop_CDD_state %>%
    #  group_by(year, subregion9) %>%
    #  summarise(value_sR9 = sum(value)) %>%
    #  ungroup ->
    #  L143.Pop_CDD_sR9

    # Population cooling degree days in subregion 13.
    #L143.Pop_CDD_state %>%
    #  group_by(year, subregion13) %>%
    #  summarise(value_sR13 = sum(value)) %>%
    #  ungroup ->
    #  L143.Pop_CDD_sR13

    # Population heating degree days in subregion 9.
    #L143.Pop_HDD_state %>%
    #  group_by(year, subregion9) %>%
    #  summarise(value_sR9 = sum(value)) %>%
    #  ungroup ->
    #  L143.Pop_HDD_sR9

    # Population heating degree days in subregion 13.
    #L143.Pop_HDD_state %>%
    #  group_by(year, subregion13) %>%
    #  summarise(value_sR13 = sum(value)) %>%
    #  ungroup ->
    #  L143.Pop_HDD_sR13

    # minyoung
    L143.Pop_CDD_state %>%
      group_by(year, state) %>%
      summarise(value_kor = sum(value))%>%
      ungroup->
      L143.Pop_CDD

    L143.Pop_HDD_state %>%
      group_by(year, state) %>%
      summarise(value_kor = sum(value))%>%
      ungroup->
      L143.Pop_HDD


    # Divide state degree days by the subregions' population degree days to calculate the
    # state's share of its subregions' population degree days.
    #
    # State share of subregion 9 population cooling degree days.
    #L143.Pop_CDD_state %>%
    #  left_join_error_no_match(L143.Pop_CDD_sR9, by = c("year", "subregion9")) %>%
    #  mutate(value = value / value_sR9, variable = "CDD") %>%
    #  select(state, subregion9, variable, year, value) %>%
    #  filter(year %in% HISTORICAL_YEARS) ->
    #  L143.share_korea_Pop_CDD_sR9

    # State share of subregion 13 population cooling degree days.
    #L143.Pop_CDD_state %>%
    #  left_join_error_no_match(L143.Pop_CDD_sR13, by = c("year", "subregion13")) %>%
    #  mutate(value = value / value_sR13, variable = "CDD") %>%
    #  select(state, subregion13, variable, year, value) %>%
    #  filter(year %in% HISTORICAL_YEARS) ->
    #  L143.share_korea_Pop_CDD_sR13

    # State share of subregion 9 population heating degree days.
    #L143.Pop_HDD_state %>%
    #  left_join_error_no_match(L143.Pop_HDD_sR9, by = c("year", "subregion9")) %>%
    #  mutate(value = value / value_sR9, variable = "HDD") %>%
    #  select(state, subregion9, variable, year, value) %>%
    #  filter(year %in% HISTORICAL_YEARS) ->
    #  L143.share_korea_Pop_HDD_sR9

    # State share of subregion 13 population heating degree days.
    #L143.Pop_HDD_state %>%
    #  left_join_error_no_match(L143.Pop_HDD_sR13, by = c("year", "subregion13")) %>%
    #  mutate(value = value / value_sR13, variable = "HDD") %>%
    #  select(state, subregion13, variable, year, value)  %>%
    #  filter(year %in% HISTORICAL_YEARS) ->
    #  L143.share_korea_Pop_HDD_sR13

    #minyoung
    L143.Pop_CDD_state %>%
      left_join_error_no_match(L143.Pop_CDD, by = c("year", "state")) %>%
      mutate(value = value / value_kor, variable = "CDD") %>%
      select(state, subregion13, variable, year, value) %>%
      filter(year %in% HISTORICAL_YEARS) ->
      L143.share_korea_Pop_CDD

    L143.Pop_HDD_state %>%
      left_join_error_no_match(L143.Pop_HDD, by = c("year", "state")) %>%
      mutate(value = value / value_kor, variable = "HDD") %>%
      select(state, variable, year, value)  %>%
      filter(year %in% HISTORICAL_YEARS) ->
      L143.share_korea_Pop_HDD

    # Part 2: Process degree days in future scenarios

    # Create a list of the future scenarios data frames to process. This method preserves the
    # file name as name of the data frames in the list.
    HDDCDD_list <- list(CDD_GFDL_A2 = CDD_GFDL_A2, HDD_GFDL_A2 = HDD_GFDL_A2)

    # Add the file name to each data frame and concatenate future scenario data frames together
    # to form a single tibble.
    HDDCDD_data <- bind_rows(HDDCDD_list, .id = 'file')

    # Parse out variable, GCM, and scenario information from the the file name in to individual columns.
    HDDCDD_data %>%
      separate(col = file, into = c("variable", "GCM", "Scen"), sep = "_") ->
      HDDCDD_data

    # Interpolate heating and cooling degree days for all historical and future years
    HDDCDD_data %>%
      gather_years %>%
      complete(nesting(variable, GCM, Scen, state), year = c(HISTORICAL_YEARS, FUTURE_YEARS)) %>%
      arrange(variable, GCM, Scen, state, year) %>%
      group_by(variable, GCM, Scen, state) %>%
      mutate(value = approx_fun(year, value, rule = 2)) %>%
      ungroup ->
      L143.HDDCDD_scen_korea

    # Set aside the scenario's final historical year because this will be used to normalize the degree days
    # latter on.
    Scen_final_historical_year <- filter(L143.HDDCDD_scen_korea, year == max(HISTORICAL_YEARS))



    # Part 3: Combine historic and future scenario degree days

    # Fill in any missing historical cooling degree days through interpolation.
    CDD_His_subregion %>%
      complete(nesting(state), year = c(HISTORICAL_YEARS)) %>%
      arrange(state, year) %>%
      #group_by(state, subregion9, subregion13) %>%
      mutate(value = approx_fun(year, value, rule = 2), variable = "CDD") %>%
      ungroup ->
      CDD_His_complete

    # Fill in any missing historical heating degree days through interpolation.
    HDD_His_subregion %>%
      complete(nesting(state), year = c(HISTORICAL_YEARS)) %>%
      arrange(state, year) %>%
      #group_by(state, subregion9, subregion13) %>%
      mutate(value = approx_fun(year, value, rule = 2), variable = "HDD") %>%
      ungroup ->
      HDD_His_complete

    # Combine the interpolated cooling and heating historic degree days into a single tibble and
    # add a column called historical value. The historical_value column will be using in the next
    # step to replace degree days for model historical years in the future scenario with historical
    # observations.
    CDD_His_complete %>%
      bind_rows(HDD_His_complete) %>%
      mutate(historical_value = TRUE)->
      DD_His


    # Replace the heating and cooling degree days values in model historic years in the future scenario
    # tibble with historic data.
    L143.HDDCDD_scen_korea %>%
      # Use left_join because we expect NAs to occur in the hist_value column when in future years.
      left_join(DD_His %>%
                  select(hist_value = value, year, variable, state, historical_value),
                by = c("state", "variable", "year")) %>%
      # When applicatble replace the GCM degree days with the histroic observations.
      mutate(value = if_else(!is.na(hist_value), hist_value, value)) %>%
      select(variable, GCM, Scen, state, year, value, historical_value) ->
      L143.HDDCDD_scen_korea_w_historical_observations


    # Normalize the future degree days by the base year value. This is necessary because sometimes there are large
    # large discrepancies between the GCM model output and the observed historical degree.

    # Start by combining the the GCM, scenario, and state degree days from the  scenario final model historical
    # year, with the scenario tibble that has the historical observation replaced values.
    L143.HDDCDD_scen_korea_w_historical_observations %>%
      # Use left_join here because we do not expect a 1:1, the GCM scenario final historical year degree day
      # will be used to normalize the future degree days in the following step.
      left_join(Scen_final_historical_year %>%
                                 select(scn_final_historical_year = value, variable, GCM, Scen, state),
                by = c("variable", "GCM", "Scen", "state")) %>%
      # Use inner join (as we do not expect a 1:1 match) to add a column of the GCM degree days in the final historical
      # year of the historical observations to the tibble.
      inner_join(L143.HDDCDD_scen_korea_w_historical_observations %>%
                  filter(year == max(HISTORICAL_YEARS)) %>%
                  select(variable, GCM, Scen, state, base_year = value),
                by = c("variable", "GCM", "Scen", "state")) %>%
      # If the historical_value is NA (aka it is a future year) then multiply by the fraction of the observation
      # degree day base year over the scenario degree day base year.
      mutate(value = if_else(is.na(historical_value), value * base_year / scn_final_historical_year, value)) ->
      L143.HDDCDD_scen_korea

    # Format for output.
    L143.HDDCDD_scen_korea %>%
      select(state, Scen, GCM, variable, year, value) ->
      L143.HDDCDD_scen_korea


    # ===================================================
    # Produce outputs
    #L143.share_korea_Pop_CDD_sR9 %>%
    #  add_title("State-level share of person cooling degree days within census division (subregion9)") %>%
    #  add_units("Unitless") %>%
    #  add_comments("Interpolate historic cooling degree data to fill in missing model historical years.") %>%
    #  add_comments("Divide state cooling degree days by the census subregion 9 total population") %>%
    #  add_legacy_name("L143.share_korea_Pop_CDD_sR9") %>%
    #  add_precursors("gcam-korea/states_subregions", "gcam-korea/Census_pop_hist",
    #                 "gcam-korea/GIS/CDD_His") ->
    #  L143.share_korea_Pop_CDD_sR9

    #L143.share_korea_Pop_CDD_sR13 %>%
    #  add_title("State-level share of person cooling degree days within subregion13") %>%
    #  add_units("Unitless") %>%
    #  add_comments("Interpolate historic cooling degree data to fill in missing model historical years.") %>%
    #  add_comments("Divide state cooling degree days by the census subregion 13 total population.") %>%
    #  add_legacy_name("L143.share_korea_Pop_CDD_sR13") %>%
    #  add_precursors("gcam-korea/states_subregions", "gcam-korea/Census_pop_hist",
    #                 "gcam-korea/GIS/CDD_His") ->
    #  L143.share_korea_Pop_CDD_sR13

    #L143.share_korea_Pop_HDD_sR9 %>%
    #  add_title("State-level share of person heating degree days within census division (subregion9)") %>%
    #  add_units("Unitless") %>%
    #  add_comments("Interpolate historic heating degree data to fill in missing model historical years.") %>%
    #  add_comments("Divide state heating degree days by the census subregion 9 total population.") %>%
    #  add_legacy_name("L143.share_korea_Pop_HDD_sR9") %>%
    #  add_precursors("gcam-korea/states_subregions", "gcam-korea/Census_pop_hist",
    #                 "gcam-korea/GIS/HDD_His") ->
    #  L143.share_korea_Pop_HDD_sR9

    #L143.share_korea_Pop_HDD_sR13 %>%
    #  add_title("State-level share of person heating degree days within census division (subregion13)") %>%
    #  add_units("Unitless") %>%
    #  add_comments("Interpolate historic heating degree data to fill in missing model historical years.") %>%
    #  add_comments("Divide state heating degree days by the census subregion 13 total population.") %>%
    #  add_legacy_name("L143.share_korea_Pop_HDD_sR13") %>%
    #  add_precursors("gcam-korea/states_subregions", "gcam-korea/Census_pop_hist",
    #                 "gcam-korea/GIS/HDD_His") ->
    #  L143.share_korea_Pop_HDD_sR13

    L143.share_korea_Pop_CDD %>%
      add_title("State-level share of person cooling degree days within subregion13") %>%
      add_units("Unitless") %>%
      add_comments("Interpolate historic cooling degree data to fill in missing model historical years.") %>%
      add_comments("Divide state cooling degree days by the census subregion 13 total population.") %>%
      add_legacy_name("L143.share_korea_Pop_CDD") %>%
      add_precursors("gcam-korea/states_subregions", "gcam-korea/Census_pop_hist",
                     "gcam-korea/GIS/CDD_His") ->
      L143.share_korea_Pop_CDD

    L143.share_korea_Pop_HDD %>%
      add_title("State-level share of person heating degree days within census division (subregion9)") %>%
      add_units("Unitless") %>%
      add_comments("Interpolate historic heating degree data to fill in missing model historical years.") %>%
      add_comments("Divide state heating degree days by the census subregion 9 total population.") %>%
      add_legacy_name("L143.share_korea_Pop_HDD") %>%
      add_precursors("gcam-korea/states_subregions", "gcam-korea/Census_pop_hist",
                     "gcam-korea/GIS/HDD_His") ->
      L143.share_korea_Pop_HDD

    L143.HDDCDD_scen_korea %>%
      add_title("Heating and cooling degree days by state and scenario") %>%
      add_units("Degree F days") %>%
      add_comments("Replace GCM degree days with historical observations.") %>%
      add_comments("Normalize future GCM degree days by the fraction of observed degree days to GCM degree days in the base year.") %>%
      add_legacy_name("L143.HDDCDD_scen_korea") %>%
      add_precursors("gcam-korea/states_subregions", "gcam-korea/Census_pop_hist", "gcam-korea/GIS/CDD_His",
                     "gcam-korea/GIS/HDD_His", "gcam-korea/GIS/HDD_GFDL_A2", "gcam-korea/GIS/CDD_GFDL_A2") ->
      L143.HDDCDD_scen_korea

    return_data(L143.share_korea_Pop_CDD,
                L143.share_korea_Pop_HDD,
                L143.HDDCDD_scen_korea)
  } else {
    stop("Unknown command")
  }
}