From 7233eed00b799e08c31ae971f997b4b3c14eaea2 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Mon, 19 Jun 2017 20:10:21 +0200 Subject: Single line of generated Step12 input file partially validated --- R/PEC_sw_focus.R | 161 +++++++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 122 insertions(+), 39 deletions(-) (limited to 'R') diff --git a/R/PEC_sw_focus.R b/R/PEC_sw_focus.R index e6f2689..77839ab 100644 --- a/R/PEC_sw_focus.R +++ b/R/PEC_sw_focus.R @@ -1,10 +1,12 @@ -#' Calculate FOCUS Step 1 PEC surface water +#' Calculate PEC surface water at FOCUS Step 1 #' -#' This is reimplementation of Step 1 of the FOCUS Step 1 and 2 calculator -#' version 3.2, authored by Michael Klein. Note that results for multiple +#' This is reimplementation of the FOCUS Step 1 and 2 calculator version 3.2, +#' authored by Michael Klein, in R. Note that results for multiple #' applications should be compared to the corresponding results for a #' single application. At current, this is not done automatically in -#' this implementation. +#' this implementation. Only Step 1 PECs are calculated. However, +#' input files are generated that are suitable as input also for Step 2 +#' to be used with the FOCUS calculator. #' #' @importFrom utils read.table #' @references FOCUS (2014) Generic guidance for Surface Water Scenarios (version 1.4). @@ -40,6 +42,12 @@ #' parent or a metabolite #' @param scenario The name of the scenario. Must be one of the scenario #' names given in \code{\link{FOCUS_Step_12_scenarios}} +#' @param region 'n' for Northern Europe or 's' for Southern Europe. If NA, only +#' Step 1 PECsw are calculated +#' @param season 'of' for October to February, 'mm' for March to May, and 'js' +#' for June to September. If NA, only step 1 PECsw are calculated +#' @param interception One of 'no interception' (default), 'minimal crop cover', +#' 'average crop cover' or 'full canopy' #' @param txt_file the name, and potentially the full path to the #' Steps.12 input text file to which the specification of the run(s) #' should be written @@ -48,11 +56,11 @@ #' @param append Should the input text file be appended? #' @examples #' # Parent only -#' dummy_1 <- chent_focus_sw(cwsat = 6000, DT50_ws = 6, Koc = 344.8) +#' dummy_1 <- chent_focus_sw("Dummy 1", cwsat = 6000, DT50_ws = 6, Koc = 344.8) #' PEC_sw_focus(dummy_1, 3000, f_drift = 0) #' #' # Metabolite -#' new_dummy <- chent_focus_sw(mw = 250, Koc = 100) +#' new_dummy <- chent_focus_sw("New Dummy", mw = 250, Koc = 100) #' M1 <- chent_focus_sw(mw = 100, cwsat = 100, DT50_ws = 100, Koc = 50, max_ws = 0, max_soil = 0.5) #' PEC_sw_focus(new_dummy, 1000, scenario = "cereals, winter", met = M1) PEC_sw_focus <- function(parent, rate, n = 1, i = NA, @@ -60,6 +68,9 @@ PEC_sw_focus <- function(parent, rate, n = 1, i = NA, met = NULL, f_drift = NA, f_rd = 0.1, scenario = FOCUS_Step_12_scenarios$names, + region = c("n", "s"), + season = c(NA, "of", "mm", "js"), + interception = "no interception", txt_file = "pesticide.txt", overwrite = FALSE, append = TRUE) { if (n > 1 & is.na(i)) stop("Please specify the interval i if n > 1") @@ -67,17 +78,17 @@ PEC_sw_focus <- function(parent, rate, n = 1, i = NA, if (is.na(f_drift)) { scenario = match.arg(scenario) f_drift = FOCUS_Step_12_scenarios$drift[scenario, "1"] / 100 - # For Step 2 we would/will select the reduced percentiles for multiple apps: - # if (n <= 8) { - # f_drift = FOCUS_Step_12_scenarios$drift[scenario, as.character(n)] / 100 - # } else { - # f_drift = FOCUS_Step_12_scenarios$drift[scenario, ">8"] / 100 - # } + # For Step 2 we would select the reduced percentiles for multiple apps: + if (n <= 8) { + f_drift_2 = FOCUS_Step_12_scenarios$drift[scenario, as.character(n)] / 100 + } else { + f_drift_2 = FOCUS_Step_12_scenarios$drift[scenario, ">8"] / 100 + } } + interception = match.arg(interception) + # Write to txt file if requested - add_line <- function(x) cat(paste0(x, "\r\n"), file = txt, append = TRUE) - add <- function(x) cat(paste(x, "\t"), file = txt, append = TRUE) if (file.exists(txt_file)) { if (append) { txt <- file(txt_file, "a") @@ -89,52 +100,111 @@ PEC_sw_focus <- function(parent, rate, n = 1, i = NA, "appending or overwriting it") } } + } else { + txt <- file(txt_file, "w") } on.exit(close(txt)) + add_line <- function(x) cat(paste0(x, "\r\n"), file = txt, append = TRUE) # Write header to txt file - header <- c("Active Substance", "Compound", "Comment", "Mol mass a.i.", - "Mol mass met.", "Water solubility", "KOC assessed compound", - "KOC parent compound", "DT50", "Max. in Water", + header <- c("Active Substance", "Compound", "Comment", + "Mol mass a.i.", "Mol mass met.", + "Water solubility", + "KOC assessed compound", "KOC parent compound", + "DT50", + "Max. in Water", "Max. in Soil asessed compound", # we reproduce the typo... - "App. Rate", "Number of App.", "Time between app.", - "App. Type", "DT50 soil parent compound", "DT50 soil", - "DT50 water", "DT50 sediment", "Region / Season", + "App. Rate", "Number of App.", "Time between app.", "App. Type", + "DT50 soil parent compound", "DT50 soil", + "DT50 water", "DT50 sediment", + "Region / Season", "Interception class") add_line(paste(header, collapse = "\t")) - if (is.null(met)) { compound = parent$name cwsat = parent$cwsat mw_ratio = 1 max_soil = 1 max_ws = 1 - Koc = parent$Koc + Koc_assessed = parent$Koc DT50_ws = parent$DT50_ws + DT50_soil = parent$DT50_soil } else { compound = met$name cwsat = met$cwsat mw_ratio = met$mw / parent$mw max_soil = met$max_soil max_ws = met$max_ws - Koc = met$Koc + Koc_assessed = met$Koc DT50_ws = met$DT50_ws + DT50_soil = met$DT50_soil } - add(parent$name) - add(compound) - add(comment) - if (is.na(parent$mw)) add("-99.00") - else add(sprintf("%.2f", parent$mw)) - if (is.na(met$mw)) add("-99.00") - else add(sprintf("%.2f", met$mw)) - add(sprintf("%.2f", cwsat)) - add(sprintf("%.2f", Koc)) - if (is.null(met)) add("0.00E+00") - else add(sprintf("%.2f", parent$Koc)) - - # Rates for a single application + # In the Step 12 input file, the application type is the index of the drift + # scenarios, with numbering starting at 0 + app_type = which(scenario == rownames(FOCUS_Step_12_scenarios$drift)) - 1 + # The interception class is a number from 1.00 (no interception) to 4.00 (full canopy) + int_class = which(interception == colnames(FOCUS_Step_12_scenarios$interception)) + + # We calculate Step 2 if a season is specified + # Interception and region/season code + if (is.na(season[1])) { + f_int = 0 + reg_sea = 0 + } else { + f_int = FOCUS_Step_12_scenarios$interception[scenario, interception] + + # the code for region / season is the sum of the region code + reg_code = switch(region, + "n" = 0, + "s" = 3) + # and the season code + sea_code = switch(season, + "of" = 0, + "mm" = 1, + "js" = 2) + reg_sea = reg_code + sea_code + } + + run_txt <- c( + ai = parent$name, compound = compound, + comment = comment) + + run_numeric <- c( + mw_ai = NA, mw_met = NA, + cwsat = cwsat, Koc_assessed = Koc_assessed, Koc_parent = 0, + DT50_ws = DT50_ws, + max_ws = 0, max_soil = 0, + rate = rate, n = n, i = if(is.na(i)) 0 else i, + app_type = app_type, + DT50_soil_parent = 0, DT50_soil = 0, DT50_water = 0, DT50_sediment = 0, + reg_sea = reg_sea, int_class = int_class) + + if (is.null(met)) { + run_numeric["DT50_soil"] = parent$DT50_soil + run_numeric["DT50_water"] = parent$DT50_water + run_numeric["DT50_sediment"] = parent$DT50_sediment + } else { + run_numeric["mw_ai"] = parent$mw + run_numeric["mw_met"] = met$mw + run_numeric["Koc_parent"] = parent$Koc + run_numeric["DT50_soil"] = met$DT50_soil + run_numeric["DT50_soil_parent"] = parent$DT50_soil + run_numeric["DT50_water"] = met$DT50_water + run_numeric["DT50_sediment"] = met$DT50_sediment + } + + print_numeric <- function(x) { + ifelse(is.na(x), " -99.00", + ifelse(x < 0.01, + sprintf("%8.2E", x), + sprintf("%8.2f", x))) + } + run_line <- paste(c(run_txt, print_numeric(run_numeric)), collapse = "\t") + add_line(run_line) + + # Rates for a single application (suffix _s) eq_rate_drift_s = mw_ratio * max_ws * rate # Parent only, or metabolite formed in soil: eq_rate_rd_s = mw_ratio * max_soil * rate @@ -156,6 +226,7 @@ PEC_sw_focus <- function(parent, rate, n = 1, i = NA, input_rd_s = f_rd * eq_rate_rd_s_tot * ratio_field_wb / 10 input_rd = n * input_rd_s + # Step 1 accumulation in the aquatic system # No accumulation between multiple applications if 3 * DT50 < i if (n > 1 && (3 * DT50_ws < i)) { input_drift = input_drift_s @@ -165,13 +236,19 @@ PEC_sw_focus <- function(parent, rate, n = 1, i = NA, input_rd = n * input_rd_s } + # Step 2 accumulation in soil + # Use interception and degradation in soil for Step 2 + k_soil = log(2)/DT50_soil + f_deg_soil = exp(- 4 * k_soil) * # Four days of degradation + (1 - exp(-n * i * k_soil)) / (1 - exp(-i * k_soil)) # multiple applications + # Fraction of compound entering the water phase via runoff/drainage depth_sw = 0.3 # m depth_sed = 0.05 # m depth_sed_eff = 0.01 # m, only relevant for sorption rho_sed = 0.8 # kg/L f_OC = 0.05 # 5% organic carbon in sediment - f_rd_sw = depth_sw / (depth_sw + (depth_sed_eff * rho_sed * f_OC * Koc)) + f_rd_sw = depth_sw / (depth_sw + (depth_sed_eff * rho_sed * f_OC * Koc_assessed)) f_rd_sed = 1 - f_rd_sw # Initial PECs @@ -237,6 +314,9 @@ PEC_sw_focus <- function(parent, rate, n = 1, i = NA, #' @param name Length one character vector containing the name #' @param cwsat Water solubility in mg/L #' @param DT50_ws Half-life in water/sediment systems in days +#' @param DT50_soil Half-life in soil in days +#' @param DT50_water Half-life in water in days (Step 2) +#' @param DT50_sediment Half-life in sediment in days (Step 2) #' @param Koc Partition coefficient between organic carbon and water #' in L/kg. #' @param mw Molar weight in g/mol. @@ -244,8 +324,11 @@ PEC_sw_focus <- function(parent, rate, n = 1, i = NA, #' @param max_ws Maximum observed fraction (dimensionless) in water/sediment #' systems #' @return A list with the substance specific properties -chent_focus_sw <- function(name, Koc, DT50_ws = NA, cwsat = 1000, mw = NA, max_soil = 1, max_ws = 1) +chent_focus_sw <- function(name, Koc, DT50_ws = NA, DT50_soil = NA, + DT50_water = NA, DT50_sediment = NA, + cwsat = 1000, mw = NA, max_soil = 1, max_ws = 1) { - list(Koc = Koc, DT50_ws = DT50_ws, cwsat = cwsat, + list(name = name, Koc = Koc, DT50_ws = DT50_ws, DT50_soil = DT50_soil, + DT50_water = DT50_water, DT50_sediment = DT50_sediment, cwsat = cwsat, mw = mw, max_soil = max_soil, max_ws = max_ws) } -- cgit v1.2.1