R/gen_lookup_tax_CO_2.R
In manirouhirad/MODSSAT: What the Package Does (One Line, Title Case)
Defines functions
gen_lookup_tax_CO_2
Documented in
gen_lookup_tax_CO_2
#' This function generates the lookup table with a tax policy for KS.
#' @param tax_amount is the amount of tax on the unit of groundwater extracted. Defaults to 1.1.
#' @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 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 maximum_well_capacity is the maximum well capacity that needs to be included in the lookup table. Defaults to 1000gpm.
#' @param well_capacity_intervals is the size of interpolation interval for well capacity. Defaults to 20gpm.
#' @param number_of_quarters is the maximum number of quarters that need to be considered in the lookup table generation.
#' @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_CO_2 = function(tax_amount = 1,
DSSAT_files = "./input_files/DSSAT_files",
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,
maximum_well_capacity = 1000,
well_capacity_intervals = 20,
number_of_quarters = 6,
num_clusters = 12
)
{
library(data.table)
tax_amount = (tax_amount - 1)/10
print(paste("this is the tax:", tax_amount, sep = " "))
# FN_optim_tax = function(jj) {
# foo_dt1 = foo_irr[Well_ID_grp == jj]
# foo_dt1[, `:=`(mean_profit_PAW, mean(profit)), by = c("Well_ID", "tot_acres", "quarter", "CR", "PAW")]
# foo_dt1[, `:=`(sd_profit_PAW, sd(profit)), by = c("Well_ID", "tot_acres", "quarter", "CR", "PAW")]
# foo_dt1[, `:=`(mean_irr_PAW, mean(irrigation)), by = c("Well_ID", "tot_acres", "quarter", "CR", "PAW")]
# foo_dt1 = unique(foo_dt1, by = c("Well_ID", "tot_acres", "quarter", "CR", "PAW"))
#
# foo_dt1[, `:=`(profit_quarter, max(mean_profit_PAW)), by = c("Well_ID", "tot_acres", "quarter")]
# foo_dt1 = foo_dt1[profit_quarter == mean_profit_PAW]
# foo_dt1[, `:=`(count, .N), by = c("Well_ID", "tot_acres", "quarter")]
# foo_dt1 = unique(foo_dt1, by = c("Well_ID", "tot_acres", "quarter"))
# foo_dt1 = foo_dt1[Well_capacity > 0 | tot_acres == 0]
#
# foo_dt1[, `:=`(profit_tot_acres, sum(profit_quarter)), by = c("Well_ID", "tot_acres")]
# foo_dt1[, `:=`(irr_tot_acres, sum(mean_irr_PAW)), by = c("Well_ID", "tot_acres")]
# foo_dt1[, `:=`(profit_Well_ID, max(profit_tot_acres)), by = c("Well_ID")]
# foo_dt1 = foo_dt1[profit_Well_ID == profit_tot_acres]
# 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")]
unique_soil = unique(KS_DSSAT[,.(SOIL_ID, WSTA)])
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()
i = 5
for (i in 1:max(KS_DSSAT_2$group)) {
soil_type = data.table(Well_ID = seq(1001, (1000+maximum_well_capacity), by=well_capacity_intervals),
Soil_Type = unique_soil[i, SOIL_ID], weather_station = unique_soil[i, WSTA])
well_capacity = data.table(Well_ID = seq(1001, (1000+maximum_well_capacity), by=well_capacity_intervals),
`Well_Capacity(gpm)` = seq(1, (maximum_well_capacity), by=well_capacity_intervals))
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")
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 + tax_amount)/1)] # this is the one in the updated version. Make sure that the units are right when using this version.
# cost_dt[, `:=`(cost_per_acre_mm, (pumping_cost + tax_amount)/25.4)] # this is the original one. Make sure that the units are right throughout!
fixed_cost = data.table::fread(fixed_cost_file)
# fixed_cost[Crop=="MZ", f_cost := 700]
fixed_cost[irr == 0, `:=`(Crop, paste("dry", Crop, sep = "-"))]
fixed_cost[, `:=`(irr, NULL)]
data.table::setkey(fixed_cost, Crop)
KS_DSSAT = KS_DSSAT_2[group == i, .(SOIL_ID, WSTA, CR,
IFREQ, PAW, HDAT, IRCM, PRCP, PRCM, HWAM)]
KS_DSSAT[, `:=`(IRCM, as.numeric(as.character(IRCM)))]
KS_DSSAT = KS_DSSAT[complete.cases(IRCM)]
KS_DSSAT[, `:=`(HDAT, substr(HDAT, 1, 4))]
cols_change = c("IFREQ", "PAW", "HDAT",
"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[SOIL_ID == unique_soil[i]]
KS_DSSAT = KS_DSSAT[IFREQ <= 25]
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, HDAT)
KS_DSSAT = KS_DSSAT[PRCP > -10 & HWAM > -10]
KS_DSSAT = KS_DSSAT[, .(SOIL_ID, WSTA, CR, PAW, HDAT,
IFREQ, IRCM, PRCP, PRCM, HWAM)]
data.table::setkey(KS_DSSAT, SOIL_ID, WSTA, CR, PAW,
HDAT, IFREQ)
KS_DSSAT[, `:=`(lead_yield, dplyr::lead(HWAM, n = 1L)),
by = c("SOIL_ID", "WSTA", "CR",
"PAW", "HDAT")]
KS_DSSAT[, `:=`(lead_irr_mm, dplyr::lead(IRCM,n = 1L)),
by = c("SOIL_ID", "WSTA", "CR",
"PAW", "HDAT")]
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,
HDAT, irr_mm = irr_int, PRCP, PRCM, yield_kg_ha = yield_int)]
# KS_DSSAT[, `:=`(yield_kg_ac, yield_kg_ac * 0.4046)]
KS_DSSAT[CR %like% "MZ", `:=`(yield_bu_ac, yield_kg_ha * 0.01885)]
KS_DSSAT[CR %like% "SG", `:=`(yield_bu_ac, yield_kg_ha * 0.01831)]
KS_DSSAT[CR %like% "WH", `:=`(yield_bu_ac, yield_kg_ha * 0.01719)]
KS_DSSAT[, `:=`(irr_inch, irr_mm * 0.0393701)]
KS_DSSAT[, yield_kg_ha := NULL]
KS_DSSAT[, irr_mm := NULL]
data.table::setkey(KS_DSSAT, SOIL_ID, WSTA, CR, IFREQ,
PAW, HDAT)
KS_DSSAT = KS_DSSAT[IFREQ == 0 | IFREQ >= IFREQ_seq]
KS_DSSAT = KS_DSSAT[IFREQ != 0 | PAW == soil_moisture_targets[1]]
KS_DSSAT[IFREQ == 0, `:=`(CR, paste("dry",
CR, sep = "-"))]
# KS_DSSAT[CR %like% "MZ", yield_bu_ac := yield_bu_ac * 1.183]
# KS_DSSAT[CR %like% "SG", yield_bu_ac := yield_bu_ac * 1.183]
# KS_DSSAT[CR %like% "WH", yield_bu_ac := yield_bu_ac * 1.156]
number_of_crops = length(KS_DSSAT[, unique(CR)])
crop_names = c(KS_DSSAT[, unique(CR)], "FA")
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)
max_tot_acres = number_of_quarters *32.5
tot_acres_list = c(0:number_of_quarters)*32.5
well_capacity_data <- well_capacity_data[rep(seq_len(nrow(well_capacity_data)),
each = length(tot_acres_list))]
well_capacity_data[, `:=`(tot_acres, rep(tot_acres_list, nrow(well_capacity_data)/length(tot_acres_list)))]
well_capacity_data <- well_capacity_data[rep(seq_len(nrow(well_capacity_data)),
each = (number_of_quarters))]
well_capacity_data[, `:=`(quarter, rep(1:number_of_quarters, nrow(well_capacity_data)/number_of_quarters))]
well_capacity_data <- well_capacity_data[rep(seq_len(nrow(well_capacity_data)), each = (number_of_crops + 1))]
well_capacity_data[, `:=`(CR, rep(crop_names, nrow(well_capacity_data)/(number_of_crops + 1)))]
well_capacity_data[, `:=`(ifreq, ceiling(((tot_acres * 1)/(Well_capacity * 0.053030303149011))*10)/10)]
# well_capacity_data[, `:=`(ifreq, round((tot_acres *
# 1)/(Well_capacity * 0.053030303149011), 1))]
well_capacity_data[ifreq > KS_DSSAT[, max(IFREQ)], ifreq := KS_DSSAT[, max(IFREQ)]]
well_capacity_data[Well_capacity == 0, ifreq := 0]
well_capacity_data[ifreq < 2 & ifreq > 0, `:=`(ifreq, 2)]
# well_capacity_data[ifreq > 10] # <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< HERE: checking whether we can assign ifreq > 10 to ifreq = 10
# 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, quarter, CR, ifreq)
well_capacity_data[, aux := tot_acres/32.5]
well_capacity_data[, CR_aux := ifelse(quarter <= aux, 1, 0)]
well_capacity_data = well_capacity_data[CR_aux == 0 | (CR != "FA" & !(CR %like% "dry"))]
well_capacity_data = well_capacity_data[CR_aux != 0 | (CR == "FA" | (CR %like% "dry"))]
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 = "HDAT"))
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],
HDAT = min(KS_DSSAT$HDAT), irr_inch = 0, PRCP = 200,
PRCM = 200, yield_bu_ac = 0)
foo = foo[rep(seq_len(nrow(foo)), each = number_of_years)]
baz = unique(KS_DSSAT, by = "HDAT")
foo[, `:=`(HDAT, rep(baz$HDAT, 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)
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)
foo_irr = 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 = foo_irr[, .(Well_ID, SOIL_ID = Soil_Type, WSTA = weather_station,
Well_capacity, tot_acres, IFREQ = ifreq, CR, quarter,
PAW, HDAT, irr_inch, PRCP, PRCM, yield_bu_ac)]
data.table::setkey(foo_irr, CR)
foo_irr = foo_irr[price_dt]
foo_irr = foo_irr[complete.cases(Well_ID)]
data.table::setkey(foo_irr, Well_ID)
cost_dt = rbind(data.table::data.table(Well_ID = 1, cost_per_acre_in = unique(cost_dt$cost_per_acre_in)),
cost_dt)
# cost_dt = rbind(data.table::data.table(Well_ID = 1, cost_per_acre_mm = unique(cost_dt$cost_per_acre_mm)),
# cost_dt)
data.table::setkey(cost_dt, Well_ID)
foo_irr = foo_irr[cost_dt]
data.table::setkey(foo_irr, CR)
data.table::setkey(fixed_cost, Crop)
foo_irr = foo_irr[fixed_cost]
foo_irr = foo_irr[complete.cases(Well_ID)]
foo_irr[, `:=`(irrigation, 32.5 * irr_inch)]
foo_irr[, `:=`(profit, 32.5 * (yield_bu_ac * price - f_cost) - irrigation * cost_per_acre_in)]
data.table::setkey(foo_irr, Well_ID, tot_acres, quarter,
HDAT)
foo_irr_2 = data.table::copy(foo_irr)
# foo_irr = data.table::copy(foo_irr_2)
# foo_irr[, `:=`(Well_ID_grp, .GRP), by = "Well_ID"]
#
# aa = max(foo_irr$Well_ID_grp)
#
# 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_optim_tax", #"Well_ID_grp",
# "setnames", "setkey", "tax_amount"), envir = environment())
# foo_dt_all <- parLapply(cl, 1:aa, FN_optim_tax)
# stopCluster(cl)
# }
# else {
# print(Sys.info()[1])
# foo_dt_all <- mclapply(X = 1:aa, FUN = FN_optim_tax,
# mc.cores = num_clusters)
# }
#
# foo_dt_all <- do.call(rbind, foo_dt_all)
# foo_dt_all = data.table(foo_dt_all)
foo_irr[, `:=`(mean_profit_PAW, mean(profit)), by = c("Well_ID", "tot_acres", "quarter", "CR", "PAW")]
# foo_irr[, `:=`(sd_profit_PAW, sd(profit)), by = c("Well_ID", "tot_acres", "quarter", "CR", "PAW")]
foo_irr[, `:=`(mean_irr_PAW, mean(irrigation)), by = c("Well_ID", "tot_acres", "quarter", "CR", "PAW")]
foo_irr = unique(foo_irr, by = c("Well_ID", "tot_acres", "quarter", "CR", "PAW"))
foo_irr[, `:=`(profit_quarter, max(mean_profit_PAW)), by = c("Well_ID", "tot_acres", "quarter")]
foo_irr = foo_irr[profit_quarter == mean_profit_PAW]
foo_irr[, `:=`(count, .N), by = c("Well_ID", "tot_acres", "quarter")]
foo_irr = unique(foo_irr, by = c("Well_ID", "tot_acres", "quarter"))
foo_irr = foo_irr[Well_capacity > 0 | tot_acres == 0]
foo_irr[, `:=`(profit_tot_acres, sum(profit_quarter)), by = c("Well_ID", "tot_acres")]
foo_irr[, `:=`(irr_tot_acres, sum(mean_irr_PAW)), by = c("Well_ID", "tot_acres")]
# foo_irr[, `:=`(profit_Well_ID, max(profit_tot_acres)), by = c("Well_ID")]
# foo_dt_all = foo_irr[profit_Well_ID == profit_tot_acres]
foo_irr[, tot_quarters := tot_acres/32.5]
# lookup_table_quarter = foo_irr[, .(Well_capacity, SOIL_ID, WSTA, tot_quarters, tot_acres, quarter, CR, PAW, mean_irr_PAW, profit_quarter)]
# data.table::setkey(lookup_table_quarter, Well_capacity, tot_quarters, quarter)
# lookup_table_well = unique(foo_irr, by = c("Well_ID", "tot_acres"))
# data.table::setkey(lookup_table_well, Well_capacity)
# lookup_table_well = lookup_table_well[, .(Well_capacity, SOIL_ID, WSTA, tot_acres, irr_tot_acres, profit_Well_ID)]
# # lookup_table_well = lookup_table_well[Well_capacity > 0, .(Well_capacity, SOIL_ID, WSTA, tot_acres, irr_tot_acres, profit_Well_ID)]
foo_irr = foo_irr[, .(Well_capacity, tot_quarters, tot_acres,
quarter, CR, PAW, id = 1)]
# baz = lookup_table_quarter[, .(Well_capacity, tot_acres,
# quarter, CR, PAW, id = 1)]
data.table::setkey(foo_irr, Well_capacity, tot_acres, quarter,
CR, PAW)
data.table::setkey(foo_irr_2, Well_capacity, tot_acres,
quarter, CR, PAW)
foo_irr = foo_irr_2[foo_irr]
lookup_table_all_years = foo_irr[, .(Well_ID, Well_capacity,
SOIL_ID, WSTA, tot_quarters, tot_acres, IFREQ, CR, quarter, PAW, HDAT,
PRCP, PRCM, yield_bu_ac, irrigation, profit)]
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,
HDAT)
# 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_CO_2
Documented in gen_lookup_tax_CO_2
#' This function generates the lookup table with a tax policy for KS.
#' @param tax_amount is the amount of tax on the unit of groundwater extracted. Defaults to 1.1.
#' @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 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 maximum_well_capacity is the maximum well capacity that needs to be included in the lookup table. Defaults to 1000gpm.
#' @param well_capacity_intervals is the size of interpolation interval for well capacity. Defaults to 20gpm.
#' @param number_of_quarters is the maximum number of quarters that need to be considered in the lookup table generation.
#' @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_CO_2 = function(tax_amount = 1,
DSSAT_files = "./input_files/DSSAT_files",
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,
maximum_well_capacity = 1000,
well_capacity_intervals = 20,
number_of_quarters = 6,
num_clusters = 12
)
{
library(data.table)
tax_amount = (tax_amount - 1)/10
print(paste("this is the tax:", tax_amount, sep = " "))
# FN_optim_tax = function(jj) {
# foo_dt1 = foo_irr[Well_ID_grp == jj]
# foo_dt1[, `:=`(mean_profit_PAW, mean(profit)), by = c("Well_ID", "tot_acres", "quarter", "CR", "PAW")]
# foo_dt1[, `:=`(sd_profit_PAW, sd(profit)), by = c("Well_ID", "tot_acres", "quarter", "CR", "PAW")]
# foo_dt1[, `:=`(mean_irr_PAW, mean(irrigation)), by = c("Well_ID", "tot_acres", "quarter", "CR", "PAW")]
# foo_dt1 = unique(foo_dt1, by = c("Well_ID", "tot_acres", "quarter", "CR", "PAW"))
#
# foo_dt1[, `:=`(profit_quarter, max(mean_profit_PAW)), by = c("Well_ID", "tot_acres", "quarter")]
# foo_dt1 = foo_dt1[profit_quarter == mean_profit_PAW]
# foo_dt1[, `:=`(count, .N), by = c("Well_ID", "tot_acres", "quarter")]
# foo_dt1 = unique(foo_dt1, by = c("Well_ID", "tot_acres", "quarter"))
# foo_dt1 = foo_dt1[Well_capacity > 0 | tot_acres == 0]
#
# foo_dt1[, `:=`(profit_tot_acres, sum(profit_quarter)), by = c("Well_ID", "tot_acres")]
# foo_dt1[, `:=`(irr_tot_acres, sum(mean_irr_PAW)), by = c("Well_ID", "tot_acres")]
# foo_dt1[, `:=`(profit_Well_ID, max(profit_tot_acres)), by = c("Well_ID")]
# foo_dt1 = foo_dt1[profit_Well_ID == profit_tot_acres]
# 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")]
unique_soil = unique(KS_DSSAT[,.(SOIL_ID, WSTA)])
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()
i = 5
for (i in 1:max(KS_DSSAT_2$group)) {
soil_type = data.table(Well_ID = seq(1001, (1000+maximum_well_capacity), by=well_capacity_intervals),
Soil_Type = unique_soil[i, SOIL_ID], weather_station = unique_soil[i, WSTA])
well_capacity = data.table(Well_ID = seq(1001, (1000+maximum_well_capacity), by=well_capacity_intervals),
`Well_Capacity(gpm)` = seq(1, (maximum_well_capacity), by=well_capacity_intervals))
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")
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 + tax_amount)/1)] # this is the one in the updated version. Make sure that the units are right when using this version.
# cost_dt[, `:=`(cost_per_acre_mm, (pumping_cost + tax_amount)/25.4)] # this is the original one. Make sure that the units are right throughout!
fixed_cost = data.table::fread(fixed_cost_file)
# fixed_cost[Crop=="MZ", f_cost := 700]
fixed_cost[irr == 0, `:=`(Crop, paste("dry", Crop, sep = "-"))]
fixed_cost[, `:=`(irr, NULL)]
data.table::setkey(fixed_cost, Crop)
KS_DSSAT = KS_DSSAT_2[group == i, .(SOIL_ID, WSTA, CR,
IFREQ, PAW, HDAT, IRCM, PRCP, PRCM, HWAM)]
KS_DSSAT[, `:=`(IRCM, as.numeric(as.character(IRCM)))]
KS_DSSAT = KS_DSSAT[complete.cases(IRCM)]
KS_DSSAT[, `:=`(HDAT, substr(HDAT, 1, 4))]
cols_change = c("IFREQ", "PAW", "HDAT",
"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[SOIL_ID == unique_soil[i]]
KS_DSSAT = KS_DSSAT[IFREQ <= 25]
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, HDAT)
KS_DSSAT = KS_DSSAT[PRCP > -10 & HWAM > -10]
KS_DSSAT = KS_DSSAT[, .(SOIL_ID, WSTA, CR, PAW, HDAT,
IFREQ, IRCM, PRCP, PRCM, HWAM)]
data.table::setkey(KS_DSSAT, SOIL_ID, WSTA, CR, PAW,
HDAT, IFREQ)
KS_DSSAT[, `:=`(lead_yield, dplyr::lead(HWAM, n = 1L)),
by = c("SOIL_ID", "WSTA", "CR",
"PAW", "HDAT")]
KS_DSSAT[, `:=`(lead_irr_mm, dplyr::lead(IRCM,n = 1L)),
by = c("SOIL_ID", "WSTA", "CR",
"PAW", "HDAT")]
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,
HDAT, irr_mm = irr_int, PRCP, PRCM, yield_kg_ha = yield_int)]
# KS_DSSAT[, `:=`(yield_kg_ac, yield_kg_ac * 0.4046)]
KS_DSSAT[CR %like% "MZ", `:=`(yield_bu_ac, yield_kg_ha * 0.01885)]
KS_DSSAT[CR %like% "SG", `:=`(yield_bu_ac, yield_kg_ha * 0.01831)]
KS_DSSAT[CR %like% "WH", `:=`(yield_bu_ac, yield_kg_ha * 0.01719)]
KS_DSSAT[, `:=`(irr_inch, irr_mm * 0.0393701)]
KS_DSSAT[, yield_kg_ha := NULL]
KS_DSSAT[, irr_mm := NULL]
data.table::setkey(KS_DSSAT, SOIL_ID, WSTA, CR, IFREQ,
PAW, HDAT)
KS_DSSAT = KS_DSSAT[IFREQ == 0 | IFREQ >= IFREQ_seq]
KS_DSSAT = KS_DSSAT[IFREQ != 0 | PAW == soil_moisture_targets[1]]
KS_DSSAT[IFREQ == 0, `:=`(CR, paste("dry",
CR, sep = "-"))]
# KS_DSSAT[CR %like% "MZ", yield_bu_ac := yield_bu_ac * 1.183]
# KS_DSSAT[CR %like% "SG", yield_bu_ac := yield_bu_ac * 1.183]
# KS_DSSAT[CR %like% "WH", yield_bu_ac := yield_bu_ac * 1.156]
number_of_crops = length(KS_DSSAT[, unique(CR)])
crop_names = c(KS_DSSAT[, unique(CR)], "FA")
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)
max_tot_acres = number_of_quarters *32.5
tot_acres_list = c(0:number_of_quarters)*32.5
well_capacity_data <- well_capacity_data[rep(seq_len(nrow(well_capacity_data)),
each = length(tot_acres_list))]
well_capacity_data[, `:=`(tot_acres, rep(tot_acres_list, nrow(well_capacity_data)/length(tot_acres_list)))]
well_capacity_data <- well_capacity_data[rep(seq_len(nrow(well_capacity_data)),
each = (number_of_quarters))]
well_capacity_data[, `:=`(quarter, rep(1:number_of_quarters, nrow(well_capacity_data)/number_of_quarters))]
well_capacity_data <- well_capacity_data[rep(seq_len(nrow(well_capacity_data)), each = (number_of_crops + 1))]
well_capacity_data[, `:=`(CR, rep(crop_names, nrow(well_capacity_data)/(number_of_crops + 1)))]
well_capacity_data[, `:=`(ifreq, ceiling(((tot_acres * 1)/(Well_capacity * 0.053030303149011))*10)/10)]
# well_capacity_data[, `:=`(ifreq, round((tot_acres *
# 1)/(Well_capacity * 0.053030303149011), 1))]
well_capacity_data[ifreq > KS_DSSAT[, max(IFREQ)], ifreq := KS_DSSAT[, max(IFREQ)]]
well_capacity_data[Well_capacity == 0, ifreq := 0]
well_capacity_data[ifreq < 2 & ifreq > 0, `:=`(ifreq, 2)]
# well_capacity_data[ifreq > 10] # <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< HERE: checking whether we can assign ifreq > 10 to ifreq = 10
# 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, quarter, CR, ifreq)
well_capacity_data[, aux := tot_acres/32.5]
well_capacity_data[, CR_aux := ifelse(quarter <= aux, 1, 0)]
well_capacity_data = well_capacity_data[CR_aux == 0 | (CR != "FA" & !(CR %like% "dry"))]
well_capacity_data = well_capacity_data[CR_aux != 0 | (CR == "FA" | (CR %like% "dry"))]
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 = "HDAT"))
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],
HDAT = min(KS_DSSAT$HDAT), irr_inch = 0, PRCP = 200,
PRCM = 200, yield_bu_ac = 0)
foo = foo[rep(seq_len(nrow(foo)), each = number_of_years)]
baz = unique(KS_DSSAT, by = "HDAT")
foo[, `:=`(HDAT, rep(baz$HDAT, 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)
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)
foo_irr = 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 = foo_irr[, .(Well_ID, SOIL_ID = Soil_Type, WSTA = weather_station,
Well_capacity, tot_acres, IFREQ = ifreq, CR, quarter,
PAW, HDAT, irr_inch, PRCP, PRCM, yield_bu_ac)]
data.table::setkey(foo_irr, CR)
foo_irr = foo_irr[price_dt]
foo_irr = foo_irr[complete.cases(Well_ID)]
data.table::setkey(foo_irr, Well_ID)
cost_dt = rbind(data.table::data.table(Well_ID = 1, cost_per_acre_in = unique(cost_dt$cost_per_acre_in)),
cost_dt)
# cost_dt = rbind(data.table::data.table(Well_ID = 1, cost_per_acre_mm = unique(cost_dt$cost_per_acre_mm)),
# cost_dt)
data.table::setkey(cost_dt, Well_ID)
foo_irr = foo_irr[cost_dt]
data.table::setkey(foo_irr, CR)
data.table::setkey(fixed_cost, Crop)
foo_irr = foo_irr[fixed_cost]
foo_irr = foo_irr[complete.cases(Well_ID)]
foo_irr[, `:=`(irrigation, 32.5 * irr_inch)]
foo_irr[, `:=`(profit, 32.5 * (yield_bu_ac * price - f_cost) - irrigation * cost_per_acre_in)]
data.table::setkey(foo_irr, Well_ID, tot_acres, quarter,
HDAT)
foo_irr_2 = data.table::copy(foo_irr)
# foo_irr = data.table::copy(foo_irr_2)
# foo_irr[, `:=`(Well_ID_grp, .GRP), by = "Well_ID"]
#
# aa = max(foo_irr$Well_ID_grp)
#
# 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_optim_tax", #"Well_ID_grp",
# "setnames", "setkey", "tax_amount"), envir = environment())
# foo_dt_all <- parLapply(cl, 1:aa, FN_optim_tax)
# stopCluster(cl)
# }
# else {
# print(Sys.info()[1])
# foo_dt_all <- mclapply(X = 1:aa, FUN = FN_optim_tax,
# mc.cores = num_clusters)
# }
#
# foo_dt_all <- do.call(rbind, foo_dt_all)
# foo_dt_all = data.table(foo_dt_all)
foo_irr[, `:=`(mean_profit_PAW, mean(profit)), by = c("Well_ID", "tot_acres", "quarter", "CR", "PAW")]
# foo_irr[, `:=`(sd_profit_PAW, sd(profit)), by = c("Well_ID", "tot_acres", "quarter", "CR", "PAW")]
foo_irr[, `:=`(mean_irr_PAW, mean(irrigation)), by = c("Well_ID", "tot_acres", "quarter", "CR", "PAW")]
foo_irr = unique(foo_irr, by = c("Well_ID", "tot_acres", "quarter", "CR", "PAW"))
foo_irr[, `:=`(profit_quarter, max(mean_profit_PAW)), by = c("Well_ID", "tot_acres", "quarter")]
foo_irr = foo_irr[profit_quarter == mean_profit_PAW]
foo_irr[, `:=`(count, .N), by = c("Well_ID", "tot_acres", "quarter")]
foo_irr = unique(foo_irr, by = c("Well_ID", "tot_acres", "quarter"))
foo_irr = foo_irr[Well_capacity > 0 | tot_acres == 0]
foo_irr[, `:=`(profit_tot_acres, sum(profit_quarter)), by = c("Well_ID", "tot_acres")]
foo_irr[, `:=`(irr_tot_acres, sum(mean_irr_PAW)), by = c("Well_ID", "tot_acres")]
# foo_irr[, `:=`(profit_Well_ID, max(profit_tot_acres)), by = c("Well_ID")]
# foo_dt_all = foo_irr[profit_Well_ID == profit_tot_acres]
foo_irr[, tot_quarters := tot_acres/32.5]
# lookup_table_quarter = foo_irr[, .(Well_capacity, SOIL_ID, WSTA, tot_quarters, tot_acres, quarter, CR, PAW, mean_irr_PAW, profit_quarter)]
# data.table::setkey(lookup_table_quarter, Well_capacity, tot_quarters, quarter)
# lookup_table_well = unique(foo_irr, by = c("Well_ID", "tot_acres"))
# data.table::setkey(lookup_table_well, Well_capacity)
# lookup_table_well = lookup_table_well[, .(Well_capacity, SOIL_ID, WSTA, tot_acres, irr_tot_acres, profit_Well_ID)]
# # lookup_table_well = lookup_table_well[Well_capacity > 0, .(Well_capacity, SOIL_ID, WSTA, tot_acres, irr_tot_acres, profit_Well_ID)]
foo_irr = foo_irr[, .(Well_capacity, tot_quarters, tot_acres,
quarter, CR, PAW, id = 1)]
# baz = lookup_table_quarter[, .(Well_capacity, tot_acres,
# quarter, CR, PAW, id = 1)]
data.table::setkey(foo_irr, Well_capacity, tot_acres, quarter,
CR, PAW)
data.table::setkey(foo_irr_2, Well_capacity, tot_acres,
quarter, CR, PAW)
foo_irr = foo_irr_2[foo_irr]
lookup_table_all_years = foo_irr[, .(Well_ID, Well_capacity,
SOIL_ID, WSTA, tot_quarters, tot_acres, IFREQ, CR, quarter, PAW, HDAT,
PRCP, PRCM, yield_bu_ac, irrigation, profit)]
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,
HDAT)
# 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.