R/gen_lookup_tax_par_win_mac_cluster.R
In manirouhirad/MODSSAT: What the Package Does (One Line, Title Case)
Defines functions
gen_lookup_tax_par_win_mac_cluster
Documented in
gen_lookup_tax_par_win_mac_cluster
#' This function is the parallelized version of the gen_lookup_subsidy function.
#' @param tax_amount is the amount of subsidy per acre-inch of groundwater extracted. Defaults to 1.
#' @param subsidy_threshold is the threshold of subsidy, i.e., the amount of groundwater extraction above which subsidy is zero. Defaults to 400.
#' @param DSSAT_files is the directory where DSSAT files are located. Defaults to "C:/Users/manirad/Dropbox/DSSAT subregions work pre-2018 annual meeting/subregion KS files/outputs_for_econ/revised".
#' @param soil_file is the file that contains base soil file. Defaults to "C:/Users/manirad/Downloads/test/Well_Soil Type_generator_07.csv".
#' @param well_capacity_file is the file that contains base well capacity. Defaults to "C:/Users/manirad/Downloads/test/Well_Capacity_ganarator.csv".
#' @param price_file is the file that includes crop prices. Defaults to "C:/Users/manirad/Dropbox/DSSAT subregions work pre-2018 annual meeting/subregion KS files/crop_prices.csv".
#' @param fixed_cost_file is the file tha includes fixed costs (per acre) costs of production. Defaults to "C:/Users/manirad/Downloads/test/fixed_cost_input.csv".
#' @param pumping_cost is the cost of pumping an acre-inch of groundwater excluding taxes. Defaults to 3.56 which is the cost of pumpin an acre-inch of groundwater in parts of Kansas.
#' @param default_well_capacity_col_name is the name of the well capacity column generated from the MODFLOW model. Defaults to 'Well_Capacity(gpm)' as this is the original column name we started with.
#' @param soil_moisture_targets is the vector of soil moisture targets that are generated by the DSSAT model. (25, 35, 45, 55, 65, 75).
#' @param IFREQ_seq is the difference between two consequtive irrigation frequencies (IFREQ) in DSSAT. Defaults to 2.
#' @param IFREQ_interpolate is the size of interpolation interval. Defaults to 0.1 so that IFREQ_seq of 2 adds 0, 0.1, 0.2, ..., 1.9.
#' @param num_clusters is the number of cores for parallelization.
#' @return returns the output table.
#' @examples
#' \dontrun{
#' gen_lookup_tax(tax_amount = 2)
#' }
#' @export
gen_lookup_tax_par_win_mac_cluster = function(tax_amount = 1,
subsidy_threshold = 1,
DSSAT_files = "./input_files/DSSAT_files",
soil_file = "./input_files/Well_Soil Type_generator_07_1000.csv",
well_capacity_file = "./input_files/Well_Capacity_ganarator_1000.csv",
price_file = "./input_files/crop_prices.csv",
fixed_cost_file = "./input_files/fixed_cost_input.csv",
pumping_cost = 3.21,
default_well_capacity_col_name = "Well_Capacity(gpm)",
soil_moisture_targets = c(25, 35, 45, 55, 65, 75),
IFREQ_seq = 2,
IFREQ_interpolate = 0.1,
num_clusters = parallel::detectCores()-2
)
{
library(data.table)
library(parallel)
tax_amount = (tax_amount - 1)/10
print(paste("this is the marginal tax rate:", tax_amount,
"per acre-inch", sep = " "))
print(paste("this is the subsidy threshold:", subsidy_threshold,
"acre-inch", sep = " "))
FN_optim2 = function(jj) {
foo_irr_3 = well_capacity_data[Well_ID_grp == jj]
foo_irr_3[ifreq > KS_DSSAT[, max(IFREQ)], ifreq := KS_DSSAT[, max(IFREQ)]]
foo_irr_3_0 = foo_irr_3[ifreq == 0]
foo_irr_3_N0 = foo_irr_3[ifreq != 0]
tryCatch({
foo_irr_3_0 = merge(foo_irr_3_0, KS_DSSAT,
by.x = c("Soil_Type", "weather_station", "CR", "ifreq"),
by.y = c("SOIL_ID", "WSTA", "CR", "IFREQ"), allow.cartesian = T)
}, error=function(e){})
foo_irr_3_N0[, diff := (KS_DSSAT[, max(IFREQ)] - ifreq) * (1/IFREQ_interpolate)]
foo_irr_3_N0[, id := 1:.N]
tryCatch({
foo_irr_3_N0 = foo_irr_3_N0[rep(1:.N, (diff+1))]
foo_irr_3_N0[, foo := 1:.N, by=c("id")]
foo_irr_3_N0[, foo := foo - 1]
foo_irr_3_N0[, ifreq := ifreq + foo * IFREQ_interpolate]
}, error=function(e){})
foo_irr_3_N0 = merge(foo_irr_3_N0, KS_DSSAT,
by.x = c("Soil_Type", "weather_station", "CR", "ifreq"),
by.y = c("SOIL_ID", "WSTA", "CR", "IFREQ"), allow.cartesian = T)
foo_irr_3_N0[, c("id", "diff", "foo") := NULL]
foo_irr_3 = rbind(foo_irr_3_N0, foo_irr_3_0)
# foo_irr_3[, ifreq_denom := ifreq/.2]
# foo_irr_3 = foo_irr_3[ifreq_denom %in% 0:200]
foo_irr_3[ifreq == KS_DSSAT[, max(IFREQ)], profit := profit/4]
rm(foo_irr_3_N0, foo_irr_3_0)
foo_irr_3[, ifreq_cap := min(ifreq), by=c("Well_capacity", "SDAT", "tot_acres", "CR", "PAW")]
foo_irr_3[, ifreq_cap := ifreq_cap + 2]
foo_irr_3 = foo_irr_3[ifreq <= ifreq_cap]
# foo_irr_3[, id := .GRP, by=c("Soil_Type", "weather_station", "CR", "Well_ID", "Well_capacity", "tot_acres", "quarter", "PAW", "SDAT")]
foo_irr_3 = unique(foo_irr_3, by=colnames(foo_irr_3))
setkey(foo_irr_3, Well_capacity, tot_acres, CR, PAW, SDAT, ifreq, quarter)
foo_irr_3[, id := .GRP, by=c("Well_capacity", "tot_acres", "ifreq")]
#----
foo_dt1_000 = foo_irr_3[quarter == 1 & tot_acres == 0, .(Well_capacity, SDAT, tot_acres, ifreq_1 = ifreq, CR_1 = CR, PAW_1 = PAW, irrigation_1 = irrigation, profit_1 = profit)]
foo_dt2_000 = foo_irr_3[quarter == 2 & tot_acres == 0, .(Well_capacity, SDAT, tot_acres, ifreq_2 = ifreq, CR_2 = CR, PAW_2 = PAW, irrigation_2 = irrigation, profit_2 = profit)]
foo_dt3_000 = foo_irr_3[quarter == 3 & tot_acres == 0, .(Well_capacity, SDAT, tot_acres, ifreq_3 = ifreq, CR_3 = CR, PAW_3 = PAW, irrigation_3 = irrigation, profit_3 = profit)]
foo_dt4_000 = foo_irr_3[quarter == 4 & tot_acres == 0, .(Well_capacity, SDAT, tot_acres, ifreq_4 = ifreq, CR_4 = CR, PAW_4 = PAW, irrigation_4 = irrigation, profit_4 = profit)]
foo_dt3_000 = merge(foo_dt3_000, foo_dt4_000, by = c("Well_capacity", "SDAT", "tot_acres"),
allow.cartesian = T)
foo_dt2_000 = merge(foo_dt2_000, foo_dt3_000, by = c("Well_capacity", "SDAT", "tot_acres"),
allow.cartesian = T)
foo_dt1_000 = merge(foo_dt1_000, foo_dt2_000, by = c("Well_capacity", "SDAT", "tot_acres"),
allow.cartesian = T)
rm(foo_dt2_000, foo_dt3_000, foo_dt4_000)
foo_dt1_000[, `:=`(irrigation_sum, irrigation_1 + irrigation_2 + irrigation_3 + irrigation_4)]
# foo_dt1_000[, `:=`(irrigation_below, ifelse(irrigation_sum < subsidy_threshold, subsidy_threshold - irrigation_sum,
# 0))]
foo_dt1_000[, `:=`(profit_sum, profit_1 + profit_2 + profit_3 + profit_4)]
# foo_dt1_000[, `:=`(subsidy_payment, irrigation_below * subsidy_amount)]
# foo_dt1_000[, `:=`(profit_sum_sub, profit_sum + irrigation_below * subsidy_amount)]
foo_dt1_000[, `:=`(profit_sum_sub, profit_sum - irrigation_sum * tax_amount)]
foo_dt1_000[, row := .GRP, by=c("CR_1", "PAW_1",
"CR_2", "PAW_2",
"CR_3", "PAW_3",
"CR_4", "PAW_4")]
foo = foo_dt1_000[, .(row, profit_sum, profit_sum_sub, irrigation_sum)]
foo[, `:=`(mean_profit_combination, mean(profit_sum)), by = c("row")]
foo[, `:=`(mean_profit_combination_sub, mean(profit_sum_sub)), by = c("row")]
foo[, `:=`(mean_irrigation_combination, mean(irrigation_sum)), by = c("row")]
foo = unique(foo, by = c("row"))
foo[, `:=`(max_p, max(mean_profit_combination_sub))]
foo = foo[mean_profit_combination_sub == max_p,.(row, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
setkey(foo, row)
setkey(foo_dt1_000, row)
foo_dt1_000 = foo_dt1_000[foo]
foo_dt1_000 = rbind(foo_dt1_000[SDAT == min(foo_dt1_000$SDAT),.(Well_capacity, tot_acres, quarter = 1, ifreq = ifreq_1, CR = CR_1, PAW = PAW_1, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_000[SDAT == min(foo_dt1_000$SDAT),.(Well_capacity, tot_acres, quarter = 2, ifreq = ifreq_2, CR = CR_2, PAW = PAW_2, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_000[SDAT == min(foo_dt1_000$SDAT),.(Well_capacity, tot_acres, quarter = 3, ifreq = ifreq_3, CR = CR_3, PAW = PAW_3, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_000[SDAT == min(foo_dt1_000$SDAT),.(Well_capacity, tot_acres, quarter = 4, ifreq = ifreq_4, CR = CR_4, PAW = PAW_4, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
)
#----------
foo_dt1_325 = foo_irr_3[quarter == 1 & tot_acres == 32.5, .(Well_capacity, SDAT, tot_acres, ifreq_1 = ifreq, CR_1 = CR, PAW_1 = PAW, irrigation_1 = irrigation, profit_1 = profit)]
foo_dt2_325 = foo_irr_3[quarter == 2 & tot_acres == 32.5, .(Well_capacity, SDAT, tot_acres, ifreq_2 = ifreq, CR_2 = CR, PAW_2 = PAW, irrigation_2 = irrigation, profit_2 = profit)]
foo_dt3_325 = foo_irr_3[quarter == 3 & tot_acres == 32.5, .(Well_capacity, SDAT, tot_acres, ifreq_3 = ifreq, CR_3 = CR, PAW_3 = PAW, irrigation_3 = irrigation, profit_3 = profit)]
foo_dt4_325 = foo_irr_3[quarter == 4 & tot_acres == 32.5, .(Well_capacity, SDAT, tot_acres, ifreq_4 = ifreq, CR_4 = CR, PAW_4 = PAW, irrigation_4 = irrigation, profit_4 = profit)]
foo_dt3_325 = merge(foo_dt3_325, foo_dt4_325, by = c("Well_capacity", "SDAT", "tot_acres"),
allow.cartesian = T)
foo_dt2_325 = merge(foo_dt2_325, foo_dt3_325, by = c("Well_capacity", "SDAT", "tot_acres"),
allow.cartesian = T)
foo_dt1_325 = merge(foo_dt1_325, foo_dt2_325, by = c("Well_capacity", "SDAT", "tot_acres"),
allow.cartesian = T)
rm(foo_dt2_325, foo_dt3_325, foo_dt4_325)
foo_dt1_325[, `:=`(irrigation_sum, irrigation_1 + irrigation_2 + irrigation_3 + irrigation_4)]
# foo_dt1_325[, `:=`(irrigation_below, ifelse(irrigation_sum < subsidy_threshold, subsidy_threshold - irrigation_sum,
# 0))]
foo_dt1_325[, `:=`(profit_sum, profit_1 + profit_2 + profit_3 + profit_4)]
# foo_dt1_325[, `:=`(subsidy_payment, irrigation_below * subsidy_amount)]
# foo_dt1_325[, `:=`(profit_sum_sub, profit_sum + irrigation_below * subsidy_amount)]
foo_dt1_325[, `:=`(profit_sum_sub, profit_sum - irrigation_sum * tax_amount)]
foo_dt1_325[, row := .GRP, by=c("CR_1", "PAW_1",
"CR_2", "PAW_2",
"CR_3", "PAW_3",
"CR_4", "PAW_4", "ifreq_1")]
foo = foo_dt1_325[, .(row, profit_sum, profit_sum_sub, irrigation_sum)]
foo[, `:=`(mean_profit_combination, mean(profit_sum)), by = c("row")]
foo[, `:=`(mean_profit_combination_sub, mean(profit_sum_sub)), by = c("row")]
foo[, `:=`(mean_irrigation_combination, mean(irrigation_sum)), by = c("row")]
foo = unique(foo, by = c("row"))
foo[, `:=`(max_p, max(mean_profit_combination_sub))]
foo = foo[mean_profit_combination_sub == max_p,.(row, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
setkey(foo, row)
setkey(foo_dt1_325, row)
foo_dt1_325 = foo_dt1_325[foo]
foo_dt1_325 = rbind(foo_dt1_325[SDAT == min(foo_dt1_325$SDAT),.(Well_capacity, tot_acres, quarter = 1, ifreq = ifreq_1, CR = CR_1, PAW = PAW_1, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_325[SDAT == min(foo_dt1_325$SDAT),.(Well_capacity, tot_acres, quarter = 2, ifreq = ifreq_2, CR = CR_2, PAW = PAW_2, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_325[SDAT == min(foo_dt1_325$SDAT),.(Well_capacity, tot_acres, quarter = 3, ifreq = ifreq_3, CR = CR_3, PAW = PAW_3, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_325[SDAT == min(foo_dt1_325$SDAT),.(Well_capacity, tot_acres, quarter = 4, ifreq = ifreq_4, CR = CR_4, PAW = PAW_4, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
)
#----------
foo_dt1_650 = foo_irr_3[quarter == 1 & tot_acres == 65, .(Well_capacity, SDAT, tot_acres, ifreq_1 = ifreq, CR_1 = CR, PAW_1 = PAW, irrigation_1 = irrigation, profit_1 = profit, id)]
foo_dt2_650 = foo_irr_3[quarter == 2 & tot_acres == 65, .(Well_capacity, SDAT, tot_acres, ifreq_2 = ifreq, CR_2 = CR, PAW_2 = PAW, irrigation_2 = irrigation, profit_2 = profit, id)]
foo_dt3_650 = foo_irr_3[quarter == 3 & tot_acres == 65, .(Well_capacity, SDAT, tot_acres, ifreq_3 = ifreq, CR_3 = CR, PAW_3 = PAW, irrigation_3 = irrigation, profit_3 = profit)]
foo_dt4_650 = foo_irr_3[quarter == 4 & tot_acres == 65, .(Well_capacity, SDAT, tot_acres, ifreq_4 = ifreq, CR_4 = CR, PAW_4 = PAW, irrigation_4 = irrigation, profit_4 = profit)]
foo_dt3_650 = merge(foo_dt3_650, foo_dt4_650, by = c("Well_capacity", "SDAT", "tot_acres"),
allow.cartesian = T)
foo_dt2_650 = merge(foo_dt2_650, foo_dt3_650, by = c("Well_capacity", "SDAT", "tot_acres"),
allow.cartesian = T)
foo_dt1_650 = merge(foo_dt1_650, foo_dt2_650, by = c("Well_capacity", "SDAT", "tot_acres", "id"),
allow.cartesian = T)
rm(foo_dt2_650, foo_dt3_650, foo_dt4_650)
foo_dt1_650[, `:=`(irrigation_sum, irrigation_1 + irrigation_2 + irrigation_3 + irrigation_4)]
# foo_dt1_650[, `:=`(irrigation_below, ifelse(irrigation_sum < subsidy_threshold, subsidy_threshold - irrigation_sum,
# 0))]
foo_dt1_650[, `:=`(profit_sum, profit_1 + profit_2 + profit_3 + profit_4)]
# foo_dt1_650[, `:=`(subsidy_payment, irrigation_below * subsidy_amount)]
# foo_dt1_650[, `:=`(profit_sum_sub, profit_sum + irrigation_below * subsidy_amount)]
foo_dt1_650[, `:=`(profit_sum_sub, profit_sum - irrigation_sum * tax_amount)]
foo_dt1_650[, row := .GRP, by=c("CR_1", "PAW_1",
"CR_2", "PAW_2",
"CR_3", "PAW_3",
"CR_4", "PAW_4", "id")]
foo = foo_dt1_650[, .(row, profit_sum, profit_sum_sub, irrigation_sum)]
foo[, `:=`(mean_profit_combination, mean(profit_sum)), by = c("row")]
foo[, `:=`(mean_profit_combination_sub, mean(profit_sum_sub)), by = c("row")]
foo[, `:=`(mean_irrigation_combination, mean(irrigation_sum)), by = c("row")]
foo = unique(foo, by = c("row"))
foo[, `:=`(max_p, max(mean_profit_combination_sub))]
foo = foo[mean_profit_combination_sub == max_p,.(row, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
setkey(foo, row)
setkey(foo_dt1_650, row)
foo_dt1_650 = foo_dt1_650[foo]
foo_dt1_650 = rbind(foo_dt1_650[SDAT == min(foo_dt1_650$SDAT),.(Well_capacity, tot_acres, quarter = 1, ifreq = ifreq_1, CR = CR_1, PAW = PAW_1, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_650[SDAT == min(foo_dt1_650$SDAT),.(Well_capacity, tot_acres, quarter = 2, ifreq = ifreq_2, CR = CR_2, PAW = PAW_2, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_650[SDAT == min(foo_dt1_650$SDAT),.(Well_capacity, tot_acres, quarter = 3, ifreq = ifreq_3, CR = CR_3, PAW = PAW_3, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_650[SDAT == min(foo_dt1_650$SDAT),.(Well_capacity, tot_acres, quarter = 4, ifreq = ifreq_4, CR = CR_4, PAW = PAW_4, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
)
#----------
foo_dt1_975 = foo_irr_3[quarter == 1 & tot_acres == 97.5, .(Well_capacity, SDAT, tot_acres, ifreq_1 = ifreq, CR_1 = CR, PAW_1 = PAW, irrigation_1 = irrigation, profit_1 = profit, id)]
foo_dt2_975 = foo_irr_3[quarter == 2 & tot_acres == 97.5, .(Well_capacity, SDAT, tot_acres, ifreq_2 = ifreq, CR_2 = CR, PAW_2 = PAW, irrigation_2 = irrigation, profit_2 = profit, id)]
foo_dt3_975 = foo_irr_3[quarter == 3 & tot_acres == 97.5, .(Well_capacity, SDAT, tot_acres, ifreq_3 = ifreq, CR_3 = CR, PAW_3 = PAW, irrigation_3 = irrigation, profit_3 = profit, id)]
foo_dt4_975 = foo_irr_3[quarter == 4 & tot_acres == 97.5, .(Well_capacity, SDAT, tot_acres, ifreq_4 = ifreq, CR_4 = CR, PAW_4 = PAW, irrigation_4 = irrigation, profit_4 = profit)]
foo_dt3_975 = merge(foo_dt3_975, foo_dt4_975, by = c("Well_capacity", "SDAT", "tot_acres"),
allow.cartesian = T)
foo_dt2_975 = merge(foo_dt2_975, foo_dt3_975, by = c("Well_capacity", "SDAT", "tot_acres", "id"),
allow.cartesian = T)
foo_dt1_975 = merge(foo_dt1_975, foo_dt2_975, by = c("Well_capacity", "SDAT", "tot_acres", "id"),
allow.cartesian = T)
rm(foo_dt2_975, foo_dt3_975, foo_dt4_975)
foo_dt1_975[, `:=`(irrigation_sum, irrigation_1 + irrigation_2 + irrigation_3 + irrigation_4)]
# foo_dt1_975[, `:=`(irrigation_below, ifelse(irrigation_sum < subsidy_threshold, subsidy_threshold - irrigation_sum,
# 0))]
foo_dt1_975[, `:=`(profit_sum, profit_1 + profit_2 + profit_3 + profit_4)]
# foo_dt1_975[, `:=`(subsidy_payment, irrigation_below * subsidy_amount)]
# foo_dt1_975[, `:=`(profit_sum_sub, profit_sum + irrigation_below * subsidy_amount)]
foo_dt1_975[, `:=`(profit_sum_sub, profit_sum - irrigation_sum * tax_amount)]
foo_dt1_975[, row := .GRP, by=c("CR_1", "PAW_1",
"CR_2", "PAW_2",
"CR_3", "PAW_3",
"CR_4", "PAW_4", "id")]
foo = foo_dt1_975[, .(row, profit_sum, profit_sum_sub, irrigation_sum)]
foo[, `:=`(mean_profit_combination, mean(profit_sum)), by = c("row")]
foo[, `:=`(mean_profit_combination_sub, mean(profit_sum_sub)), by = c("row")]
foo[, `:=`(mean_irrigation_combination, mean(irrigation_sum)), by = c("row")]
foo = unique(foo, by = c("row"))
foo[, `:=`(max_p, max(mean_profit_combination_sub))]
foo = foo[mean_profit_combination_sub == max_p,.(row, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
setkey(foo, row)
setkey(foo_dt1_975, row)
foo_dt1_975 = foo_dt1_975[foo]
foo_dt1_975 = rbind(foo_dt1_975[SDAT == min(foo_dt1_975$SDAT),.(Well_capacity, tot_acres, quarter = 1, ifreq = ifreq_1, CR = CR_1, PAW = PAW_1, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_975[SDAT == min(foo_dt1_975$SDAT),.(Well_capacity, tot_acres, quarter = 2, ifreq = ifreq_2, CR = CR_2, PAW = PAW_2, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_975[SDAT == min(foo_dt1_975$SDAT),.(Well_capacity, tot_acres, quarter = 3, ifreq = ifreq_3, CR = CR_3, PAW = PAW_3, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_975[SDAT == min(foo_dt1_975$SDAT),.(Well_capacity, tot_acres, quarter = 4, ifreq = ifreq_4, CR = CR_4, PAW = PAW_4, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
)
#----------
foo_dt1_130 = foo_irr_3[quarter == 1 & tot_acres == 130, .(Well_capacity, SDAT, tot_acres, ifreq_1 = ifreq, CR_1 = CR, PAW_1 = PAW, irrigation_1 = irrigation, profit_1 = profit, id)]
foo_dt2_130 = foo_irr_3[quarter == 2 & tot_acres == 130, .(Well_capacity, SDAT, tot_acres, ifreq_2 = ifreq, CR_2 = CR, PAW_2 = PAW, irrigation_2 = irrigation, profit_2 = profit, id)]
foo_dt3_130 = foo_irr_3[quarter == 3 & tot_acres == 130, .(Well_capacity, SDAT, tot_acres, ifreq_3 = ifreq, CR_3 = CR, PAW_3 = PAW, irrigation_3 = irrigation, profit_3 = profit, id)]
foo_dt4_130 = foo_irr_3[quarter == 4 & tot_acres == 130, .(Well_capacity, SDAT, tot_acres, ifreq_4 = ifreq, CR_4 = CR, PAW_4 = PAW, irrigation_4 = irrigation, profit_4 = profit, id)]
foo_dt3_130 = merge(foo_dt3_130, foo_dt4_130, by = c("Well_capacity", "SDAT", "tot_acres", "id"),
allow.cartesian = T)
foo_dt2_130 = merge(foo_dt2_130, foo_dt3_130, by = c("Well_capacity", "SDAT", "tot_acres", "id"),
allow.cartesian = T)
foo_dt1_130 = merge(foo_dt1_130, foo_dt2_130, by = c("Well_capacity", "SDAT", "tot_acres", "id"),
allow.cartesian = T)
rm(foo_dt2_130, foo_dt3_130, foo_dt4_130)
foo_dt1_130[, `:=`(irrigation_sum, irrigation_1 + irrigation_2 + irrigation_3 + irrigation_4)]
# foo_dt1_130[, `:=`(irrigation_below, ifelse(irrigation_sum < subsidy_threshold, subsidy_threshold - irrigation_sum,
# 0))]
foo_dt1_130[, `:=`(profit_sum, profit_1 + profit_2 + profit_3 + profit_4)]
# foo_dt1_130[, `:=`(subsidy_payment, irrigation_below * subsidy_amount)]
# foo_dt1_130[, `:=`(profit_sum_sub, profit_sum + irrigation_below * subsidy_amount)]
foo_dt1_130[, `:=`(profit_sum_sub, profit_sum - irrigation_sum * tax_amount)]
foo_dt1_130[, row := .GRP, by=c("CR_1", "PAW_1",
"CR_2", "PAW_2",
"CR_3", "PAW_3",
"CR_4", "PAW_4", "id")]
foo = foo_dt1_130[, .(row, profit_sum, profit_sum_sub, irrigation_sum)]
foo[, `:=`(mean_profit_combination, mean(profit_sum)), by = c("row")]
foo[, `:=`(mean_profit_combination_sub, mean(profit_sum_sub)), by = c("row")]
foo[, `:=`(mean_irrigation_combination, mean(irrigation_sum)), by = c("row")]
foo = unique(foo, by = c("row"))
foo[, `:=`(max_p, max(mean_profit_combination_sub))]
foo = foo[mean_profit_combination_sub == max_p,.(row, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
setkey(foo, row)
setkey(foo_dt1_130, row)
foo_dt1_130 = foo_dt1_130[foo]
foo_dt1_130 = rbind(foo_dt1_130[SDAT == min(foo_dt1_130$SDAT),.(Well_capacity, tot_acres, quarter = 1, ifreq = ifreq_1, CR = CR_1, PAW = PAW_1, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_130[SDAT == min(foo_dt1_130$SDAT),.(Well_capacity, tot_acres, quarter = 2, ifreq = ifreq_2, CR = CR_2, PAW = PAW_2, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_130[SDAT == min(foo_dt1_130$SDAT),.(Well_capacity, tot_acres, quarter = 3, ifreq = ifreq_3, CR = CR_3, PAW = PAW_3, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_130[SDAT == min(foo_dt1_130$SDAT),.(Well_capacity, tot_acres, quarter = 4, ifreq = ifreq_4, CR = CR_4, PAW = PAW_4, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
)
#----------
foo_dt1 = rbind(foo_dt1_000, foo_dt1_325, foo_dt1_650, foo_dt1_975, foo_dt1_130)
rm(foo_dt1_000, foo_dt1_325, foo_dt1_650, foo_dt1_975, foo_dt1_130)
foo_dt1[, max_mean_p := max(mean_profit_combination_sub, na.rm = T)]
foo_dt1 = foo_dt1[max_mean_p == mean_profit_combination_sub]
return(foo_dt1)
}
col_new = c("RUNNO", "TRNO", "R_pound", "O_pound", "C_pound",
"CR", "MODEL", "EXNAME", "FNAM", "WSTA", "SOIL_ID", "SDAT",
"PDAT", "EDAT", "ADAT", "MDAT", "HDAT", "DWAP", "CWAM",
"HWAM", "HWAH", "BWAH", "PWAM", "HWUM", "H_AM", "H_UM",
"HIAM", "LAIX", "IR_M", "IRCM", "PRCM", "ETCM", "EPCM",
"ESCM", "ROCM", "DRCM", "SWXM", "NI_M", "NICM", "NFXM",
"NUCM", "NLCM", "NIAM", "CNAM", "GNAM", "N2OEC", "PI_M",
"PICM", "PUPC", "SPAM", "KI_M", "KICM", "KUPC", "SKAM",
"RECM", "ONTAM", "ONAM", "OPTAM", "OPAM", "OCTAM", "OCAM",
"CO2EC", "DMPPM", "DMPEM", "DMPTM", "DMPIM", "YPPM",
"YPEM", "YPTM", "YPIM", "DPNAM", "DPNUM", "YPNAM", "YPNUM",
"NDCH", "TMAXA", "TMINA", "SRADA", "DAYLA", "CO2A", "PRCP",
"ETCP", "ESCP", "EPCP", "PAW", "IFREQ")
filenames = list.files(path = DSSAT_files, pattern = "*.OSU",
full.names = TRUE, recursive = TRUE)
ldf <- lapply(filenames, read.table, fill = T)
ldf <- lapply(ldf, function(x) x[-2, ])
ldf <- mapply(cbind, ldf, file_name = filenames, SIMPLIFY = F)
KS_DSSAT = data.table::rbindlist(ldf)
KS_DSSAT = data.table(KS_DSSAT)
KS_DSSAT = KS_DSSAT[complete.cases(V85)]
KS_DSSAT[, `:=`(file_name, NULL)]
rm(ldf)
KS_DSSAT[, `:=`(foo, nchar(as.character(V9)))]
KS_DSSAT[foo < 4, `:=`(V9, paste(V9, V11, V10, sep = "_"))]
KS_DSSAT[foo < 4, `:=`(V10, NA)]
KS_DSSAT[foo < 4, `:=`(V11, NA)]
KS_DSSAT = KS_DSSAT[, !sapply(KS_DSSAT, function(x) all(is.na(x))),
with = FALSE]
KS_DSSAT[, `:=`(V9, as.character(V9))]
KS_DSSAT[, `:=`(PAW, substr(V9, 1, regexpr("_", V9) - 2))]
KS_DSSAT[, `:=`(IFREQ, substr(V9, regexpr("Q", V9) + 2, nchar(V9)))]
KS_DSSAT[, `:=`(V9, NULL)]
KS_DSSAT[, `:=`(foo, NULL)]
data.table::setnames(KS_DSSAT, old = colnames(KS_DSSAT),
new = col_new)
KS_DSSAT[, `:=`(WSTA, substr(WSTA, 1, 4))]
KS_DSSAT[, `:=`(group, .GRP), by = c("WSTA", "SOIL_ID")]
KS_DSSAT_2 = data.table::copy(KS_DSSAT)
lookup_table_all_years_2 = data.table::data.table()
lookup_table_quarter_2 = data.table::data.table()
lookup_table_well_2 = data.table::data.table()
for (i in 1:max(KS_DSSAT_2$group)) {
soil_type = data.table::fread(soil_file)
soil_type[, `:=`(Soil_Type, KS_DSSAT_2[group == i, unique(SOIL_ID)])]
soil_type[, `:=`(Soil_Type, gsub("KSFC00000", "KS0000000",
Soil_Type))]
well_capacity = data.table::fread(well_capacity_file)
well_capacity = well_capacity[`Well_Capacity(gpm)` <= 801]
data.table::setkey(soil_type, Well_ID)
data.table::setkey(well_capacity, Well_ID)
well_capacity_data = soil_type[well_capacity]
data.table::setnames(well_capacity_data, old = default_well_capacity_col_name,
"Well_capacity")
WSTAT = data.table(Well_ID = well_capacity_data$Well_ID,
weather_station = KS_DSSAT_2[group == i, unique(WSTA)])
data.table::setkey(WSTAT, Well_ID)
data.table::setkey(well_capacity_data, Well_ID)
well_capacity_data = WSTAT[well_capacity_data]
well_capacity_data = well_capacity_data[, .(Well_ID,
Soil_Type, weather_station, Well_capacity)]
price_dt = data.table::fread(price_file)
data.table::setkey(price_dt, CR)
cost_dt = well_capacity_data[, .(Well_ID)]
data.table::setkey(cost_dt, Well_ID)
cost_dt[, `:=`(cost_per_acre_in, (pumping_cost)/1)]
fixed_cost = data.table::fread(fixed_cost_file)
fixed_cost[irr == 0, `:=`(Crop, paste("dry", Crop, sep = "-"))]
fixed_cost[, `:=`(irr, NULL)]
fixed_cost = rbind(fixed_cost, data.table::data.table(Crop = "FA",
f_cost = 0))
data.table::setkey(fixed_cost, Crop)
KS_DSSAT = KS_DSSAT_2[group == i, .(SOIL_ID, WSTA, CR,
IFREQ, PAW, SDAT, IRCM, PRCP, PRCM, HWAM)]
KS_DSSAT[, `:=`(IRCM, as.numeric(as.character(IRCM)))]
KS_DSSAT = KS_DSSAT[complete.cases(IRCM)]
KS_DSSAT[, `:=`(SDAT, substr(SDAT, 1, 4))]
cols_change = c("IFREQ", "PAW", "SDAT", "IRCM", "PRCP",
"PRCM", "HWAM")
KS_DSSAT[, `:=`((cols_change), lapply(.SD, function(x) as.numeric(as.character(x)))),
.SDcols = cols_change]
cols_change = c("SOIL_ID", "WSTA", "CR")
KS_DSSAT[, `:=`((cols_change), lapply(.SD, as.character)),
.SDcols = cols_change]
KS_DSSAT = KS_DSSAT[IFREQ < 17]
qux1 = KS_DSSAT[IFREQ == 6 & PAW == 75]
qux2 = KS_DSSAT[IFREQ == 8 & PAW == 75]
qux1[, `:=`(IFREQ, 8)]
qux2[, `:=`(IFREQ, 6)]
KS_DSSAT = KS_DSSAT[!((IFREQ == 6 | IFREQ == 8) & PAW ==
75)]
KS_DSSAT = rbind(KS_DSSAT, qux1)
KS_DSSAT = rbind(KS_DSSAT, qux2)
KS_DSSAT = KS_DSSAT[PAW %in% soil_moisture_targets |
IFREQ == 0]
KS_DSSAT_0 = KS_DSSAT[IFREQ == 0]
KS_DSSAT = KS_DSSAT[IFREQ != 0]
KS_DSSAT_0 <- KS_DSSAT_0[rep(seq_len(nrow(KS_DSSAT_0)),
each = length(soil_moisture_targets))]
KS_DSSAT_0[, `:=`(PAW, rep(soil_moisture_targets, nrow(KS_DSSAT_0)/length(soil_moisture_targets)))]
KS_DSSAT = rbind(KS_DSSAT, KS_DSSAT_0)
data.table::setkey(KS_DSSAT, SOIL_ID, WSTA, CR, IFREQ,
PAW, SDAT)
KS_DSSAT = KS_DSSAT[PRCP > -10 & HWAM > -10]
KS_DSSAT = KS_DSSAT[, .(SOIL_ID, WSTA, CR, PAW, SDAT,
IFREQ, IRCM, PRCP, PRCM, HWAM)]
data.table::setkey(KS_DSSAT, SOIL_ID, WSTA, CR, PAW,
SDAT, IFREQ)
# foo = copy(KS_DSSAT)
# foo[, mean_HWAM := mean(HWAM), by=c("SOIL_ID", "WSTA", "CR", "IFREQ", "PAW")]
# foo[, mean_IRCM := mean(IRCM), by=c("SOIL_ID", "WSTA", "CR", "IFREQ", "PAW")]
# foo = unique(foo, by=c("SOIL_ID", "WSTA", "CR", "IFREQ", "PAW"))
#
# ggplot(foo[IFREQ == 4], aes(x=mean_IRCM, y= mean_HWAM, group=CR, col=CR))+
# geom_point()+
# geom_line()
# ggplot(KS_DSSAT[IFREQ == 4], aes(x=IRCM, y= HWAM, group=CR, col=CR))+
# geom_point()+
# geom_smooth()
#
# ggplot(KS_DSSAT[CR == "MZ"], aes(x=IRCM, y= HWAM, group=factor(IFREQ), col=factor(IFREQ)))+
# geom_point()
KS_DSSAT[, `:=`(lead_yield, dplyr::lead(HWAM, n = 1L)),
by = c("SOIL_ID", "WSTA", "CR", "PAW", "SDAT")]
KS_DSSAT[, `:=`(lead_irr_mm, dplyr::lead(IRCM, n = 1L)),
by = c("SOIL_ID", "WSTA", "CR", "PAW", "SDAT")]
KS_DSSAT <- KS_DSSAT[rep(seq_len(nrow(KS_DSSAT)), each = IFREQ_seq *
10)]
df_foo = data.table::data.table(x = seq(0, IFREQ_seq -
0.1, IFREQ_interpolate))
df_foo = do.call(rbind, replicate(nrow(KS_DSSAT)/nrow(df_foo),
df_foo, simplify = F))
KS_DSSAT = data.table::data.table(KS_DSSAT, IFREQ_int = df_foo$x)
KS_DSSAT[, `:=`(foo, max(IFREQ)), by = c("SOIL_ID", "WSTA",
"CR")]
KS_DSSAT[, `:=`(IFREQ, IFREQ + IFREQ_int)]
KS_DSSAT = KS_DSSAT[IFREQ <= foo]
KS_DSSAT[, `:=`(foo, NULL)]
KS_DSSAT = KS_DSSAT[complete.cases(lead_yield)]
KS_DSSAT[, `:=`(yield_int, HWAM + (lead_yield - HWAM)/IFREQ_seq *
IFREQ_int)]
KS_DSSAT[, `:=`(irr_int, IRCM + (lead_irr_mm - IRCM)/IFREQ_seq *
IFREQ_int)]
KS_DSSAT = KS_DSSAT[, .(SOIL_ID, WSTA, CR, IFREQ, PAW,
SDAT, irr_mm = irr_int, PRCP, PRCM, yield_kg_ac = yield_int)]
KS_DSSAT[, `:=`(yield_kg_ac, yield_kg_ac * 0.4046)]
data.table::setkey(KS_DSSAT, SOIL_ID, WSTA, CR, IFREQ,
PAW, SDAT)
KS_DSSAT = KS_DSSAT[IFREQ == 0 | IFREQ >= IFREQ_seq]
KS_DSSAT = KS_DSSAT[IFREQ != 0 | PAW == soil_moisture_targets[1]]
KS_DSSAT = KS_DSSAT[SDAT < 2016]
KS_DSSAT[IFREQ == 0, `:=`(CR, paste("dry", CR, sep = "-"))]
KS_DSSAT[CR %like% "MZ", `:=`(yield_kg_ac, yield_kg_ac *
1.183)]
KS_DSSAT[CR %like% "SG", `:=`(yield_kg_ac, yield_kg_ac *
1.183)]
KS_DSSAT[CR %like% "WH", `:=`(yield_kg_ac, yield_kg_ac *
1.156)]
number_of_crops = length(KS_DSSAT[, unique(CR)])
well_capacity_data = rbind(data.table::data.table(Well_ID = 1,
Soil_Type = unique(well_capacity_data$Soil_Type),
weather_station = unique(well_capacity_data$weather_station),
Well_capacity = 0), well_capacity_data)
well_capacity_data[, `:=`(quarter_1, 0)]
well_capacity_data[, `:=`(quarter_2, 0)]
well_capacity_data[, `:=`(quarter_3, 0)]
well_capacity_data[, `:=`(quarter_4, 0)]
well_capacity_data <- well_capacity_data[rep(seq_len(nrow(well_capacity_data)),
each = 5)]
well_capacity_data[, `:=`(tot_acres, rep(c(0, 32.5, 65,
97.5, 130), nrow(well_capacity_data)/5))]
well_capacity_data <- well_capacity_data[rep(seq_len(nrow(well_capacity_data)),
each = (number_of_crops + 1))]
well_capacity_data[, `:=`(quarter_4, rep(0:number_of_crops,
nrow(well_capacity_data)/(number_of_crops + 1)))]
well_capacity_data[, `:=`(quarter_3, rep(0:number_of_crops,
nrow(well_capacity_data)/(number_of_crops + 1)))]
well_capacity_data[, `:=`(quarter_2, rep(0:number_of_crops,
nrow(well_capacity_data)/(number_of_crops + 1)))]
well_capacity_data[, `:=`(quarter_1, rep(0:number_of_crops,
nrow(well_capacity_data)/(number_of_crops + 1)))]
well_capacity_data[, `:=`(ifreq, round((tot_acres * 1)/(Well_capacity *
0.053030303149011), 1))]
# well_capacity_data[!(complete.cases(ifreq) & ifreq < 100), `:=`(ifreq, 0)]
well_capacity_data[!complete.cases(ifreq), `:=`(ifreq, 0)]
well_capacity_data[ifreq < 2 & ifreq > 0, `:=`(ifreq,
2)]
well_capacity_data = reshape2::melt(well_capacity_data,
id = c("Well_ID", "Soil_Type", "weather_station",
"Well_capacity", "tot_acres", "ifreq"))
well_capacity_data = data.table::data.table(well_capacity_data)
data.table::setkey(well_capacity_data, Well_ID, Soil_Type,
weather_station, tot_acres)
quz = data.table(CR = KS_DSSAT[, unique(CR)])
foo = data.table(value = 0:number_of_crops, rbind(data.table(CR = "FA"),
quz))
setkey(well_capacity_data, value)
setkey(foo, value)
well_capacity_data = well_capacity_data[foo]
well_capacity_data[, `:=`(value, NULL)]
well_capacity_data[, `:=`(quarter, ifelse(variable ==
"quarter_1", 1, ifelse(variable == "quarter_2", 2,
ifelse(variable == "quarter_3", 3, 4))))]
well_capacity_data[, `:=`(variable, NULL)]
tot_acres_0 = well_capacity_data[tot_acres == 0]
tot_acres_0 = tot_acres_0[data.table::like(CR, "FA") |
data.table::like(CR, "dry")]
tot_acres_325 = well_capacity_data[tot_acres == 32.5]
tot_acres_325 = tot_acres_325[(quarter == 1 & !(data.table::like(CR,
"FA") | data.table::like(CR, "dry"))) | (quarter !=
1 & (data.table::like(CR, "FA") | data.table::like(CR,
"dry")))]
tot_acres_65 = well_capacity_data[tot_acres == 65]
tot_acres_65 = tot_acres_65[(quarter < 3 & !(data.table::like(CR,
"FA") | data.table::like(CR, "dry"))) | (quarter >
2 & (data.table::like(CR, "FA") | data.table::like(CR,
"dry")))]
tot_acres_975 = well_capacity_data[tot_acres == 97.5]
tot_acres_975 = tot_acres_975[(quarter < 4 & !(data.table::like(CR,
"FA") | data.table::like(CR, "dry"))) | (quarter >
3 & (data.table::like(CR, "FA") | data.table::like(CR,
"dry")))]
tot_acres_130 = well_capacity_data[tot_acres == 130]
tot_acres_130 = tot_acres_130[!(data.table::like(CR,
"FA") | data.table::like(CR, "dry"))]
well_capacity_data = rbind(tot_acres_0, tot_acres_325,
tot_acres_65, tot_acres_975, tot_acres_130)
rm(tot_acres_0, tot_acres_325, tot_acres_65, tot_acres_975,
tot_acres_130)
data.table::setkey(well_capacity_data, Well_ID, tot_acres,
quarter)
KS_DSSAT[, `:=`(IFREQ, round(IFREQ, 1))]
number_of_years = nrow(unique(KS_DSSAT, by = "SDAT"))
foo = data.table::data.table(SOIL_ID = KS_DSSAT_2[group ==
i, unique(SOIL_ID)], WSTA = KS_DSSAT_2[group == i,
unique(WSTA)], CR = "FA", IFREQ = 0, PAW = soil_moisture_targets[1],
SDAT = min(KS_DSSAT$SDAT), irr_mm = 0, PRCP = 200,
PRCM = 200, yield_kg_ac = 0)
foo = foo[rep(seq_len(nrow(foo)), each = number_of_years)]
baz = unique(KS_DSSAT, by = "SDAT")
foo[, `:=`(SDAT, rep(baz$SDAT, nrow(foo)/nrow(baz)))]
foo[, `:=`(PRCP, rep(baz$PRCP, nrow(foo)/nrow(baz)))]
KS_DSSAT = rbind(KS_DSSAT, foo)
data.table::setkey(KS_DSSAT, SOIL_ID, WSTA, CR, IFREQ)
data.table::setkey(KS_DSSAT, CR)
KS_DSSAT = KS_DSSAT[price_dt]
KS_DSSAT[, cost_per_acre_in := unique(cost_dt$cost_per_acre_in)] # fix this one
data.table::setkey(KS_DSSAT, CR)
fixed_cost = unique(fixed_cost, by=colnames(fixed_cost))
data.table::setkey(fixed_cost, Crop)
KS_DSSAT = KS_DSSAT[fixed_cost]
KS_DSSAT[, `:=`(irrigation, 32.5 * irr_mm * 0.0393701)]
KS_DSSAT[, `:=`(profit, 32.5 * (yield_kg_ac * price -
f_cost) - irrigation * cost_per_acre_in)]
KS_DSSAT = KS_DSSAT[complete.cases(IFREQ)]
well_capacity_data[data.table::like(CR, "dry"), `:=`(ifreq,
0)]
well_capacity_data[CR == "FA", `:=`(ifreq, 0)]
cols = colnames(well_capacity_data)
well_capacity_data = unique(well_capacity_data, by = cols)
data.table::setkey(well_capacity_data, Well_ID, tot_acres, quarter)
well_capacity_data[, `:=`(Well_ID_grp, .GRP), by = "Well_ID"]
# well_capacity_data[, `:=`(group_1, .GRP), by = c("Well_ID", "tot_acres", "SDAT")]
# well_capacity_data[, `:=`(group_2, 1:.N), by = c("Well_ID", "tot_acres", "SDAT", "quarter")]
aa = max(well_capacity_data$Well_ID_grp)
foo_irr_2 = merge(well_capacity_data, KS_DSSAT, by.x = c("Soil_Type",
"weather_station", "CR", "ifreq"), by.y = c("SOIL_ID",
"WSTA", "CR", "IFREQ"), allow.cartesian = T)
foo_irr_2 = foo_irr_2[, .(Well_ID, SOIL_ID = Soil_Type, WSTA = weather_station,
Well_capacity, tot_acres, IFREQ = ifreq, CR, quarter,
PAW, SDAT, irr_mm, PRCP, PRCM, irrigation, yield_kg_ac, profit)]
if (Sys.info()[1] == "Windows") {
library(snow)
cl <- makeCluster(num_clusters)
print(Sys.info()[1])
parallel::clusterExport(cl, varlist = c("foo_irr", "data.table", ".", "aa", "FN_optim2",
"setnames", "setkey", "tax_amount",
"subsidy_threshold"), envir = environment())
foo_dt_all_1 <- parLapply(cl, 1:floor(aa/4), FN_optim2)
foo_dt_all_2 <- parLapply(cl, (floor(aa/4)+1):(2*floor(aa/4)), FN_optim2)
foo_dt_all_3 <- parLapply(cl, (2*floor(aa/4)+1):(3*floor(aa/4)), FN_optim2)
foo_dt_all_4 <- parLapply(cl, (3*floor(aa/4)+1):(4*floor(aa/4)), FN_optim2)
stopCluster(cl)
}
else {
print(Sys.info()[1])
foo_dt_all_1 <- mclapply(X = 1:floor(aa/4), FUN = FN_optim2,
mc.cores = num_clusters)
foo_dt_all_2 <- mclapply(X = (floor(aa/4)+1):(2*floor(aa/4)), FUN = FN_optim2,
mc.cores = num_clusters)
foo_dt_all_3 <- mclapply(X = (2*floor(aa/4)+1):(3*floor(aa/4)), FUN = FN_optim2,
mc.cores = num_clusters)
foo_dt_all_4 <- mclapply(X = (3*floor(aa/4)+1):aa, FUN = FN_optim2,
mc.cores = num_clusters)
}
foo_dt_all_1 <- do.call(rbind, foo_dt_all_1)
foo_dt_all_1 = data.table(foo_dt_all_1)
foo_dt_all_2 <- do.call(rbind, foo_dt_all_2)
foo_dt_all_2 = data.table(foo_dt_all_2)
foo_dt_all_3 <- do.call(rbind, foo_dt_all_3)
foo_dt_all_3 = data.table(foo_dt_all_3)
foo_dt_all_4 <- do.call(rbind, foo_dt_all_4)
foo_dt_all_4 = data.table(foo_dt_all_4)
foo_dt_all = rbind(foo_dt_all_1, foo_dt_all_2, foo_dt_all_3, foo_dt_all_4)
foo_dt_all = unique(foo_dt_all, by=c("Well_capacity", "quarter"))
foo_dt_all = foo_dt_all[, .(Well_capacity, tot_acres, quarter, ifreq,
CR, PAW, mean_irrigation_combination, mean_profit_combination,
mean_profit_combination_sub)]
setkey(foo_dt_all, Well_capacity, tot_acres, quarter, CR, PAW)
setkey(foo_irr_2, Well_capacity, tot_acres, quarter, CR, PAW)
lookup_table_quarter = foo_irr_2[foo_dt_all]
lookup_table_quarter = unique(lookup_table_quarter, by=colnames(lookup_table_quarter))
lookup_table_quarter = unique(lookup_table_quarter, by=c("SOIL_ID", "WSTA", "Well_capacity",
"tot_acres", "IFREQ", "CR", "quarter", "PAW", "SDAT"))
lookup_table_well = unique(lookup_table_quarter, by = "Well_capacity")
lookup_table_quarter = lookup_table_quarter[, .(Well_capacity,
SOIL_ID, WSTA, tot_acres, quarter, ifreq, SDAT, CR, PAW,
irrigation_quarter = irrigation, yield_kg_ac, profit_quarter = profit, #subsidy_amount
tax_amount)]
data.table::setkey(lookup_table_well, Well_capacity)
lookup_table_well = lookup_table_well[, .(Well_capacity,
SOIL_ID, WSTA, tot_acres, irr_tot_acres = mean_irrigation_combination,
profit_Well_ID = mean_profit_combination, profit_Well_ID_subsidy = mean_profit_combination_sub, #subsidy_amount
tax_amount)]
lookup_table_quarter = unique(lookup_table_quarter, by=colnames(lookup_table_quarter))
lookup_table_well = unique(lookup_table_well, by=colnames(lookup_table_well))
lookup_table_all_years = copy(lookup_table_quarter)
lookup_table_all_years[, `:=`(irr_tot_acres, sum(irrigation_quarter)),
by = c("Well_capacity", "SDAT")]
lookup_table_all_years[, `:=`(profit_Well_ID, sum(profit_quarter)),
by = c("Well_capacity", "SDAT")]
lookup_table_all_years[, `:=`(irr_below, ifelse(irr_tot_acres <
subsidy_threshold, subsidy_threshold - irr_tot_acres,
0))]
lookup_table_all_years = unique(lookup_table_all_years,
by = c("Well_capacity", "SDAT"))
lookup_table_all_years[, `:=`(profit_Well_ID_sub,
profit_Well_ID + irr_below * tax_amount)]
lookup_table_all_years = lookup_table_all_years[, .(Well_capacity,
SOIL_ID, WSTA, tot_acres, SDAT, irr_tot_acres, profit_Well_ID,
irr_below, profit_Well_ID_sub)]
lookup_table_all_years_2 = rbind(lookup_table_all_years_2,
lookup_table_all_years)
lookup_table_quarter_2 = rbind(lookup_table_quarter_2,
lookup_table_quarter)
lookup_table_well_2 = rbind(lookup_table_well_2, lookup_table_well)
print(i)
}
data.table::setkey(lookup_table_all_years_2, SOIL_ID, Well_capacity,
SDAT)
data.table::setkey(lookup_table_quarter_2, SOIL_ID, Well_capacity,
quarter)
data.table::setkey(lookup_table_well_2, SOIL_ID, Well_capacity)
write.csv(lookup_table_well_2, "lookup_table_well_2.csv")
saveRDS(lookup_table_quarter_2, "lookup_table_quarter_2.rds")
saveRDS(lookup_table_all_years_2, "lookup_table_all_years_2.rds")
}
manirouhirad/MODSSAT documentation built on April 15, 2024, 11:31 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions gen_lookup_tax_par_win_mac_cluster
Documented in gen_lookup_tax_par_win_mac_cluster
#' This function is the parallelized version of the gen_lookup_subsidy function.
#' @param tax_amount is the amount of subsidy per acre-inch of groundwater extracted. Defaults to 1.
#' @param subsidy_threshold is the threshold of subsidy, i.e., the amount of groundwater extraction above which subsidy is zero. Defaults to 400.
#' @param DSSAT_files is the directory where DSSAT files are located. Defaults to "C:/Users/manirad/Dropbox/DSSAT subregions work pre-2018 annual meeting/subregion KS files/outputs_for_econ/revised".
#' @param soil_file is the file that contains base soil file. Defaults to "C:/Users/manirad/Downloads/test/Well_Soil Type_generator_07.csv".
#' @param well_capacity_file is the file that contains base well capacity. Defaults to "C:/Users/manirad/Downloads/test/Well_Capacity_ganarator.csv".
#' @param price_file is the file that includes crop prices. Defaults to "C:/Users/manirad/Dropbox/DSSAT subregions work pre-2018 annual meeting/subregion KS files/crop_prices.csv".
#' @param fixed_cost_file is the file tha includes fixed costs (per acre) costs of production. Defaults to "C:/Users/manirad/Downloads/test/fixed_cost_input.csv".
#' @param pumping_cost is the cost of pumping an acre-inch of groundwater excluding taxes. Defaults to 3.56 which is the cost of pumpin an acre-inch of groundwater in parts of Kansas.
#' @param default_well_capacity_col_name is the name of the well capacity column generated from the MODFLOW model. Defaults to 'Well_Capacity(gpm)' as this is the original column name we started with.
#' @param soil_moisture_targets is the vector of soil moisture targets that are generated by the DSSAT model. (25, 35, 45, 55, 65, 75).
#' @param IFREQ_seq is the difference between two consequtive irrigation frequencies (IFREQ) in DSSAT. Defaults to 2.
#' @param IFREQ_interpolate is the size of interpolation interval. Defaults to 0.1 so that IFREQ_seq of 2 adds 0, 0.1, 0.2, ..., 1.9.
#' @param num_clusters is the number of cores for parallelization.
#' @return returns the output table.
#' @examples
#' \dontrun{
#' gen_lookup_tax(tax_amount = 2)
#' }
#' @export
gen_lookup_tax_par_win_mac_cluster = function(tax_amount = 1,
subsidy_threshold = 1,
DSSAT_files = "./input_files/DSSAT_files",
soil_file = "./input_files/Well_Soil Type_generator_07_1000.csv",
well_capacity_file = "./input_files/Well_Capacity_ganarator_1000.csv",
price_file = "./input_files/crop_prices.csv",
fixed_cost_file = "./input_files/fixed_cost_input.csv",
pumping_cost = 3.21,
default_well_capacity_col_name = "Well_Capacity(gpm)",
soil_moisture_targets = c(25, 35, 45, 55, 65, 75),
IFREQ_seq = 2,
IFREQ_interpolate = 0.1,
num_clusters = parallel::detectCores()-2
)
{
library(data.table)
library(parallel)
tax_amount = (tax_amount - 1)/10
print(paste("this is the marginal tax rate:", tax_amount,
"per acre-inch", sep = " "))
print(paste("this is the subsidy threshold:", subsidy_threshold,
"acre-inch", sep = " "))
FN_optim2 = function(jj) {
foo_irr_3 = well_capacity_data[Well_ID_grp == jj]
foo_irr_3[ifreq > KS_DSSAT[, max(IFREQ)], ifreq := KS_DSSAT[, max(IFREQ)]]
foo_irr_3_0 = foo_irr_3[ifreq == 0]
foo_irr_3_N0 = foo_irr_3[ifreq != 0]
tryCatch({
foo_irr_3_0 = merge(foo_irr_3_0, KS_DSSAT,
by.x = c("Soil_Type", "weather_station", "CR", "ifreq"),
by.y = c("SOIL_ID", "WSTA", "CR", "IFREQ"), allow.cartesian = T)
}, error=function(e){})
foo_irr_3_N0[, diff := (KS_DSSAT[, max(IFREQ)] - ifreq) * (1/IFREQ_interpolate)]
foo_irr_3_N0[, id := 1:.N]
tryCatch({
foo_irr_3_N0 = foo_irr_3_N0[rep(1:.N, (diff+1))]
foo_irr_3_N0[, foo := 1:.N, by=c("id")]
foo_irr_3_N0[, foo := foo - 1]
foo_irr_3_N0[, ifreq := ifreq + foo * IFREQ_interpolate]
}, error=function(e){})
foo_irr_3_N0 = merge(foo_irr_3_N0, KS_DSSAT,
by.x = c("Soil_Type", "weather_station", "CR", "ifreq"),
by.y = c("SOIL_ID", "WSTA", "CR", "IFREQ"), allow.cartesian = T)
foo_irr_3_N0[, c("id", "diff", "foo") := NULL]
foo_irr_3 = rbind(foo_irr_3_N0, foo_irr_3_0)
# foo_irr_3[, ifreq_denom := ifreq/.2]
# foo_irr_3 = foo_irr_3[ifreq_denom %in% 0:200]
foo_irr_3[ifreq == KS_DSSAT[, max(IFREQ)], profit := profit/4]
rm(foo_irr_3_N0, foo_irr_3_0)
foo_irr_3[, ifreq_cap := min(ifreq), by=c("Well_capacity", "SDAT", "tot_acres", "CR", "PAW")]
foo_irr_3[, ifreq_cap := ifreq_cap + 2]
foo_irr_3 = foo_irr_3[ifreq <= ifreq_cap]
# foo_irr_3[, id := .GRP, by=c("Soil_Type", "weather_station", "CR", "Well_ID", "Well_capacity", "tot_acres", "quarter", "PAW", "SDAT")]
foo_irr_3 = unique(foo_irr_3, by=colnames(foo_irr_3))
setkey(foo_irr_3, Well_capacity, tot_acres, CR, PAW, SDAT, ifreq, quarter)
foo_irr_3[, id := .GRP, by=c("Well_capacity", "tot_acres", "ifreq")]
#----
foo_dt1_000 = foo_irr_3[quarter == 1 & tot_acres == 0, .(Well_capacity, SDAT, tot_acres, ifreq_1 = ifreq, CR_1 = CR, PAW_1 = PAW, irrigation_1 = irrigation, profit_1 = profit)]
foo_dt2_000 = foo_irr_3[quarter == 2 & tot_acres == 0, .(Well_capacity, SDAT, tot_acres, ifreq_2 = ifreq, CR_2 = CR, PAW_2 = PAW, irrigation_2 = irrigation, profit_2 = profit)]
foo_dt3_000 = foo_irr_3[quarter == 3 & tot_acres == 0, .(Well_capacity, SDAT, tot_acres, ifreq_3 = ifreq, CR_3 = CR, PAW_3 = PAW, irrigation_3 = irrigation, profit_3 = profit)]
foo_dt4_000 = foo_irr_3[quarter == 4 & tot_acres == 0, .(Well_capacity, SDAT, tot_acres, ifreq_4 = ifreq, CR_4 = CR, PAW_4 = PAW, irrigation_4 = irrigation, profit_4 = profit)]
foo_dt3_000 = merge(foo_dt3_000, foo_dt4_000, by = c("Well_capacity", "SDAT", "tot_acres"),
allow.cartesian = T)
foo_dt2_000 = merge(foo_dt2_000, foo_dt3_000, by = c("Well_capacity", "SDAT", "tot_acres"),
allow.cartesian = T)
foo_dt1_000 = merge(foo_dt1_000, foo_dt2_000, by = c("Well_capacity", "SDAT", "tot_acres"),
allow.cartesian = T)
rm(foo_dt2_000, foo_dt3_000, foo_dt4_000)
foo_dt1_000[, `:=`(irrigation_sum, irrigation_1 + irrigation_2 + irrigation_3 + irrigation_4)]
# foo_dt1_000[, `:=`(irrigation_below, ifelse(irrigation_sum < subsidy_threshold, subsidy_threshold - irrigation_sum,
# 0))]
foo_dt1_000[, `:=`(profit_sum, profit_1 + profit_2 + profit_3 + profit_4)]
# foo_dt1_000[, `:=`(subsidy_payment, irrigation_below * subsidy_amount)]
# foo_dt1_000[, `:=`(profit_sum_sub, profit_sum + irrigation_below * subsidy_amount)]
foo_dt1_000[, `:=`(profit_sum_sub, profit_sum - irrigation_sum * tax_amount)]
foo_dt1_000[, row := .GRP, by=c("CR_1", "PAW_1",
"CR_2", "PAW_2",
"CR_3", "PAW_3",
"CR_4", "PAW_4")]
foo = foo_dt1_000[, .(row, profit_sum, profit_sum_sub, irrigation_sum)]
foo[, `:=`(mean_profit_combination, mean(profit_sum)), by = c("row")]
foo[, `:=`(mean_profit_combination_sub, mean(profit_sum_sub)), by = c("row")]
foo[, `:=`(mean_irrigation_combination, mean(irrigation_sum)), by = c("row")]
foo = unique(foo, by = c("row"))
foo[, `:=`(max_p, max(mean_profit_combination_sub))]
foo = foo[mean_profit_combination_sub == max_p,.(row, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
setkey(foo, row)
setkey(foo_dt1_000, row)
foo_dt1_000 = foo_dt1_000[foo]
foo_dt1_000 = rbind(foo_dt1_000[SDAT == min(foo_dt1_000$SDAT),.(Well_capacity, tot_acres, quarter = 1, ifreq = ifreq_1, CR = CR_1, PAW = PAW_1, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_000[SDAT == min(foo_dt1_000$SDAT),.(Well_capacity, tot_acres, quarter = 2, ifreq = ifreq_2, CR = CR_2, PAW = PAW_2, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_000[SDAT == min(foo_dt1_000$SDAT),.(Well_capacity, tot_acres, quarter = 3, ifreq = ifreq_3, CR = CR_3, PAW = PAW_3, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_000[SDAT == min(foo_dt1_000$SDAT),.(Well_capacity, tot_acres, quarter = 4, ifreq = ifreq_4, CR = CR_4, PAW = PAW_4, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
)
#----------
foo_dt1_325 = foo_irr_3[quarter == 1 & tot_acres == 32.5, .(Well_capacity, SDAT, tot_acres, ifreq_1 = ifreq, CR_1 = CR, PAW_1 = PAW, irrigation_1 = irrigation, profit_1 = profit)]
foo_dt2_325 = foo_irr_3[quarter == 2 & tot_acres == 32.5, .(Well_capacity, SDAT, tot_acres, ifreq_2 = ifreq, CR_2 = CR, PAW_2 = PAW, irrigation_2 = irrigation, profit_2 = profit)]
foo_dt3_325 = foo_irr_3[quarter == 3 & tot_acres == 32.5, .(Well_capacity, SDAT, tot_acres, ifreq_3 = ifreq, CR_3 = CR, PAW_3 = PAW, irrigation_3 = irrigation, profit_3 = profit)]
foo_dt4_325 = foo_irr_3[quarter == 4 & tot_acres == 32.5, .(Well_capacity, SDAT, tot_acres, ifreq_4 = ifreq, CR_4 = CR, PAW_4 = PAW, irrigation_4 = irrigation, profit_4 = profit)]
foo_dt3_325 = merge(foo_dt3_325, foo_dt4_325, by = c("Well_capacity", "SDAT", "tot_acres"),
allow.cartesian = T)
foo_dt2_325 = merge(foo_dt2_325, foo_dt3_325, by = c("Well_capacity", "SDAT", "tot_acres"),
allow.cartesian = T)
foo_dt1_325 = merge(foo_dt1_325, foo_dt2_325, by = c("Well_capacity", "SDAT", "tot_acres"),
allow.cartesian = T)
rm(foo_dt2_325, foo_dt3_325, foo_dt4_325)
foo_dt1_325[, `:=`(irrigation_sum, irrigation_1 + irrigation_2 + irrigation_3 + irrigation_4)]
# foo_dt1_325[, `:=`(irrigation_below, ifelse(irrigation_sum < subsidy_threshold, subsidy_threshold - irrigation_sum,
# 0))]
foo_dt1_325[, `:=`(profit_sum, profit_1 + profit_2 + profit_3 + profit_4)]
# foo_dt1_325[, `:=`(subsidy_payment, irrigation_below * subsidy_amount)]
# foo_dt1_325[, `:=`(profit_sum_sub, profit_sum + irrigation_below * subsidy_amount)]
foo_dt1_325[, `:=`(profit_sum_sub, profit_sum - irrigation_sum * tax_amount)]
foo_dt1_325[, row := .GRP, by=c("CR_1", "PAW_1",
"CR_2", "PAW_2",
"CR_3", "PAW_3",
"CR_4", "PAW_4", "ifreq_1")]
foo = foo_dt1_325[, .(row, profit_sum, profit_sum_sub, irrigation_sum)]
foo[, `:=`(mean_profit_combination, mean(profit_sum)), by = c("row")]
foo[, `:=`(mean_profit_combination_sub, mean(profit_sum_sub)), by = c("row")]
foo[, `:=`(mean_irrigation_combination, mean(irrigation_sum)), by = c("row")]
foo = unique(foo, by = c("row"))
foo[, `:=`(max_p, max(mean_profit_combination_sub))]
foo = foo[mean_profit_combination_sub == max_p,.(row, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
setkey(foo, row)
setkey(foo_dt1_325, row)
foo_dt1_325 = foo_dt1_325[foo]
foo_dt1_325 = rbind(foo_dt1_325[SDAT == min(foo_dt1_325$SDAT),.(Well_capacity, tot_acres, quarter = 1, ifreq = ifreq_1, CR = CR_1, PAW = PAW_1, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_325[SDAT == min(foo_dt1_325$SDAT),.(Well_capacity, tot_acres, quarter = 2, ifreq = ifreq_2, CR = CR_2, PAW = PAW_2, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_325[SDAT == min(foo_dt1_325$SDAT),.(Well_capacity, tot_acres, quarter = 3, ifreq = ifreq_3, CR = CR_3, PAW = PAW_3, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_325[SDAT == min(foo_dt1_325$SDAT),.(Well_capacity, tot_acres, quarter = 4, ifreq = ifreq_4, CR = CR_4, PAW = PAW_4, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
)
#----------
foo_dt1_650 = foo_irr_3[quarter == 1 & tot_acres == 65, .(Well_capacity, SDAT, tot_acres, ifreq_1 = ifreq, CR_1 = CR, PAW_1 = PAW, irrigation_1 = irrigation, profit_1 = profit, id)]
foo_dt2_650 = foo_irr_3[quarter == 2 & tot_acres == 65, .(Well_capacity, SDAT, tot_acres, ifreq_2 = ifreq, CR_2 = CR, PAW_2 = PAW, irrigation_2 = irrigation, profit_2 = profit, id)]
foo_dt3_650 = foo_irr_3[quarter == 3 & tot_acres == 65, .(Well_capacity, SDAT, tot_acres, ifreq_3 = ifreq, CR_3 = CR, PAW_3 = PAW, irrigation_3 = irrigation, profit_3 = profit)]
foo_dt4_650 = foo_irr_3[quarter == 4 & tot_acres == 65, .(Well_capacity, SDAT, tot_acres, ifreq_4 = ifreq, CR_4 = CR, PAW_4 = PAW, irrigation_4 = irrigation, profit_4 = profit)]
foo_dt3_650 = merge(foo_dt3_650, foo_dt4_650, by = c("Well_capacity", "SDAT", "tot_acres"),
allow.cartesian = T)
foo_dt2_650 = merge(foo_dt2_650, foo_dt3_650, by = c("Well_capacity", "SDAT", "tot_acres"),
allow.cartesian = T)
foo_dt1_650 = merge(foo_dt1_650, foo_dt2_650, by = c("Well_capacity", "SDAT", "tot_acres", "id"),
allow.cartesian = T)
rm(foo_dt2_650, foo_dt3_650, foo_dt4_650)
foo_dt1_650[, `:=`(irrigation_sum, irrigation_1 + irrigation_2 + irrigation_3 + irrigation_4)]
# foo_dt1_650[, `:=`(irrigation_below, ifelse(irrigation_sum < subsidy_threshold, subsidy_threshold - irrigation_sum,
# 0))]
foo_dt1_650[, `:=`(profit_sum, profit_1 + profit_2 + profit_3 + profit_4)]
# foo_dt1_650[, `:=`(subsidy_payment, irrigation_below * subsidy_amount)]
# foo_dt1_650[, `:=`(profit_sum_sub, profit_sum + irrigation_below * subsidy_amount)]
foo_dt1_650[, `:=`(profit_sum_sub, profit_sum - irrigation_sum * tax_amount)]
foo_dt1_650[, row := .GRP, by=c("CR_1", "PAW_1",
"CR_2", "PAW_2",
"CR_3", "PAW_3",
"CR_4", "PAW_4", "id")]
foo = foo_dt1_650[, .(row, profit_sum, profit_sum_sub, irrigation_sum)]
foo[, `:=`(mean_profit_combination, mean(profit_sum)), by = c("row")]
foo[, `:=`(mean_profit_combination_sub, mean(profit_sum_sub)), by = c("row")]
foo[, `:=`(mean_irrigation_combination, mean(irrigation_sum)), by = c("row")]
foo = unique(foo, by = c("row"))
foo[, `:=`(max_p, max(mean_profit_combination_sub))]
foo = foo[mean_profit_combination_sub == max_p,.(row, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
setkey(foo, row)
setkey(foo_dt1_650, row)
foo_dt1_650 = foo_dt1_650[foo]
foo_dt1_650 = rbind(foo_dt1_650[SDAT == min(foo_dt1_650$SDAT),.(Well_capacity, tot_acres, quarter = 1, ifreq = ifreq_1, CR = CR_1, PAW = PAW_1, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_650[SDAT == min(foo_dt1_650$SDAT),.(Well_capacity, tot_acres, quarter = 2, ifreq = ifreq_2, CR = CR_2, PAW = PAW_2, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_650[SDAT == min(foo_dt1_650$SDAT),.(Well_capacity, tot_acres, quarter = 3, ifreq = ifreq_3, CR = CR_3, PAW = PAW_3, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_650[SDAT == min(foo_dt1_650$SDAT),.(Well_capacity, tot_acres, quarter = 4, ifreq = ifreq_4, CR = CR_4, PAW = PAW_4, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
)
#----------
foo_dt1_975 = foo_irr_3[quarter == 1 & tot_acres == 97.5, .(Well_capacity, SDAT, tot_acres, ifreq_1 = ifreq, CR_1 = CR, PAW_1 = PAW, irrigation_1 = irrigation, profit_1 = profit, id)]
foo_dt2_975 = foo_irr_3[quarter == 2 & tot_acres == 97.5, .(Well_capacity, SDAT, tot_acres, ifreq_2 = ifreq, CR_2 = CR, PAW_2 = PAW, irrigation_2 = irrigation, profit_2 = profit, id)]
foo_dt3_975 = foo_irr_3[quarter == 3 & tot_acres == 97.5, .(Well_capacity, SDAT, tot_acres, ifreq_3 = ifreq, CR_3 = CR, PAW_3 = PAW, irrigation_3 = irrigation, profit_3 = profit, id)]
foo_dt4_975 = foo_irr_3[quarter == 4 & tot_acres == 97.5, .(Well_capacity, SDAT, tot_acres, ifreq_4 = ifreq, CR_4 = CR, PAW_4 = PAW, irrigation_4 = irrigation, profit_4 = profit)]
foo_dt3_975 = merge(foo_dt3_975, foo_dt4_975, by = c("Well_capacity", "SDAT", "tot_acres"),
allow.cartesian = T)
foo_dt2_975 = merge(foo_dt2_975, foo_dt3_975, by = c("Well_capacity", "SDAT", "tot_acres", "id"),
allow.cartesian = T)
foo_dt1_975 = merge(foo_dt1_975, foo_dt2_975, by = c("Well_capacity", "SDAT", "tot_acres", "id"),
allow.cartesian = T)
rm(foo_dt2_975, foo_dt3_975, foo_dt4_975)
foo_dt1_975[, `:=`(irrigation_sum, irrigation_1 + irrigation_2 + irrigation_3 + irrigation_4)]
# foo_dt1_975[, `:=`(irrigation_below, ifelse(irrigation_sum < subsidy_threshold, subsidy_threshold - irrigation_sum,
# 0))]
foo_dt1_975[, `:=`(profit_sum, profit_1 + profit_2 + profit_3 + profit_4)]
# foo_dt1_975[, `:=`(subsidy_payment, irrigation_below * subsidy_amount)]
# foo_dt1_975[, `:=`(profit_sum_sub, profit_sum + irrigation_below * subsidy_amount)]
foo_dt1_975[, `:=`(profit_sum_sub, profit_sum - irrigation_sum * tax_amount)]
foo_dt1_975[, row := .GRP, by=c("CR_1", "PAW_1",
"CR_2", "PAW_2",
"CR_3", "PAW_3",
"CR_4", "PAW_4", "id")]
foo = foo_dt1_975[, .(row, profit_sum, profit_sum_sub, irrigation_sum)]
foo[, `:=`(mean_profit_combination, mean(profit_sum)), by = c("row")]
foo[, `:=`(mean_profit_combination_sub, mean(profit_sum_sub)), by = c("row")]
foo[, `:=`(mean_irrigation_combination, mean(irrigation_sum)), by = c("row")]
foo = unique(foo, by = c("row"))
foo[, `:=`(max_p, max(mean_profit_combination_sub))]
foo = foo[mean_profit_combination_sub == max_p,.(row, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
setkey(foo, row)
setkey(foo_dt1_975, row)
foo_dt1_975 = foo_dt1_975[foo]
foo_dt1_975 = rbind(foo_dt1_975[SDAT == min(foo_dt1_975$SDAT),.(Well_capacity, tot_acres, quarter = 1, ifreq = ifreq_1, CR = CR_1, PAW = PAW_1, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_975[SDAT == min(foo_dt1_975$SDAT),.(Well_capacity, tot_acres, quarter = 2, ifreq = ifreq_2, CR = CR_2, PAW = PAW_2, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_975[SDAT == min(foo_dt1_975$SDAT),.(Well_capacity, tot_acres, quarter = 3, ifreq = ifreq_3, CR = CR_3, PAW = PAW_3, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_975[SDAT == min(foo_dt1_975$SDAT),.(Well_capacity, tot_acres, quarter = 4, ifreq = ifreq_4, CR = CR_4, PAW = PAW_4, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
)
#----------
foo_dt1_130 = foo_irr_3[quarter == 1 & tot_acres == 130, .(Well_capacity, SDAT, tot_acres, ifreq_1 = ifreq, CR_1 = CR, PAW_1 = PAW, irrigation_1 = irrigation, profit_1 = profit, id)]
foo_dt2_130 = foo_irr_3[quarter == 2 & tot_acres == 130, .(Well_capacity, SDAT, tot_acres, ifreq_2 = ifreq, CR_2 = CR, PAW_2 = PAW, irrigation_2 = irrigation, profit_2 = profit, id)]
foo_dt3_130 = foo_irr_3[quarter == 3 & tot_acres == 130, .(Well_capacity, SDAT, tot_acres, ifreq_3 = ifreq, CR_3 = CR, PAW_3 = PAW, irrigation_3 = irrigation, profit_3 = profit, id)]
foo_dt4_130 = foo_irr_3[quarter == 4 & tot_acres == 130, .(Well_capacity, SDAT, tot_acres, ifreq_4 = ifreq, CR_4 = CR, PAW_4 = PAW, irrigation_4 = irrigation, profit_4 = profit, id)]
foo_dt3_130 = merge(foo_dt3_130, foo_dt4_130, by = c("Well_capacity", "SDAT", "tot_acres", "id"),
allow.cartesian = T)
foo_dt2_130 = merge(foo_dt2_130, foo_dt3_130, by = c("Well_capacity", "SDAT", "tot_acres", "id"),
allow.cartesian = T)
foo_dt1_130 = merge(foo_dt1_130, foo_dt2_130, by = c("Well_capacity", "SDAT", "tot_acres", "id"),
allow.cartesian = T)
rm(foo_dt2_130, foo_dt3_130, foo_dt4_130)
foo_dt1_130[, `:=`(irrigation_sum, irrigation_1 + irrigation_2 + irrigation_3 + irrigation_4)]
# foo_dt1_130[, `:=`(irrigation_below, ifelse(irrigation_sum < subsidy_threshold, subsidy_threshold - irrigation_sum,
# 0))]
foo_dt1_130[, `:=`(profit_sum, profit_1 + profit_2 + profit_3 + profit_4)]
# foo_dt1_130[, `:=`(subsidy_payment, irrigation_below * subsidy_amount)]
# foo_dt1_130[, `:=`(profit_sum_sub, profit_sum + irrigation_below * subsidy_amount)]
foo_dt1_130[, `:=`(profit_sum_sub, profit_sum - irrigation_sum * tax_amount)]
foo_dt1_130[, row := .GRP, by=c("CR_1", "PAW_1",
"CR_2", "PAW_2",
"CR_3", "PAW_3",
"CR_4", "PAW_4", "id")]
foo = foo_dt1_130[, .(row, profit_sum, profit_sum_sub, irrigation_sum)]
foo[, `:=`(mean_profit_combination, mean(profit_sum)), by = c("row")]
foo[, `:=`(mean_profit_combination_sub, mean(profit_sum_sub)), by = c("row")]
foo[, `:=`(mean_irrigation_combination, mean(irrigation_sum)), by = c("row")]
foo = unique(foo, by = c("row"))
foo[, `:=`(max_p, max(mean_profit_combination_sub))]
foo = foo[mean_profit_combination_sub == max_p,.(row, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
setkey(foo, row)
setkey(foo_dt1_130, row)
foo_dt1_130 = foo_dt1_130[foo]
foo_dt1_130 = rbind(foo_dt1_130[SDAT == min(foo_dt1_130$SDAT),.(Well_capacity, tot_acres, quarter = 1, ifreq = ifreq_1, CR = CR_1, PAW = PAW_1, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_130[SDAT == min(foo_dt1_130$SDAT),.(Well_capacity, tot_acres, quarter = 2, ifreq = ifreq_2, CR = CR_2, PAW = PAW_2, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_130[SDAT == min(foo_dt1_130$SDAT),.(Well_capacity, tot_acres, quarter = 3, ifreq = ifreq_3, CR = CR_3, PAW = PAW_3, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
, foo_dt1_130[SDAT == min(foo_dt1_130$SDAT),.(Well_capacity, tot_acres, quarter = 4, ifreq = ifreq_4, CR = CR_4, PAW = PAW_4, mean_irrigation_combination, mean_profit_combination, mean_profit_combination_sub)]
)
#----------
foo_dt1 = rbind(foo_dt1_000, foo_dt1_325, foo_dt1_650, foo_dt1_975, foo_dt1_130)
rm(foo_dt1_000, foo_dt1_325, foo_dt1_650, foo_dt1_975, foo_dt1_130)
foo_dt1[, max_mean_p := max(mean_profit_combination_sub, na.rm = T)]
foo_dt1 = foo_dt1[max_mean_p == mean_profit_combination_sub]
return(foo_dt1)
}
col_new = c("RUNNO", "TRNO", "R_pound", "O_pound", "C_pound",
"CR", "MODEL", "EXNAME", "FNAM", "WSTA", "SOIL_ID", "SDAT",
"PDAT", "EDAT", "ADAT", "MDAT", "HDAT", "DWAP", "CWAM",
"HWAM", "HWAH", "BWAH", "PWAM", "HWUM", "H_AM", "H_UM",
"HIAM", "LAIX", "IR_M", "IRCM", "PRCM", "ETCM", "EPCM",
"ESCM", "ROCM", "DRCM", "SWXM", "NI_M", "NICM", "NFXM",
"NUCM", "NLCM", "NIAM", "CNAM", "GNAM", "N2OEC", "PI_M",
"PICM", "PUPC", "SPAM", "KI_M", "KICM", "KUPC", "SKAM",
"RECM", "ONTAM", "ONAM", "OPTAM", "OPAM", "OCTAM", "OCAM",
"CO2EC", "DMPPM", "DMPEM", "DMPTM", "DMPIM", "YPPM",
"YPEM", "YPTM", "YPIM", "DPNAM", "DPNUM", "YPNAM", "YPNUM",
"NDCH", "TMAXA", "TMINA", "SRADA", "DAYLA", "CO2A", "PRCP",
"ETCP", "ESCP", "EPCP", "PAW", "IFREQ")
filenames = list.files(path = DSSAT_files, pattern = "*.OSU",
full.names = TRUE, recursive = TRUE)
ldf <- lapply(filenames, read.table, fill = T)
ldf <- lapply(ldf, function(x) x[-2, ])
ldf <- mapply(cbind, ldf, file_name = filenames, SIMPLIFY = F)
KS_DSSAT = data.table::rbindlist(ldf)
KS_DSSAT = data.table(KS_DSSAT)
KS_DSSAT = KS_DSSAT[complete.cases(V85)]
KS_DSSAT[, `:=`(file_name, NULL)]
rm(ldf)
KS_DSSAT[, `:=`(foo, nchar(as.character(V9)))]
KS_DSSAT[foo < 4, `:=`(V9, paste(V9, V11, V10, sep = "_"))]
KS_DSSAT[foo < 4, `:=`(V10, NA)]
KS_DSSAT[foo < 4, `:=`(V11, NA)]
KS_DSSAT = KS_DSSAT[, !sapply(KS_DSSAT, function(x) all(is.na(x))),
with = FALSE]
KS_DSSAT[, `:=`(V9, as.character(V9))]
KS_DSSAT[, `:=`(PAW, substr(V9, 1, regexpr("_", V9) - 2))]
KS_DSSAT[, `:=`(IFREQ, substr(V9, regexpr("Q", V9) + 2, nchar(V9)))]
KS_DSSAT[, `:=`(V9, NULL)]
KS_DSSAT[, `:=`(foo, NULL)]
data.table::setnames(KS_DSSAT, old = colnames(KS_DSSAT),
new = col_new)
KS_DSSAT[, `:=`(WSTA, substr(WSTA, 1, 4))]
KS_DSSAT[, `:=`(group, .GRP), by = c("WSTA", "SOIL_ID")]
KS_DSSAT_2 = data.table::copy(KS_DSSAT)
lookup_table_all_years_2 = data.table::data.table()
lookup_table_quarter_2 = data.table::data.table()
lookup_table_well_2 = data.table::data.table()
for (i in 1:max(KS_DSSAT_2$group)) {
soil_type = data.table::fread(soil_file)
soil_type[, `:=`(Soil_Type, KS_DSSAT_2[group == i, unique(SOIL_ID)])]
soil_type[, `:=`(Soil_Type, gsub("KSFC00000", "KS0000000",
Soil_Type))]
well_capacity = data.table::fread(well_capacity_file)
well_capacity = well_capacity[`Well_Capacity(gpm)` <= 801]
data.table::setkey(soil_type, Well_ID)
data.table::setkey(well_capacity, Well_ID)
well_capacity_data = soil_type[well_capacity]
data.table::setnames(well_capacity_data, old = default_well_capacity_col_name,
"Well_capacity")
WSTAT = data.table(Well_ID = well_capacity_data$Well_ID,
weather_station = KS_DSSAT_2[group == i, unique(WSTA)])
data.table::setkey(WSTAT, Well_ID)
data.table::setkey(well_capacity_data, Well_ID)
well_capacity_data = WSTAT[well_capacity_data]
well_capacity_data = well_capacity_data[, .(Well_ID,
Soil_Type, weather_station, Well_capacity)]
price_dt = data.table::fread(price_file)
data.table::setkey(price_dt, CR)
cost_dt = well_capacity_data[, .(Well_ID)]
data.table::setkey(cost_dt, Well_ID)
cost_dt[, `:=`(cost_per_acre_in, (pumping_cost)/1)]
fixed_cost = data.table::fread(fixed_cost_file)
fixed_cost[irr == 0, `:=`(Crop, paste("dry", Crop, sep = "-"))]
fixed_cost[, `:=`(irr, NULL)]
fixed_cost = rbind(fixed_cost, data.table::data.table(Crop = "FA",
f_cost = 0))
data.table::setkey(fixed_cost, Crop)
KS_DSSAT = KS_DSSAT_2[group == i, .(SOIL_ID, WSTA, CR,
IFREQ, PAW, SDAT, IRCM, PRCP, PRCM, HWAM)]
KS_DSSAT[, `:=`(IRCM, as.numeric(as.character(IRCM)))]
KS_DSSAT = KS_DSSAT[complete.cases(IRCM)]
KS_DSSAT[, `:=`(SDAT, substr(SDAT, 1, 4))]
cols_change = c("IFREQ", "PAW", "SDAT", "IRCM", "PRCP",
"PRCM", "HWAM")
KS_DSSAT[, `:=`((cols_change), lapply(.SD, function(x) as.numeric(as.character(x)))),
.SDcols = cols_change]
cols_change = c("SOIL_ID", "WSTA", "CR")
KS_DSSAT[, `:=`((cols_change), lapply(.SD, as.character)),
.SDcols = cols_change]
KS_DSSAT = KS_DSSAT[IFREQ < 17]
qux1 = KS_DSSAT[IFREQ == 6 & PAW == 75]
qux2 = KS_DSSAT[IFREQ == 8 & PAW == 75]
qux1[, `:=`(IFREQ, 8)]
qux2[, `:=`(IFREQ, 6)]
KS_DSSAT = KS_DSSAT[!((IFREQ == 6 | IFREQ == 8) & PAW ==
75)]
KS_DSSAT = rbind(KS_DSSAT, qux1)
KS_DSSAT = rbind(KS_DSSAT, qux2)
KS_DSSAT = KS_DSSAT[PAW %in% soil_moisture_targets |
IFREQ == 0]
KS_DSSAT_0 = KS_DSSAT[IFREQ == 0]
KS_DSSAT = KS_DSSAT[IFREQ != 0]
KS_DSSAT_0 <- KS_DSSAT_0[rep(seq_len(nrow(KS_DSSAT_0)),
each = length(soil_moisture_targets))]
KS_DSSAT_0[, `:=`(PAW, rep(soil_moisture_targets, nrow(KS_DSSAT_0)/length(soil_moisture_targets)))]
KS_DSSAT = rbind(KS_DSSAT, KS_DSSAT_0)
data.table::setkey(KS_DSSAT, SOIL_ID, WSTA, CR, IFREQ,
PAW, SDAT)
KS_DSSAT = KS_DSSAT[PRCP > -10 & HWAM > -10]
KS_DSSAT = KS_DSSAT[, .(SOIL_ID, WSTA, CR, PAW, SDAT,
IFREQ, IRCM, PRCP, PRCM, HWAM)]
data.table::setkey(KS_DSSAT, SOIL_ID, WSTA, CR, PAW,
SDAT, IFREQ)
# foo = copy(KS_DSSAT)
# foo[, mean_HWAM := mean(HWAM), by=c("SOIL_ID", "WSTA", "CR", "IFREQ", "PAW")]
# foo[, mean_IRCM := mean(IRCM), by=c("SOIL_ID", "WSTA", "CR", "IFREQ", "PAW")]
# foo = unique(foo, by=c("SOIL_ID", "WSTA", "CR", "IFREQ", "PAW"))
#
# ggplot(foo[IFREQ == 4], aes(x=mean_IRCM, y= mean_HWAM, group=CR, col=CR))+
# geom_point()+
# geom_line()
# ggplot(KS_DSSAT[IFREQ == 4], aes(x=IRCM, y= HWAM, group=CR, col=CR))+
# geom_point()+
# geom_smooth()
#
# ggplot(KS_DSSAT[CR == "MZ"], aes(x=IRCM, y= HWAM, group=factor(IFREQ), col=factor(IFREQ)))+
# geom_point()
KS_DSSAT[, `:=`(lead_yield, dplyr::lead(HWAM, n = 1L)),
by = c("SOIL_ID", "WSTA", "CR", "PAW", "SDAT")]
KS_DSSAT[, `:=`(lead_irr_mm, dplyr::lead(IRCM, n = 1L)),
by = c("SOIL_ID", "WSTA", "CR", "PAW", "SDAT")]
KS_DSSAT <- KS_DSSAT[rep(seq_len(nrow(KS_DSSAT)), each = IFREQ_seq *
10)]
df_foo = data.table::data.table(x = seq(0, IFREQ_seq -
0.1, IFREQ_interpolate))
df_foo = do.call(rbind, replicate(nrow(KS_DSSAT)/nrow(df_foo),
df_foo, simplify = F))
KS_DSSAT = data.table::data.table(KS_DSSAT, IFREQ_int = df_foo$x)
KS_DSSAT[, `:=`(foo, max(IFREQ)), by = c("SOIL_ID", "WSTA",
"CR")]
KS_DSSAT[, `:=`(IFREQ, IFREQ + IFREQ_int)]
KS_DSSAT = KS_DSSAT[IFREQ <= foo]
KS_DSSAT[, `:=`(foo, NULL)]
KS_DSSAT = KS_DSSAT[complete.cases(lead_yield)]
KS_DSSAT[, `:=`(yield_int, HWAM + (lead_yield - HWAM)/IFREQ_seq *
IFREQ_int)]
KS_DSSAT[, `:=`(irr_int, IRCM + (lead_irr_mm - IRCM)/IFREQ_seq *
IFREQ_int)]
KS_DSSAT = KS_DSSAT[, .(SOIL_ID, WSTA, CR, IFREQ, PAW,
SDAT, irr_mm = irr_int, PRCP, PRCM, yield_kg_ac = yield_int)]
KS_DSSAT[, `:=`(yield_kg_ac, yield_kg_ac * 0.4046)]
data.table::setkey(KS_DSSAT, SOIL_ID, WSTA, CR, IFREQ,
PAW, SDAT)
KS_DSSAT = KS_DSSAT[IFREQ == 0 | IFREQ >= IFREQ_seq]
KS_DSSAT = KS_DSSAT[IFREQ != 0 | PAW == soil_moisture_targets[1]]
KS_DSSAT = KS_DSSAT[SDAT < 2016]
KS_DSSAT[IFREQ == 0, `:=`(CR, paste("dry", CR, sep = "-"))]
KS_DSSAT[CR %like% "MZ", `:=`(yield_kg_ac, yield_kg_ac *
1.183)]
KS_DSSAT[CR %like% "SG", `:=`(yield_kg_ac, yield_kg_ac *
1.183)]
KS_DSSAT[CR %like% "WH", `:=`(yield_kg_ac, yield_kg_ac *
1.156)]
number_of_crops = length(KS_DSSAT[, unique(CR)])
well_capacity_data = rbind(data.table::data.table(Well_ID = 1,
Soil_Type = unique(well_capacity_data$Soil_Type),
weather_station = unique(well_capacity_data$weather_station),
Well_capacity = 0), well_capacity_data)
well_capacity_data[, `:=`(quarter_1, 0)]
well_capacity_data[, `:=`(quarter_2, 0)]
well_capacity_data[, `:=`(quarter_3, 0)]
well_capacity_data[, `:=`(quarter_4, 0)]
well_capacity_data <- well_capacity_data[rep(seq_len(nrow(well_capacity_data)),
each = 5)]
well_capacity_data[, `:=`(tot_acres, rep(c(0, 32.5, 65,
97.5, 130), nrow(well_capacity_data)/5))]
well_capacity_data <- well_capacity_data[rep(seq_len(nrow(well_capacity_data)),
each = (number_of_crops + 1))]
well_capacity_data[, `:=`(quarter_4, rep(0:number_of_crops,
nrow(well_capacity_data)/(number_of_crops + 1)))]
well_capacity_data[, `:=`(quarter_3, rep(0:number_of_crops,
nrow(well_capacity_data)/(number_of_crops + 1)))]
well_capacity_data[, `:=`(quarter_2, rep(0:number_of_crops,
nrow(well_capacity_data)/(number_of_crops + 1)))]
well_capacity_data[, `:=`(quarter_1, rep(0:number_of_crops,
nrow(well_capacity_data)/(number_of_crops + 1)))]
well_capacity_data[, `:=`(ifreq, round((tot_acres * 1)/(Well_capacity *
0.053030303149011), 1))]
# well_capacity_data[!(complete.cases(ifreq) & ifreq < 100), `:=`(ifreq, 0)]
well_capacity_data[!complete.cases(ifreq), `:=`(ifreq, 0)]
well_capacity_data[ifreq < 2 & ifreq > 0, `:=`(ifreq,
2)]
well_capacity_data = reshape2::melt(well_capacity_data,
id = c("Well_ID", "Soil_Type", "weather_station",
"Well_capacity", "tot_acres", "ifreq"))
well_capacity_data = data.table::data.table(well_capacity_data)
data.table::setkey(well_capacity_data, Well_ID, Soil_Type,
weather_station, tot_acres)
quz = data.table(CR = KS_DSSAT[, unique(CR)])
foo = data.table(value = 0:number_of_crops, rbind(data.table(CR = "FA"),
quz))
setkey(well_capacity_data, value)
setkey(foo, value)
well_capacity_data = well_capacity_data[foo]
well_capacity_data[, `:=`(value, NULL)]
well_capacity_data[, `:=`(quarter, ifelse(variable ==
"quarter_1", 1, ifelse(variable == "quarter_2", 2,
ifelse(variable == "quarter_3", 3, 4))))]
well_capacity_data[, `:=`(variable, NULL)]
tot_acres_0 = well_capacity_data[tot_acres == 0]
tot_acres_0 = tot_acres_0[data.table::like(CR, "FA") |
data.table::like(CR, "dry")]
tot_acres_325 = well_capacity_data[tot_acres == 32.5]
tot_acres_325 = tot_acres_325[(quarter == 1 & !(data.table::like(CR,
"FA") | data.table::like(CR, "dry"))) | (quarter !=
1 & (data.table::like(CR, "FA") | data.table::like(CR,
"dry")))]
tot_acres_65 = well_capacity_data[tot_acres == 65]
tot_acres_65 = tot_acres_65[(quarter < 3 & !(data.table::like(CR,
"FA") | data.table::like(CR, "dry"))) | (quarter >
2 & (data.table::like(CR, "FA") | data.table::like(CR,
"dry")))]
tot_acres_975 = well_capacity_data[tot_acres == 97.5]
tot_acres_975 = tot_acres_975[(quarter < 4 & !(data.table::like(CR,
"FA") | data.table::like(CR, "dry"))) | (quarter >
3 & (data.table::like(CR, "FA") | data.table::like(CR,
"dry")))]
tot_acres_130 = well_capacity_data[tot_acres == 130]
tot_acres_130 = tot_acres_130[!(data.table::like(CR,
"FA") | data.table::like(CR, "dry"))]
well_capacity_data = rbind(tot_acres_0, tot_acres_325,
tot_acres_65, tot_acres_975, tot_acres_130)
rm(tot_acres_0, tot_acres_325, tot_acres_65, tot_acres_975,
tot_acres_130)
data.table::setkey(well_capacity_data, Well_ID, tot_acres,
quarter)
KS_DSSAT[, `:=`(IFREQ, round(IFREQ, 1))]
number_of_years = nrow(unique(KS_DSSAT, by = "SDAT"))
foo = data.table::data.table(SOIL_ID = KS_DSSAT_2[group ==
i, unique(SOIL_ID)], WSTA = KS_DSSAT_2[group == i,
unique(WSTA)], CR = "FA", IFREQ = 0, PAW = soil_moisture_targets[1],
SDAT = min(KS_DSSAT$SDAT), irr_mm = 0, PRCP = 200,
PRCM = 200, yield_kg_ac = 0)
foo = foo[rep(seq_len(nrow(foo)), each = number_of_years)]
baz = unique(KS_DSSAT, by = "SDAT")
foo[, `:=`(SDAT, rep(baz$SDAT, nrow(foo)/nrow(baz)))]
foo[, `:=`(PRCP, rep(baz$PRCP, nrow(foo)/nrow(baz)))]
KS_DSSAT = rbind(KS_DSSAT, foo)
data.table::setkey(KS_DSSAT, SOIL_ID, WSTA, CR, IFREQ)
data.table::setkey(KS_DSSAT, CR)
KS_DSSAT = KS_DSSAT[price_dt]
KS_DSSAT[, cost_per_acre_in := unique(cost_dt$cost_per_acre_in)] # fix this one
data.table::setkey(KS_DSSAT, CR)
fixed_cost = unique(fixed_cost, by=colnames(fixed_cost))
data.table::setkey(fixed_cost, Crop)
KS_DSSAT = KS_DSSAT[fixed_cost]
KS_DSSAT[, `:=`(irrigation, 32.5 * irr_mm * 0.0393701)]
KS_DSSAT[, `:=`(profit, 32.5 * (yield_kg_ac * price -
f_cost) - irrigation * cost_per_acre_in)]
KS_DSSAT = KS_DSSAT[complete.cases(IFREQ)]
well_capacity_data[data.table::like(CR, "dry"), `:=`(ifreq,
0)]
well_capacity_data[CR == "FA", `:=`(ifreq, 0)]
cols = colnames(well_capacity_data)
well_capacity_data = unique(well_capacity_data, by = cols)
data.table::setkey(well_capacity_data, Well_ID, tot_acres, quarter)
well_capacity_data[, `:=`(Well_ID_grp, .GRP), by = "Well_ID"]
# well_capacity_data[, `:=`(group_1, .GRP), by = c("Well_ID", "tot_acres", "SDAT")]
# well_capacity_data[, `:=`(group_2, 1:.N), by = c("Well_ID", "tot_acres", "SDAT", "quarter")]
aa = max(well_capacity_data$Well_ID_grp)
foo_irr_2 = merge(well_capacity_data, KS_DSSAT, by.x = c("Soil_Type",
"weather_station", "CR", "ifreq"), by.y = c("SOIL_ID",
"WSTA", "CR", "IFREQ"), allow.cartesian = T)
foo_irr_2 = foo_irr_2[, .(Well_ID, SOIL_ID = Soil_Type, WSTA = weather_station,
Well_capacity, tot_acres, IFREQ = ifreq, CR, quarter,
PAW, SDAT, irr_mm, PRCP, PRCM, irrigation, yield_kg_ac, profit)]
if (Sys.info()[1] == "Windows") {
library(snow)
cl <- makeCluster(num_clusters)
print(Sys.info()[1])
parallel::clusterExport(cl, varlist = c("foo_irr", "data.table", ".", "aa", "FN_optim2",
"setnames", "setkey", "tax_amount",
"subsidy_threshold"), envir = environment())
foo_dt_all_1 <- parLapply(cl, 1:floor(aa/4), FN_optim2)
foo_dt_all_2 <- parLapply(cl, (floor(aa/4)+1):(2*floor(aa/4)), FN_optim2)
foo_dt_all_3 <- parLapply(cl, (2*floor(aa/4)+1):(3*floor(aa/4)), FN_optim2)
foo_dt_all_4 <- parLapply(cl, (3*floor(aa/4)+1):(4*floor(aa/4)), FN_optim2)
stopCluster(cl)
}
else {
print(Sys.info()[1])
foo_dt_all_1 <- mclapply(X = 1:floor(aa/4), FUN = FN_optim2,
mc.cores = num_clusters)
foo_dt_all_2 <- mclapply(X = (floor(aa/4)+1):(2*floor(aa/4)), FUN = FN_optim2,
mc.cores = num_clusters)
foo_dt_all_3 <- mclapply(X = (2*floor(aa/4)+1):(3*floor(aa/4)), FUN = FN_optim2,
mc.cores = num_clusters)
foo_dt_all_4 <- mclapply(X = (3*floor(aa/4)+1):aa, FUN = FN_optim2,
mc.cores = num_clusters)
}
foo_dt_all_1 <- do.call(rbind, foo_dt_all_1)
foo_dt_all_1 = data.table(foo_dt_all_1)
foo_dt_all_2 <- do.call(rbind, foo_dt_all_2)
foo_dt_all_2 = data.table(foo_dt_all_2)
foo_dt_all_3 <- do.call(rbind, foo_dt_all_3)
foo_dt_all_3 = data.table(foo_dt_all_3)
foo_dt_all_4 <- do.call(rbind, foo_dt_all_4)
foo_dt_all_4 = data.table(foo_dt_all_4)
foo_dt_all = rbind(foo_dt_all_1, foo_dt_all_2, foo_dt_all_3, foo_dt_all_4)
foo_dt_all = unique(foo_dt_all, by=c("Well_capacity", "quarter"))
foo_dt_all = foo_dt_all[, .(Well_capacity, tot_acres, quarter, ifreq,
CR, PAW, mean_irrigation_combination, mean_profit_combination,
mean_profit_combination_sub)]
setkey(foo_dt_all, Well_capacity, tot_acres, quarter, CR, PAW)
setkey(foo_irr_2, Well_capacity, tot_acres, quarter, CR, PAW)
lookup_table_quarter = foo_irr_2[foo_dt_all]
lookup_table_quarter = unique(lookup_table_quarter, by=colnames(lookup_table_quarter))
lookup_table_quarter = unique(lookup_table_quarter, by=c("SOIL_ID", "WSTA", "Well_capacity",
"tot_acres", "IFREQ", "CR", "quarter", "PAW", "SDAT"))
lookup_table_well = unique(lookup_table_quarter, by = "Well_capacity")
lookup_table_quarter = lookup_table_quarter[, .(Well_capacity,
SOIL_ID, WSTA, tot_acres, quarter, ifreq, SDAT, CR, PAW,
irrigation_quarter = irrigation, yield_kg_ac, profit_quarter = profit, #subsidy_amount
tax_amount)]
data.table::setkey(lookup_table_well, Well_capacity)
lookup_table_well = lookup_table_well[, .(Well_capacity,
SOIL_ID, WSTA, tot_acres, irr_tot_acres = mean_irrigation_combination,
profit_Well_ID = mean_profit_combination, profit_Well_ID_subsidy = mean_profit_combination_sub, #subsidy_amount
tax_amount)]
lookup_table_quarter = unique(lookup_table_quarter, by=colnames(lookup_table_quarter))
lookup_table_well = unique(lookup_table_well, by=colnames(lookup_table_well))
lookup_table_all_years = copy(lookup_table_quarter)
lookup_table_all_years[, `:=`(irr_tot_acres, sum(irrigation_quarter)),
by = c("Well_capacity", "SDAT")]
lookup_table_all_years[, `:=`(profit_Well_ID, sum(profit_quarter)),
by = c("Well_capacity", "SDAT")]
lookup_table_all_years[, `:=`(irr_below, ifelse(irr_tot_acres <
subsidy_threshold, subsidy_threshold - irr_tot_acres,
0))]
lookup_table_all_years = unique(lookup_table_all_years,
by = c("Well_capacity", "SDAT"))
lookup_table_all_years[, `:=`(profit_Well_ID_sub,
profit_Well_ID + irr_below * tax_amount)]
lookup_table_all_years = lookup_table_all_years[, .(Well_capacity,
SOIL_ID, WSTA, tot_acres, SDAT, irr_tot_acres, profit_Well_ID,
irr_below, profit_Well_ID_sub)]
lookup_table_all_years_2 = rbind(lookup_table_all_years_2,
lookup_table_all_years)
lookup_table_quarter_2 = rbind(lookup_table_quarter_2,
lookup_table_quarter)
lookup_table_well_2 = rbind(lookup_table_well_2, lookup_table_well)
print(i)
}
data.table::setkey(lookup_table_all_years_2, SOIL_ID, Well_capacity,
SDAT)
data.table::setkey(lookup_table_quarter_2, SOIL_ID, Well_capacity,
quarter)
data.table::setkey(lookup_table_well_2, SOIL_ID, Well_capacity)
write.csv(lookup_table_well_2, "lookup_table_well_2.csv")
saveRDS(lookup_table_quarter_2, "lookup_table_quarter_2.rds")
saveRDS(lookup_table_all_years_2, "lookup_table_all_years_2.rds")
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.