gcam-korea-release: Build the GCAM-Korea Input Datasets

Documented in module_energy_LA114.wind

#' module_energy_LA114.wind
#'
#' Calculate regional supply curves for wind using country-level supply curves.
#'
#' @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{L114.RsrcCurves_EJ_R_wind}. The corresponding file in the
#' original data system was \code{LA114.wind.R} (energy level1).
#' @details Regional supply curves are generated by summing up country-level supply curves and then smoothing them out.
#' Country-level supply curves are generated using data from Zhou. Every country within a region shares the
#' the same price points for which to calculate supply. The minimum price is determined to be the lowest base price
#' within that region; maximum prices are calculated to be the price where supply is 95% that of the available resource,
#' and the highest maximum price within the region is then used. The rest of the price points are filled out in equal increments.
#' @importFrom assertthat assert_that
#' @importFrom dplyr filter mutate select
#' @importFrom tidyr gather spread
#' @author AJS June 2017
module_energy_LA114.wind <- function(command, ...) {
  if(command == driver.DECLARE_INPUTS) {
    return(c(FILE = "common/iso_GCAM_regID",
             FILE = "energy/Zhou_wind_supply_ctry_EJ"))
  } else if(command == driver.DECLARE_OUTPUTS) {
    return(c("L114.RsrcCurves_EJ_R_wind"))
  } else if(command == driver.MAKE) {

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

    # Load required inputs
    iso_GCAM_regID <- get_data(all_data, "common/iso_GCAM_regID")
    Zhou_wind_supply_ctry_EJ <- get_data(all_data, "energy/Zhou_wind_supply_ctry_EJ")

    # This routine takes country level wind resource information and aggregates that to GCAM regions.
    # The country level information is provided as GCAM smooth renewable resource parameters
    # (e.g. maxSubResource | mid-price | curve-exponent | base-price). Then:
    # 1) Each country level curve is expressed numerically by evaluating each country's curve
    # 2) The resulting numerical resource curves for each country within a GCAM region are added together
    # 3) And a new aggregate smooth renewable resource parameters are derived for that region by re-fitting
    #    the smooth renewable resource curve.

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

    # Silence package notes
    `mid-price` <- `curve-exponent` <- `base-price` <- base.price <- iso <- GCAM_region_ID <- curve.exponent <-
      max.price <- maxSubResource <- p <- supply <- optimize <- resource <- subresource <- maxSubResource <-
      optimize <- mid.price <- . <- price <- NULL

    # First, define three useful functions that will be used later: evaluate_smooth_res_curve,
    # smooth_res_curve_approx_error, and generate_max_prices.

    # The function, evaluate_smooth_res_curve, computes the smooth renewable resource function
    # supply = (p - base.price) ^ curve.exponent / (mid.price ^ curve.exponent +
    # (p - base.price) ^ curve.exponent * maxSubResource
    # Note that all of these can be vectors
    evaluate_smooth_res_curve <- function(curve.exponent, mid.price, base.price, maxSubResource, p) {
      supply <- ((p - base.price) ^ curve.exponent) / (mid.price ^ curve.exponent + ((p - base.price) ^ curve.exponent)) * maxSubResource
      # zero out the supply where the price was less than the base.price
      supply[p < base.price] <- 0
      supply
    }

    # The function, smooth_res_curve_approx_error, checks how well the given smooth renewable curve matches the given supply-points.
    # Note that the first argument is the one that is changed by optimize when trying to minimize the error
    smooth_res_curve_approx_error <- function(curve.exponent, mid.price, base.price, maxSubResource, supply_points) {
      f_p <- evaluate_smooth_res_curve(curve.exponent, mid.price, base.price, maxSubResource, supply_points$price)
      error <- f_p - supply_points$supply
      crossprod(error, error)
    }

    # The function, generate_max_prices, calculates the maximum price for each country supply curve where supply is 95% of the maximum resource.
    generate_max_prices <- function(base.price, mid.price, curve.exponent) {
      # Max price is defined to be the price at which supply is 95% of the maximum resource.
      # This range should cover the entire curve with sufficient density.
      # Note that further down the code, defining the max price to be where supply is 95% of max resource instead of 99%
      # saves a lot of data (31K rows for "regional_price_supply_points" instead of 1.6 million rows)
      L114.supply_tol <- 0.95
      (L114.supply_tol / (1-L114.supply_tol)) ^ (1 / curve.exponent) * mid.price + base.price # max price
    }

    # Country-level supply curves from the Zhou file are used to build up the regional supply curves
    Zhou_wind_supply_ctry_EJ %>%
      rename(mid.price = `mid-price`, curve.exponent = `curve-exponent`, base.price = `base-price`) %>% # Modifying title names
      left_join_error_no_match(iso_GCAM_regID, by = "iso") %>% # Attach region IDs to the countries
      # Convert prices per unit of energy from 2007 USD/kWh to 1975 USD/GJ
      mutate(mid.price = mid.price * gdp_deflator(1975, 2007) / CONV_KWH_GJ,
             base.price = base.price * gdp_deflator(1975, 2007) / CONV_KWH_GJ) ->
      L114.RsrcCurves_EJ_ctry_wind

    # Evaluate all supply curves by GCAM_region_ID to generate (price, supply) points
    # All countries within each GCAM region ID will share the same price axis. The minimum and maximum prices
    # will be from the countries within the region with the lowest minimum and highest maximum prices, respectively.
    # The minimum prices are taken from the Zhou data. The maximum prices are calculated using the function,
    # generate_max_prices, which will determine the price at which supply is 95% of the max resource.

    # The price increment between min and max price is set at 0.1 USD
    L114.price_inc <- 0.1 # 1975 USD

    # The min and max price points for each region will be determined. Remember that the minimum prices are taken
    # from the Zhou data. The maximum prices are calculated using the function, generate_max_prices.
    L114.RsrcCurves_EJ_ctry_wind %>%
      # Calculate max price using the function, generate_max_prices
      mutate(max.price = generate_max_prices(base.price, mid.price, curve.exponent)) %>%
      # Select minimum and maxium price for each GCAM region
      group_by(GCAM_region_ID) %>% # Note that when only one group is selected, no need for "ungroup" afterwards
      summarise(min.price = min(base.price), max.price = max(max.price)) %>%
      ungroup ->
      regional_min_max_price_points

    # Supply curves for each GCAM region ID will be generated.
    # First, a set of price points will be generated for each region and applied to each country.
    # Supply will then be calculated at each price point for each country.
    # Finally, supply will be summed for each GCAM region ID.
    regional_price_supply_points <- lapply(regional_min_max_price_points$GCAM_region_ID, function(region_ID) {

      # Data will be evaulated for each GCAM region ID
      data <- filter(regional_min_max_price_points, GCAM_region_ID == region_ID)

      # Generates all the price points by establishing points at 0.1 USD increments from min to max
      prices <- tibble(p = seq(data$min.price, data$max.price, L114.price_inc))

      # Applies the price points to each country (maintaining essential curve data)
      L114.RsrcCurves_EJ_ctry_wind %>%
        filter(GCAM_region_ID == region_ID) %>%
        select(curve.exponent, mid.price, base.price, maxSubResource) %>%
        repeat_add_columns(prices) ->
        ctry_curves

      # Calculate supply at each price point
      supply_points <- do.call(evaluate_smooth_res_curve, ctry_curves)

      # Sum up the supply at each price point for all countries within the region
      ctry_curves %>%
        select(p) %>%
        mutate(GCAM_region_ID = region_ID, supply = supply_points) %>%
        group_by(GCAM_region_ID, p) %>%
        summarise(supply = sum(supply)) %>%
        ungroup() %>%
        rename(price = p)

    }) # End of lapply

    # One data frame was created for each GCAM region. Need to combine them all into one.
    regional_price_supply_points <- bind_rows(regional_price_supply_points)


    # Refit the supply curve to the supply points we just calculated by evaluating the ctry level supply curves
    # so that we can aggregate them to the regional level.
    L114.RsrcCurves_EJ_R_wind <- lapply(unique(regional_price_supply_points$GCAM_region_ID), function(region_ID) {

      # Data will be evaulated for each GCAM region ID
      L114.SupplyPoints.currR <- filter(regional_price_supply_points, GCAM_region_ID == region_ID)

      # Creating table with region ID, base-price, and supply
      L114.RsrcCurves_EJ_R_wind.currR <- tibble(GCAM_region_ID = region_ID, base.price = min(L114.SupplyPoints.currR$price),
                                                maxSubResource = max(L114.SupplyPoints.currR$supply))

      # The points were calculated in increasing price order with upward sloping curves. Thus we can
      # approximate the mid.price by looking at the first price point of the supply that is 50% of the max.
      half_max_resource = 0.5 * L114.RsrcCurves_EJ_R_wind.currR$maxSubResource

      L114.SupplyPoints.currR %>%
        arrange(price) %>%
        filter(supply >= half_max_resource) %>%
        # The first price point most closely corresponds to the mid.price, which is the price at which
        # supply is half of the maximum resource. Select the first row to get the mid.price.
        filter(row_number() == 1) %>%
        pull(price) ->
        mid.price

      # Add mid-price to table
      L114.RsrcCurves_EJ_R_wind.currR %>%
        repeat_add_columns(tibble::tibble(mid.price = mid.price)) %>%
        # GCAM assumes that the base.price has already been taken out of the mid.price
        mutate(mid.price = mid.price - base.price) ->
        L114.RsrcCurves_EJ_R_wind.currR

      # We must solve for the curve exponent that best fits the supply-points
      L114.error_min_curve.exp <- optimize(f = smooth_res_curve_approx_error, interval = c(1.0, 15.0),
                                           L114.RsrcCurves_EJ_R_wind.currR$mid.price,
                                           L114.RsrcCurves_EJ_R_wind.currR$base.price,
                                           L114.RsrcCurves_EJ_R_wind.currR$maxSubResource,
                                           L114.SupplyPoints.currR)

      # Add curve exponent to table
      L114.RsrcCurves_EJ_R_wind.currR$curve.exponent <- L114.error_min_curve.exp$minimum

      L114.RsrcCurves_EJ_R_wind.currR

    }) # End of lapply


  L114.RsrcCurves_EJ_R_wind %>%
    # One data frame was created for each GCAM region. Need to combine them all into one.
    bind_rows() %>%
    # Creating "resource" and "subresource" columns that report "onshore wind resource"
    mutate(resource = "onshore wind resource",
           subresource = "onshore wind resource") %>%
    # Reorder columns
    select(GCAM_region_ID, resource, subresource, maxSubResource, mid.price, curve.exponent, base.price) %>%

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

    add_title("Wind resource curves by GCAM region") %>%
    add_units("Energy in EJ, price per unit energy in $1975/GJ") %>%
    add_comments("Regional supply curves were generated by summing up country-level supply curves and then smoothing them out.") %>%
    add_comments("Country-level supply curves were generated using data from Zhou.
                   All countries within a region were evaluated along the same price points.") %>%
    add_legacy_name("L114.RsrcCurves_EJ_R_wind") %>%
    add_precursors("common/iso_GCAM_regID", "energy/Zhou_wind_supply_ctry_EJ") ->
    L114.RsrcCurves_EJ_R_wind

  return_data(L114.RsrcCurves_EJ_R_wind)
} else {
  stop("Unknown command")
}
}