R/gen_lookup_crop_insurance.R
In manirouhirad/MODSSAT: What the Package Does (One Line, Title Case)
Defines functions
gen_lookup_crop_insurance
Documented in
gen_lookup_crop_insurance
#' This function generates outputs for different crop insurance parameters.
#' @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 well_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_files 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 insparamfile is the value of insurance parameters that are used to estimate the premium rate.
#' @param parcel_input_APH_file is the input APH file at the parcel-level.
#' @param county_input_APH_file is the input APH file at the county-level.
#' @param parcel_output_APH_folder is the output APH file at the parcel-level.
#' @param county_output_APH_folder is the output APH file at the county-level.
#' @param cover_file is the file that includes different coverage levels and their respective subsidies.
#' @param econ_output_folder is the folder where outputs will be written in.
#' @param hydrology_file_year is the file that communicates daily groundwater use for each well. This file is read by the MODFLOW model.
#' @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 missing_soil_types is the name of the soil type that will be assigned to wells that end up with an NA soil type. This is to avoid dropping wells from the analysis.
#' @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 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 minimum_well_capacity is the minimum well capacity in the model. If well capacity falls below this capacity, it is set to this minimum. Defaults to 100 gallons per minute.
#' @param maximum_well_capacity is the maximum well capacity in the model. If well capacity falls above this capacity, it is set to this maximum. Defaults to 3000 gallons per minute.
#' @param first_year_of_GW is the first year that the GW model exists. This may be different than the first year of simulation. Defaults to 1997.
#' @param last_year_of_GW is the last year that the GW model exists. This may be different than the last year of simulation. Defaults to 2007.
#' @param irrigation_season_days Number of days in an irrigation season. Defaults to 70.
#' @param first_year_of_simulation is the first year that the hydro-economic simulation starts. Defaults to 2000.
#' @param insurance_subsidy_increase is the change in insurance subsidy. This can be positive or negative.
#' @param no_policy_indicator is an indicator that is equal to 0 if the policy is no crop insurance policy. Otherwise it's equal to 1.
#' @param num_clusters is the number of cores for parallelization.
#' @return returns the output table.
#' @examples
#' \dontrun{
#' gen_lookup_tax(subsidy_amount = 2)
#' }
#' @export
gen_lookup_crop_insurance = function(DSSAT_files = "./input_files/DSSAT_files",
well_soil_file = "./input_files/Well_Soil Type.csv",
well_capacity_files = "./Well Capacity",
price_file = "./input_files/crop_prices.csv",
fixed_cost_file = "./input_files/fixed_cost_input.csv",
insparamfile = "./input_files/insparam.csv",
parcel_input_APH_file = "./parcel_APH.csv",
county_input_APH_file = "./county_ref.csv",
parcel_output_APH_folder = "./Econ_output/APH_history/APH_parcel/",
county_output_APH_folder = "./Econ_output/APH_history/APH_county/",
cover_file = "./input_files/coverage_subsidy.csv",
econ_output_folder = "./Econ_output/results_with_insurance/",
hydrology_file_year = "./KS_DSSAT_output.csv",
default_well_capacity_col_name = "Well_Capacity(gpm)",
missing_soil_types = "KS00000007",
pumping_cost = 3,
soil_moisture_targets = c(25, 35, 45, 55, 65, 75),
IFREQ_seq = 2,
IFREQ_interpolate = 0.1,
minimum_well_capacity = 0,
maximum_well_capacity = 1000,
first_year_of_GW = 1997,
last_year_of_GW = 2007,
irrigation_season_days = 70,
first_year_of_simulation = 1999,
insurance_subsidy_increase = 0,
no_policy_indicator = 0,
num_clusters = parallel::detectCores()-2
)
{
#............................................................................#
# load necessary packages #
#............................................................................#
library(data.table)
#............................................................................#
# read input files #
#............................................................................#
# DSSAT
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)
# wells: soil type and well capacity <--------- read every year. Fine which year it is.
filenames = list.files(path = well_capacity_files, pattern = "*.csv",
full.names = TRUE)
ldf <- lapply(filenames, fread, fill = T)
ldf <- mapply(cbind, ldf, file_name = filenames, SIMPLIFY = F)
year_dt = rbindlist(ldf)
foo = data.table(Well_ID = 1, V1 = 1, file_name = paste0(well_capacity_files,
"/", first_year_of_simulation, "_Capacity.csv"))
setnames(foo, old = "V1", new = default_well_capacity_col_name)
year_dt = rbind(year_dt, foo)
year_dt[, `:=`(file_name, substr(file_name, nchar(file_name) -
16, nchar(file_name) - 13))]
year_dt[, `:=`(file_name, as.integer(file_name))]
setkey(year_dt, file_name)
year_dt = year_dt[nrow(year_dt)]
year_2 = year_dt$file_name
year_dt[, file_name := ifelse(file_name <= (last_year_of_GW - 1), file_name + 1,
ifelse(file_name > (last_year_of_GW - 1) & file_name <= (last_year_of_GW - 1 + last_year_of_GW - first_year_of_GW + 1), file_name - (-1 + last_year_of_GW - first_year_of_GW + 1),
ifelse(file_name > (last_year_of_GW - 1 + last_year_of_GW - first_year_of_GW + 1) & file_name <= (last_year_of_GW - 1 + 2 * (last_year_of_GW - first_year_of_GW + 1)), file_name - (-1 + 2 * (last_year_of_GW - first_year_of_GW + 1)),
ifelse(file_name > (last_year_of_GW - 1 + 2 * (last_year_of_GW - first_year_of_GW + 1)) & file_name <= (last_year_of_GW - 1 + 3 * (last_year_of_GW - first_year_of_GW + 1)), file_name - (-1 + 3 * (last_year_of_GW - first_year_of_GW + 1)),
file_name - (-1 + 4 * (last_year_of_GW - first_year_of_GW + 1))))))]
year_dt = year_dt$file_name
print(paste0("year_dt is equal to ", year_dt))
print(paste0("year_2 is equal to ", year_2))
# year_dt = 2000
well_capacity_data = rbindlist(ldf)
well_capacity_data[, `:=`(file_name, substr(file_name, nchar(file_name) -
16, nchar(file_name) - 13))]
well_capacity_data[, `:=`(file_name, as.integer(file_name))]
well_capacity_data = well_capacity_data[file_name == year_2,.(Well_ID, `Well_Capacity(gpm)`)]
soil_type = fread(well_soil_file)
soil_type[, V1 := NULL]
soil_type[, `:=`(Soil_Type, gsub("KSFC00000", "KS0000000", Soil_Type))]
data.table::setkey(soil_type, Well_ID)
data.table::setkey(well_capacity_data, Well_ID)
well_capacity_data = soil_type[well_capacity_data]
data.table::setnames(well_capacity_data, old = default_well_capacity_col_name, "Well_capacity")
well_capacity_data[is.na(Soil_Type), Soil_Type:= missing_soil_types]
well_capacity_data[Well_capacity>1300, Well_capacity := 1300]
well_capacity_data[Well_capacity< 5 , Well_capacity := 5]
# parameters: costs, prices
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[Crop == "MZ", `:=`(f_cost, 500)]
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)
# parameters: insurance
inpa = fread(insparamfile)
inpa[, crop := c("MZ", "dry-MZ", "WH", "dry-WH", "SG", "dry-SG")]
inpa = rbind(inpa, data.table(crop="FA", Irrigation = "N", refrate = 0, exp = 1, fixrate = 0, rdfact = 0, unitfact = 1, refamnt = 0))
# parameters: insurance subsidy rate. This is key here.
cover = fread(cover_file)
#variables: APH for fields and county <------------------------------------------ get back to this one
# what needs to be done here is to read the APH and county ref in and estimate APH ratio
parcel_APH = fread(parcel_input_APH_file)
county_APH = fread(county_input_APH_file)
setkey(parcel_APH, CR)
setkey(county_APH, CR)
parcel_APH = parcel_APH[county_APH]
parcel_APH[, APH_ratio := APH/APH_cntyref]
# APH_ratio<-fread(APH_ratiofile) #read in APH
# APH_ratio = unique(APH_ratio, by="CR")
#............................................................................#
# read input files #
#............................................................................#
# Clean DSSAT output data
KS_DSSAT = KS_DSSAT[, .(SOIL_ID, 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 = colnames(KS_DSSAT)[c(3:9)]
KS_DSSAT[, `:=`((cols_change), lapply(.SD, function(x) as.numeric(as.character(x)))), .SDcols = cols_change]
cols_change = colnames(KS_DSSAT)[c(1, 2)]
KS_DSSAT[, `:=`((cols_change), lapply(.SD, as.character)), .SDcols = cols_change]
# KS_DSSAT = KS_DSSAT[SOIL_ID == unique_soil[i]]
qux = KS_DSSAT[(IFREQ == 16 | IFREQ == 18 | IFREQ ==
20) & CR == "MZ"]
qux[, `:=`(HWAM_6, ifelse(IFREQ == 16, HWAM, 0))]
qux[, `:=`(HWAM_10, ifelse(IFREQ == 20, HWAM, 0))]
qux[, `:=`(HWAM_6, max(HWAM_6)), by = c("SOIL_ID", "CR",
"PAW", "SDAT")]
qux[, `:=`(HWAM_10, max(HWAM_10)), by = c("SOIL_ID",
"CR", "PAW", "SDAT")]
qux[, `:=`(HWAM, ifelse(IFREQ == 18, (HWAM_6 + HWAM_10)/2,
HWAM))]
qux[, `:=`(IRCM_6, ifelse(IFREQ == 16, IRCM, 0))]
qux[, `:=`(IRCM_10, ifelse(IFREQ == 20, IRCM, 0))]
qux[, `:=`(IRCM_6, max(IRCM_6)), by = c("SOIL_ID", "CR",
"PAW", "SDAT")]
qux[, `:=`(IRCM_10, max(IRCM_10)), by = c("SOIL_ID",
"CR", "PAW", "SDAT")]
qux[, `:=`(IRCM, ifelse(IFREQ == 18, (IRCM_6 + IRCM_10)/2,
IRCM))]
qux = qux[IFREQ == 18, .(SOIL_ID, CR, IFREQ, PAW, SDAT,
IRCM, PRCP, PRCM, HWAM)]
KS_DSSAT = KS_DSSAT[IFREQ != 18 | CR != "MZ"]
KS_DSSAT = rbind(KS_DSSAT, qux)
setkey(KS_DSSAT, SOIL_ID, CR, PAW, SDAT, IFREQ)
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_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 = 6)]
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, CR, IFREQ, PAW,
SDAT)
KS_DSSAT = KS_DSSAT[PRCP > -10 & HWAM > -10]
KS_DSSAT = KS_DSSAT[, .(SOIL_ID, CR, PAW, SDAT, IFREQ,
IRCM, PRCP, PRCM, HWAM)]
data.table::setkey(KS_DSSAT, SOIL_ID, CR, PAW, SDAT,
IFREQ)
KS_DSSAT[, `:=`(lead_yield, dplyr::lead(HWAM, n = 1L)),
by = c("SOIL_ID", "CR", "PAW", "SDAT")]
KS_DSSAT[, `:=`(lead_irr_mm, dplyr::lead(IRCM, n = 1L)),
by = c("SOIL_ID", "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", "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, 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, 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[IFREQ == 0, `:=`(CR, paste("dry", CR, sep = "-"))]
# add crop well capacity combinations to existing wells.
number_of_crops = length(KS_DSSAT[, unique(CR)])
setkey(well_capacity_data, Well_ID, Soil_Type)
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[ifreq < 2 & ifreq > 0, `:=`(ifreq,
2)]
well_capacity_data = reshape2::melt(well_capacity_data,
id = c("Well_ID", "Soil_Type", "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,
tot_acres)
well_capacity_data[, `:=`(CR, ifelse(value == 1, "MZ",
ifelse(value == 2, "WH", ifelse(value == 3, "SG",
ifelse(value == 4, "dry-MZ", ifelse(value ==
5, "dry-WH", ifelse(value == 6, "dry-SG", "FA")))))))]
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)
## Add Fallow
KS_DSSAT[, `:=`(IFREQ, round(IFREQ, 1))]
number_of_years = nrow(unique(KS_DSSAT, by = "SDAT"))
foo = data.table(SOIL_ID = data.table(SOIL_ID = c("KS00000001", "KS00000002", "KS00000003", "KS00000004",
"KS00000005", "KS00000006", "KS00000007"),
CR = c("FA", "FA", "FA", "FA", "FA", "FA", "FA"), IFREQ = 0,
PAW = c("25", "25", "25", "25", "25", "25", "25"),
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)))]
foo[, `:=`(PRCM, rep(baz$PRCM, nrow(foo)/nrow(baz)))]
foo = foo[,.(SOIL_ID = SOIL_ID.SOIL_ID, CR = SOIL_ID.CR, IFREQ= SOIL_ID.IFREQ, PAW= SOIL_ID.PAW, SDAT,
irr_mm = SOIL_ID.irr_mm, PRCP, PRCM, yield_kg_ac= SOIL_ID.yield_kg_ac)]
KS_DSSAT = rbind(KS_DSSAT, foo)
data.table::setkey(KS_DSSAT, SOIL_ID, 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, Soil_Type, CR,
ifreq)
# merge DSSAT and well data
foo_irr = merge(well_capacity_data, KS_DSSAT, by.x = c("Soil_Type", "CR", "ifreq"), by.y = c("SOIL_ID", "CR", "IFREQ"),
allow.cartesian = T)
foo_irr = foo_irr[, .(Well_ID, SOIL_ID = Soil_Type, Well_capacity,
tot_acres, IFREQ = ifreq, CR, quarter, PAW, SDAT,
irr_mm, PRCP, PRCM, yield_kg_ac)]
data.table::setkey(foo_irr, Well_capacity)
setkey(well_capacity_data, Well_capacity)
# then merge it with other parameters like prices
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)
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)]
# add different coverage levels
foo_irr[, cover := 0]
foo_irr[, subs := 0]
dt = data.table()
for(j in 1:nrow(cover)){ #loop over coverage and subsidy
fooadd <- copy(foo_irr) #start with the base array
qux = as.numeric(cover[j,1])
fooadd[, cover := rep(qux, nrow(fooadd))] #add coverage level
qux = as.numeric(cover[j,2])
fooadd[, subs := rep(qux, nrow(fooadd))] #add subsidy
dt<-rbind(dt,fooadd) #add onto foo_irr
print(j)
}
foo_irr = rbind(foo_irr, dt)
rm(fooadd) #clean out fooadd
rm(dt) #clean out fooadd
#Generate insurance payout here
setkey(foo_irr,CR)
setkey(inpa, crop)
foo_irr=foo_irr[inpa]
foo_irr=foo_irr[complete.cases(Well_ID)] # remove NA's
foo = unique(parcel_APH, by="Well_ID")
foo[, CR := "FA"]
foo[, APH := 0]
foo[, APH_cntyref := 0]
foo[, APH_ratio := 1]
parcel_APH = rbind(parcel_APH, foo)
setkey(foo_irr, Well_ID, CR)
setkey(parcel_APH, Well_ID, CR)
foo_irr = foo_irr[parcel_APH]
#............................................................................#
# calculate premium rate #
#............................................................................#
foo_irr[, pr := APH_cntyref^exp]
foo_irr[, pr := pr * refrate + fixrate]
foo_irr[, pr := pr * unitfact * rdfact]
foo_irr[CR == "FA", pr := 0]
foo_irr[, prevpr := pr]
foo_irr[, pr := ifelse(pr < prevpr*1.2, pr, prevpr)]
foo_irr[ pr > .999, pr := .999]
foo_irr[, prem := pr * cover * APH * (1-subs)]
foo_irr = foo_irr[, .(Well_ID, SOIL_ID, Well_capacity,
SDAT, PRCP, PRCM, price, cost_per_acre_in, f_cost,
tot_acres, IFREQ, quarter, CR, PAW,
irr_mm, yield_kg_ac,
cover, subs, refrate, exp, fixrate, rdfact, unitfact, refamnt, APH, APH_cntyref, APH_ratio, prem)]
#............................................................................#
# calculate profit #
#............................................................................#
foo_irr[, yield_ins := ifelse(yield_kg_ac > cover * APH * no_policy_indicator, yield_kg_ac, cover * APH)]
foo_irr[, liabpay := ifelse(yield_kg_ac < cover * APH * no_policy_indicator, cover * APH - yield_kg_ac, 0)]
foo_irr[, `:=`(irrigation, 32.5 * irr_mm * 0.0393701)]
foo_irr[, `:=`(profit, 32.5 * (price * yield_ins - f_cost - prem) - irrigation * cost_per_acre_in)]
rm(df_foo)
rm(KS_DSSAT)
rm(KS_DSSAT_0)
rm(well_capacity_data)
foo_irr_2 = foo_irr[SDAT == year_dt]
foo_irr[, row:= 1:.N]
#............................................................................#
# expected profit-max #
#............................................................................#
foo_irr = foo_irr[,.(row, Well_ID, Well_capacity, tot_acres, quarter, CR, PAW, cover, irrigation, yield_kg_ac, profit, liabpay)]
foo_irr[, mean_yield := mean(yield_kg_ac), by=c("Well_ID", "tot_acres", "quarter", "CR", "PAW", "cover")]
foo_irr[, mean_profit_practice := mean(profit), by = c("Well_ID", "tot_acres", "quarter", "CR", "PAW", "cover")]
foo_irr[, mean_irrigation_practice := mean(irrigation), by = c("Well_ID", "tot_acres", "quarter", "CR", "PAW", "cover")]
foo_irr = unique(foo_irr, by = c("Well_ID", "tot_acres", "quarter", "CR", "PAW", "cover"))
foo_irr[, profit_quarter := max(mean_profit_practice), by = c("Well_ID", "tot_acres", "quarter")]
foo_irr = foo_irr[profit_quarter == mean_profit_practice]
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_irrigation_practice)), by = c("Well_ID", "tot_acres")]
foo_irr[, `:=`(profit_Well_ID, max(profit_tot_acres)), by = c("Well_ID")]
setkey(foo_irr, Well_ID, tot_acres)
foo_irr = foo_irr[profit_Well_ID == profit_tot_acres]
# so here you want to see what each producer will do at the beginning of the season
# and then see what year it is and then see the outcome
### add year of the simulation and well capacity here:
### also create a separate output that includes crop yield (APH) and county mean
lookup_table_quarter = foo_irr[, .(Well_ID, tot_acres, quarter, CR, PAW, cover,
exp_irrigation_quarter = mean_irrigation_practice,
exp_yield_quarter = mean_yield,
exp_profit_quarter= profit_quarter,
exp_liabpay = liabpay)]
lookup_table_quarter = unique(lookup_table_quarter)
setkey(lookup_table_quarter, Well_ID, tot_acres, quarter, CR, PAW, cover)
setkey(foo_irr_2, Well_ID, tot_acres, quarter, CR, PAW, cover)
lookup_table_year = foo_irr_2[lookup_table_quarter]
setkey(lookup_table_year, Well_capacity, Well_ID, quarter)
lookup_table_year[, year := year_2]
#............................................................................#
# outputs #
#............................................................................#
# 1. output for the hydrology model:
dt = lookup_table_year[,.(Well_ID, irrigation)]
dt[, irrigation := sum(irrigation), by="Well_ID"]
dt = unique(dt, by="Well_ID")
dt[, `:=`(output_rate_acin_day, irrigation/irrigation_season_days)]
output_hydrology = dt[,.(Well_ID, output_rate_acin_day)]
write.csv(output_hydrology, hydrology_file_year, row.names = FALSE)
# 2. output of the econ model:
# econ_output_in = lookup_table_year[1]
# econ_output_in[, year := 0]
# econ_output_in[, Well_ID := 0]
# write.csv(econ_output_in, "./Econ_output/KS_DSSAT_output_ins.csv", row.names = F)
# write.csv(econ_output_in, "./input_files/KS_DSSAT_output_ins.csv", row.names = F)
econ_output_in = fread("./Econ_output/KS_DSSAT_output_ins.csv") ## <------------- fix the year here
econ_output = rbind(econ_output_in, lookup_table_year)
write.csv(econ_output, paste0(econ_output_folder, "Econ_output_",
insurance_subsidy_increase, "_", year_2+1, ".csv"), row.names = FALSE)
# 3. APH for each well.
parcel_APH_output = lookup_table_year[,.(Well_ID, quarter, q1=1, CR, yield_kg_ac)]
parcel_APH_output[, yield_kg_q := yield_kg_ac * 32.5]
parcel_APH_output[, yield_kg_q := sum(yield_kg_q), by=c("Well_ID", "CR")]
parcel_APH_output[, q1 := sum(q1), by=c("Well_ID", "CR")]
parcel_APH_output = unique(parcel_APH_output, by =c("Well_ID", "CR"))
parcel_APH_output[, yield_kg_ac := yield_kg_q/(q1*32.5)]
parcel_APH_output = parcel_APH_output[,.(Well_ID, year = year_2+1, CR, yield_kg_ac)]
# parcel_APH = parcel_APH[,.(Well_ID, CR, APH)]
#
# setkey(parcel_APH, Well_ID, CR)
# setkey(parcel_APH_output, Well_ID, CR)
#
# parcel_APH_output = parcel_APH_output[parcel_APH]
# parcel_APH_output[is.na(yield_kg_ac), yield_kg_ac := APH]
# parcel_APH_output[, APH := .9*APH + .1*yield_kg_ac]
# parcel_APH_output = parcel_APH_output[,.(Well_ID, year = year_2+1, CR, APH)]
write.csv(parcel_APH_output, paste0(parcel_output_APH_folder, "parcel_APH_", year_2+1, ".csv"), row.names = F)
# 4. APH for the county.
county_APH_output = lookup_table_year[,.(Well_ID, quarter, q1=1, CR, yield_kg_ac)]
county_APH_output[, yield_kg_q := yield_kg_ac * 32.5]
county_APH_output[, yield_kg_q := sum(yield_kg_q), by="CR"]
county_APH_output[, q1 := sum(q1), by="CR"]
county_APH_output = unique(county_APH_output, by ="CR")
county_APH_output[, yield_kg_ac := yield_kg_q/(q1*32.5)]
county_APH_output = county_APH_output[,.(year = year_2+1, CR, yield_kg_ac)] # <~~~~~~~~~~~~~~~~ call it cntyref
# setkey(county_APH, CR)
# setkey(county_APH_output, CR)
# county_APH_output = county_APH_output[county_APH]
# county_APH_output[is.na(yield_kg_ac), yield_kg_ac := APH_cntyref]
# county_APH_output[, APH_cntyref := .9*APH_cntyref + .1*yield_kg_ac]
# county_APH_output = county_APH_output[,.(year = year_2+1, CR, APH_cntyref)]
write.csv(county_APH_output, paste0(county_output_APH_folder, "county_ref_", year_2+1, ".csv"), row.names = F)
}
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_crop_insurance
Documented in gen_lookup_crop_insurance
#' This function generates outputs for different crop insurance parameters.
#' @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 well_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_files 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 insparamfile is the value of insurance parameters that are used to estimate the premium rate.
#' @param parcel_input_APH_file is the input APH file at the parcel-level.
#' @param county_input_APH_file is the input APH file at the county-level.
#' @param parcel_output_APH_folder is the output APH file at the parcel-level.
#' @param county_output_APH_folder is the output APH file at the county-level.
#' @param cover_file is the file that includes different coverage levels and their respective subsidies.
#' @param econ_output_folder is the folder where outputs will be written in.
#' @param hydrology_file_year is the file that communicates daily groundwater use for each well. This file is read by the MODFLOW model.
#' @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 missing_soil_types is the name of the soil type that will be assigned to wells that end up with an NA soil type. This is to avoid dropping wells from the analysis.
#' @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 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 minimum_well_capacity is the minimum well capacity in the model. If well capacity falls below this capacity, it is set to this minimum. Defaults to 100 gallons per minute.
#' @param maximum_well_capacity is the maximum well capacity in the model. If well capacity falls above this capacity, it is set to this maximum. Defaults to 3000 gallons per minute.
#' @param first_year_of_GW is the first year that the GW model exists. This may be different than the first year of simulation. Defaults to 1997.
#' @param last_year_of_GW is the last year that the GW model exists. This may be different than the last year of simulation. Defaults to 2007.
#' @param irrigation_season_days Number of days in an irrigation season. Defaults to 70.
#' @param first_year_of_simulation is the first year that the hydro-economic simulation starts. Defaults to 2000.
#' @param insurance_subsidy_increase is the change in insurance subsidy. This can be positive or negative.
#' @param no_policy_indicator is an indicator that is equal to 0 if the policy is no crop insurance policy. Otherwise it's equal to 1.
#' @param num_clusters is the number of cores for parallelization.
#' @return returns the output table.
#' @examples
#' \dontrun{
#' gen_lookup_tax(subsidy_amount = 2)
#' }
#' @export
gen_lookup_crop_insurance = function(DSSAT_files = "./input_files/DSSAT_files",
well_soil_file = "./input_files/Well_Soil Type.csv",
well_capacity_files = "./Well Capacity",
price_file = "./input_files/crop_prices.csv",
fixed_cost_file = "./input_files/fixed_cost_input.csv",
insparamfile = "./input_files/insparam.csv",
parcel_input_APH_file = "./parcel_APH.csv",
county_input_APH_file = "./county_ref.csv",
parcel_output_APH_folder = "./Econ_output/APH_history/APH_parcel/",
county_output_APH_folder = "./Econ_output/APH_history/APH_county/",
cover_file = "./input_files/coverage_subsidy.csv",
econ_output_folder = "./Econ_output/results_with_insurance/",
hydrology_file_year = "./KS_DSSAT_output.csv",
default_well_capacity_col_name = "Well_Capacity(gpm)",
missing_soil_types = "KS00000007",
pumping_cost = 3,
soil_moisture_targets = c(25, 35, 45, 55, 65, 75),
IFREQ_seq = 2,
IFREQ_interpolate = 0.1,
minimum_well_capacity = 0,
maximum_well_capacity = 1000,
first_year_of_GW = 1997,
last_year_of_GW = 2007,
irrigation_season_days = 70,
first_year_of_simulation = 1999,
insurance_subsidy_increase = 0,
no_policy_indicator = 0,
num_clusters = parallel::detectCores()-2
)
{
#............................................................................#
# load necessary packages #
#............................................................................#
library(data.table)
#............................................................................#
# read input files #
#............................................................................#
# DSSAT
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)
# wells: soil type and well capacity <--------- read every year. Fine which year it is.
filenames = list.files(path = well_capacity_files, pattern = "*.csv",
full.names = TRUE)
ldf <- lapply(filenames, fread, fill = T)
ldf <- mapply(cbind, ldf, file_name = filenames, SIMPLIFY = F)
year_dt = rbindlist(ldf)
foo = data.table(Well_ID = 1, V1 = 1, file_name = paste0(well_capacity_files,
"/", first_year_of_simulation, "_Capacity.csv"))
setnames(foo, old = "V1", new = default_well_capacity_col_name)
year_dt = rbind(year_dt, foo)
year_dt[, `:=`(file_name, substr(file_name, nchar(file_name) -
16, nchar(file_name) - 13))]
year_dt[, `:=`(file_name, as.integer(file_name))]
setkey(year_dt, file_name)
year_dt = year_dt[nrow(year_dt)]
year_2 = year_dt$file_name
year_dt[, file_name := ifelse(file_name <= (last_year_of_GW - 1), file_name + 1,
ifelse(file_name > (last_year_of_GW - 1) & file_name <= (last_year_of_GW - 1 + last_year_of_GW - first_year_of_GW + 1), file_name - (-1 + last_year_of_GW - first_year_of_GW + 1),
ifelse(file_name > (last_year_of_GW - 1 + last_year_of_GW - first_year_of_GW + 1) & file_name <= (last_year_of_GW - 1 + 2 * (last_year_of_GW - first_year_of_GW + 1)), file_name - (-1 + 2 * (last_year_of_GW - first_year_of_GW + 1)),
ifelse(file_name > (last_year_of_GW - 1 + 2 * (last_year_of_GW - first_year_of_GW + 1)) & file_name <= (last_year_of_GW - 1 + 3 * (last_year_of_GW - first_year_of_GW + 1)), file_name - (-1 + 3 * (last_year_of_GW - first_year_of_GW + 1)),
file_name - (-1 + 4 * (last_year_of_GW - first_year_of_GW + 1))))))]
year_dt = year_dt$file_name
print(paste0("year_dt is equal to ", year_dt))
print(paste0("year_2 is equal to ", year_2))
# year_dt = 2000
well_capacity_data = rbindlist(ldf)
well_capacity_data[, `:=`(file_name, substr(file_name, nchar(file_name) -
16, nchar(file_name) - 13))]
well_capacity_data[, `:=`(file_name, as.integer(file_name))]
well_capacity_data = well_capacity_data[file_name == year_2,.(Well_ID, `Well_Capacity(gpm)`)]
soil_type = fread(well_soil_file)
soil_type[, V1 := NULL]
soil_type[, `:=`(Soil_Type, gsub("KSFC00000", "KS0000000", Soil_Type))]
data.table::setkey(soil_type, Well_ID)
data.table::setkey(well_capacity_data, Well_ID)
well_capacity_data = soil_type[well_capacity_data]
data.table::setnames(well_capacity_data, old = default_well_capacity_col_name, "Well_capacity")
well_capacity_data[is.na(Soil_Type), Soil_Type:= missing_soil_types]
well_capacity_data[Well_capacity>1300, Well_capacity := 1300]
well_capacity_data[Well_capacity< 5 , Well_capacity := 5]
# parameters: costs, prices
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[Crop == "MZ", `:=`(f_cost, 500)]
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)
# parameters: insurance
inpa = fread(insparamfile)
inpa[, crop := c("MZ", "dry-MZ", "WH", "dry-WH", "SG", "dry-SG")]
inpa = rbind(inpa, data.table(crop="FA", Irrigation = "N", refrate = 0, exp = 1, fixrate = 0, rdfact = 0, unitfact = 1, refamnt = 0))
# parameters: insurance subsidy rate. This is key here.
cover = fread(cover_file)
#variables: APH for fields and county <------------------------------------------ get back to this one
# what needs to be done here is to read the APH and county ref in and estimate APH ratio
parcel_APH = fread(parcel_input_APH_file)
county_APH = fread(county_input_APH_file)
setkey(parcel_APH, CR)
setkey(county_APH, CR)
parcel_APH = parcel_APH[county_APH]
parcel_APH[, APH_ratio := APH/APH_cntyref]
# APH_ratio<-fread(APH_ratiofile) #read in APH
# APH_ratio = unique(APH_ratio, by="CR")
#............................................................................#
# read input files #
#............................................................................#
# Clean DSSAT output data
KS_DSSAT = KS_DSSAT[, .(SOIL_ID, 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 = colnames(KS_DSSAT)[c(3:9)]
KS_DSSAT[, `:=`((cols_change), lapply(.SD, function(x) as.numeric(as.character(x)))), .SDcols = cols_change]
cols_change = colnames(KS_DSSAT)[c(1, 2)]
KS_DSSAT[, `:=`((cols_change), lapply(.SD, as.character)), .SDcols = cols_change]
# KS_DSSAT = KS_DSSAT[SOIL_ID == unique_soil[i]]
qux = KS_DSSAT[(IFREQ == 16 | IFREQ == 18 | IFREQ ==
20) & CR == "MZ"]
qux[, `:=`(HWAM_6, ifelse(IFREQ == 16, HWAM, 0))]
qux[, `:=`(HWAM_10, ifelse(IFREQ == 20, HWAM, 0))]
qux[, `:=`(HWAM_6, max(HWAM_6)), by = c("SOIL_ID", "CR",
"PAW", "SDAT")]
qux[, `:=`(HWAM_10, max(HWAM_10)), by = c("SOIL_ID",
"CR", "PAW", "SDAT")]
qux[, `:=`(HWAM, ifelse(IFREQ == 18, (HWAM_6 + HWAM_10)/2,
HWAM))]
qux[, `:=`(IRCM_6, ifelse(IFREQ == 16, IRCM, 0))]
qux[, `:=`(IRCM_10, ifelse(IFREQ == 20, IRCM, 0))]
qux[, `:=`(IRCM_6, max(IRCM_6)), by = c("SOIL_ID", "CR",
"PAW", "SDAT")]
qux[, `:=`(IRCM_10, max(IRCM_10)), by = c("SOIL_ID",
"CR", "PAW", "SDAT")]
qux[, `:=`(IRCM, ifelse(IFREQ == 18, (IRCM_6 + IRCM_10)/2,
IRCM))]
qux = qux[IFREQ == 18, .(SOIL_ID, CR, IFREQ, PAW, SDAT,
IRCM, PRCP, PRCM, HWAM)]
KS_DSSAT = KS_DSSAT[IFREQ != 18 | CR != "MZ"]
KS_DSSAT = rbind(KS_DSSAT, qux)
setkey(KS_DSSAT, SOIL_ID, CR, PAW, SDAT, IFREQ)
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_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 = 6)]
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, CR, IFREQ, PAW,
SDAT)
KS_DSSAT = KS_DSSAT[PRCP > -10 & HWAM > -10]
KS_DSSAT = KS_DSSAT[, .(SOIL_ID, CR, PAW, SDAT, IFREQ,
IRCM, PRCP, PRCM, HWAM)]
data.table::setkey(KS_DSSAT, SOIL_ID, CR, PAW, SDAT,
IFREQ)
KS_DSSAT[, `:=`(lead_yield, dplyr::lead(HWAM, n = 1L)),
by = c("SOIL_ID", "CR", "PAW", "SDAT")]
KS_DSSAT[, `:=`(lead_irr_mm, dplyr::lead(IRCM, n = 1L)),
by = c("SOIL_ID", "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", "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, 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, 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[IFREQ == 0, `:=`(CR, paste("dry", CR, sep = "-"))]
# add crop well capacity combinations to existing wells.
number_of_crops = length(KS_DSSAT[, unique(CR)])
setkey(well_capacity_data, Well_ID, Soil_Type)
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[ifreq < 2 & ifreq > 0, `:=`(ifreq,
2)]
well_capacity_data = reshape2::melt(well_capacity_data,
id = c("Well_ID", "Soil_Type", "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,
tot_acres)
well_capacity_data[, `:=`(CR, ifelse(value == 1, "MZ",
ifelse(value == 2, "WH", ifelse(value == 3, "SG",
ifelse(value == 4, "dry-MZ", ifelse(value ==
5, "dry-WH", ifelse(value == 6, "dry-SG", "FA")))))))]
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)
## Add Fallow
KS_DSSAT[, `:=`(IFREQ, round(IFREQ, 1))]
number_of_years = nrow(unique(KS_DSSAT, by = "SDAT"))
foo = data.table(SOIL_ID = data.table(SOIL_ID = c("KS00000001", "KS00000002", "KS00000003", "KS00000004",
"KS00000005", "KS00000006", "KS00000007"),
CR = c("FA", "FA", "FA", "FA", "FA", "FA", "FA"), IFREQ = 0,
PAW = c("25", "25", "25", "25", "25", "25", "25"),
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)))]
foo[, `:=`(PRCM, rep(baz$PRCM, nrow(foo)/nrow(baz)))]
foo = foo[,.(SOIL_ID = SOIL_ID.SOIL_ID, CR = SOIL_ID.CR, IFREQ= SOIL_ID.IFREQ, PAW= SOIL_ID.PAW, SDAT,
irr_mm = SOIL_ID.irr_mm, PRCP, PRCM, yield_kg_ac= SOIL_ID.yield_kg_ac)]
KS_DSSAT = rbind(KS_DSSAT, foo)
data.table::setkey(KS_DSSAT, SOIL_ID, 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, Soil_Type, CR,
ifreq)
# merge DSSAT and well data
foo_irr = merge(well_capacity_data, KS_DSSAT, by.x = c("Soil_Type", "CR", "ifreq"), by.y = c("SOIL_ID", "CR", "IFREQ"),
allow.cartesian = T)
foo_irr = foo_irr[, .(Well_ID, SOIL_ID = Soil_Type, Well_capacity,
tot_acres, IFREQ = ifreq, CR, quarter, PAW, SDAT,
irr_mm, PRCP, PRCM, yield_kg_ac)]
data.table::setkey(foo_irr, Well_capacity)
setkey(well_capacity_data, Well_capacity)
# then merge it with other parameters like prices
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)
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)]
# add different coverage levels
foo_irr[, cover := 0]
foo_irr[, subs := 0]
dt = data.table()
for(j in 1:nrow(cover)){ #loop over coverage and subsidy
fooadd <- copy(foo_irr) #start with the base array
qux = as.numeric(cover[j,1])
fooadd[, cover := rep(qux, nrow(fooadd))] #add coverage level
qux = as.numeric(cover[j,2])
fooadd[, subs := rep(qux, nrow(fooadd))] #add subsidy
dt<-rbind(dt,fooadd) #add onto foo_irr
print(j)
}
foo_irr = rbind(foo_irr, dt)
rm(fooadd) #clean out fooadd
rm(dt) #clean out fooadd
#Generate insurance payout here
setkey(foo_irr,CR)
setkey(inpa, crop)
foo_irr=foo_irr[inpa]
foo_irr=foo_irr[complete.cases(Well_ID)] # remove NA's
foo = unique(parcel_APH, by="Well_ID")
foo[, CR := "FA"]
foo[, APH := 0]
foo[, APH_cntyref := 0]
foo[, APH_ratio := 1]
parcel_APH = rbind(parcel_APH, foo)
setkey(foo_irr, Well_ID, CR)
setkey(parcel_APH, Well_ID, CR)
foo_irr = foo_irr[parcel_APH]
#............................................................................#
# calculate premium rate #
#............................................................................#
foo_irr[, pr := APH_cntyref^exp]
foo_irr[, pr := pr * refrate + fixrate]
foo_irr[, pr := pr * unitfact * rdfact]
foo_irr[CR == "FA", pr := 0]
foo_irr[, prevpr := pr]
foo_irr[, pr := ifelse(pr < prevpr*1.2, pr, prevpr)]
foo_irr[ pr > .999, pr := .999]
foo_irr[, prem := pr * cover * APH * (1-subs)]
foo_irr = foo_irr[, .(Well_ID, SOIL_ID, Well_capacity,
SDAT, PRCP, PRCM, price, cost_per_acre_in, f_cost,
tot_acres, IFREQ, quarter, CR, PAW,
irr_mm, yield_kg_ac,
cover, subs, refrate, exp, fixrate, rdfact, unitfact, refamnt, APH, APH_cntyref, APH_ratio, prem)]
#............................................................................#
# calculate profit #
#............................................................................#
foo_irr[, yield_ins := ifelse(yield_kg_ac > cover * APH * no_policy_indicator, yield_kg_ac, cover * APH)]
foo_irr[, liabpay := ifelse(yield_kg_ac < cover * APH * no_policy_indicator, cover * APH - yield_kg_ac, 0)]
foo_irr[, `:=`(irrigation, 32.5 * irr_mm * 0.0393701)]
foo_irr[, `:=`(profit, 32.5 * (price * yield_ins - f_cost - prem) - irrigation * cost_per_acre_in)]
rm(df_foo)
rm(KS_DSSAT)
rm(KS_DSSAT_0)
rm(well_capacity_data)
foo_irr_2 = foo_irr[SDAT == year_dt]
foo_irr[, row:= 1:.N]
#............................................................................#
# expected profit-max #
#............................................................................#
foo_irr = foo_irr[,.(row, Well_ID, Well_capacity, tot_acres, quarter, CR, PAW, cover, irrigation, yield_kg_ac, profit, liabpay)]
foo_irr[, mean_yield := mean(yield_kg_ac), by=c("Well_ID", "tot_acres", "quarter", "CR", "PAW", "cover")]
foo_irr[, mean_profit_practice := mean(profit), by = c("Well_ID", "tot_acres", "quarter", "CR", "PAW", "cover")]
foo_irr[, mean_irrigation_practice := mean(irrigation), by = c("Well_ID", "tot_acres", "quarter", "CR", "PAW", "cover")]
foo_irr = unique(foo_irr, by = c("Well_ID", "tot_acres", "quarter", "CR", "PAW", "cover"))
foo_irr[, profit_quarter := max(mean_profit_practice), by = c("Well_ID", "tot_acres", "quarter")]
foo_irr = foo_irr[profit_quarter == mean_profit_practice]
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_irrigation_practice)), by = c("Well_ID", "tot_acres")]
foo_irr[, `:=`(profit_Well_ID, max(profit_tot_acres)), by = c("Well_ID")]
setkey(foo_irr, Well_ID, tot_acres)
foo_irr = foo_irr[profit_Well_ID == profit_tot_acres]
# so here you want to see what each producer will do at the beginning of the season
# and then see what year it is and then see the outcome
### add year of the simulation and well capacity here:
### also create a separate output that includes crop yield (APH) and county mean
lookup_table_quarter = foo_irr[, .(Well_ID, tot_acres, quarter, CR, PAW, cover,
exp_irrigation_quarter = mean_irrigation_practice,
exp_yield_quarter = mean_yield,
exp_profit_quarter= profit_quarter,
exp_liabpay = liabpay)]
lookup_table_quarter = unique(lookup_table_quarter)
setkey(lookup_table_quarter, Well_ID, tot_acres, quarter, CR, PAW, cover)
setkey(foo_irr_2, Well_ID, tot_acres, quarter, CR, PAW, cover)
lookup_table_year = foo_irr_2[lookup_table_quarter]
setkey(lookup_table_year, Well_capacity, Well_ID, quarter)
lookup_table_year[, year := year_2]
#............................................................................#
# outputs #
#............................................................................#
# 1. output for the hydrology model:
dt = lookup_table_year[,.(Well_ID, irrigation)]
dt[, irrigation := sum(irrigation), by="Well_ID"]
dt = unique(dt, by="Well_ID")
dt[, `:=`(output_rate_acin_day, irrigation/irrigation_season_days)]
output_hydrology = dt[,.(Well_ID, output_rate_acin_day)]
write.csv(output_hydrology, hydrology_file_year, row.names = FALSE)
# 2. output of the econ model:
# econ_output_in = lookup_table_year[1]
# econ_output_in[, year := 0]
# econ_output_in[, Well_ID := 0]
# write.csv(econ_output_in, "./Econ_output/KS_DSSAT_output_ins.csv", row.names = F)
# write.csv(econ_output_in, "./input_files/KS_DSSAT_output_ins.csv", row.names = F)
econ_output_in = fread("./Econ_output/KS_DSSAT_output_ins.csv") ## <------------- fix the year here
econ_output = rbind(econ_output_in, lookup_table_year)
write.csv(econ_output, paste0(econ_output_folder, "Econ_output_",
insurance_subsidy_increase, "_", year_2+1, ".csv"), row.names = FALSE)
# 3. APH for each well.
parcel_APH_output = lookup_table_year[,.(Well_ID, quarter, q1=1, CR, yield_kg_ac)]
parcel_APH_output[, yield_kg_q := yield_kg_ac * 32.5]
parcel_APH_output[, yield_kg_q := sum(yield_kg_q), by=c("Well_ID", "CR")]
parcel_APH_output[, q1 := sum(q1), by=c("Well_ID", "CR")]
parcel_APH_output = unique(parcel_APH_output, by =c("Well_ID", "CR"))
parcel_APH_output[, yield_kg_ac := yield_kg_q/(q1*32.5)]
parcel_APH_output = parcel_APH_output[,.(Well_ID, year = year_2+1, CR, yield_kg_ac)]
# parcel_APH = parcel_APH[,.(Well_ID, CR, APH)]
#
# setkey(parcel_APH, Well_ID, CR)
# setkey(parcel_APH_output, Well_ID, CR)
#
# parcel_APH_output = parcel_APH_output[parcel_APH]
# parcel_APH_output[is.na(yield_kg_ac), yield_kg_ac := APH]
# parcel_APH_output[, APH := .9*APH + .1*yield_kg_ac]
# parcel_APH_output = parcel_APH_output[,.(Well_ID, year = year_2+1, CR, APH)]
write.csv(parcel_APH_output, paste0(parcel_output_APH_folder, "parcel_APH_", year_2+1, ".csv"), row.names = F)
# 4. APH for the county.
county_APH_output = lookup_table_year[,.(Well_ID, quarter, q1=1, CR, yield_kg_ac)]
county_APH_output[, yield_kg_q := yield_kg_ac * 32.5]
county_APH_output[, yield_kg_q := sum(yield_kg_q), by="CR"]
county_APH_output[, q1 := sum(q1), by="CR"]
county_APH_output = unique(county_APH_output, by ="CR")
county_APH_output[, yield_kg_ac := yield_kg_q/(q1*32.5)]
county_APH_output = county_APH_output[,.(year = year_2+1, CR, yield_kg_ac)] # <~~~~~~~~~~~~~~~~ call it cntyref
# setkey(county_APH, CR)
# setkey(county_APH_output, CR)
# county_APH_output = county_APH_output[county_APH]
# county_APH_output[is.na(yield_kg_ac), yield_kg_ac := APH_cntyref]
# county_APH_output[, APH_cntyref := .9*APH_cntyref + .1*yield_kg_ac]
# county_APH_output = county_APH_output[,.(year = year_2+1, CR, APH_cntyref)]
write.csv(county_APH_output, paste0(county_output_APH_folder, "county_ref_", year_2+1, ".csv"), row.names = F)
}
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.