diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/GUS.R | 4 | ||||
| -rw-r--r-- | R/PEC_soil.R | 5 | ||||
| -rw-r--r-- | R/PEC_sw_drainage_UK.R | 106 | ||||
| -rw-r--r-- | R/PEC_sw_drift.R | 218 | ||||
| -rw-r--r-- | R/PEC_sw_exposit.R (renamed from R/PEC_sw_exposit_runoff.R) | 88 | ||||
| -rw-r--r-- | R/PEC_sw_sed.R | 22 | ||||
| -rw-r--r-- | R/TOXSWA_cwa.R | 10 | ||||
| -rw-r--r-- | R/endpoint.R | 2 | ||||
| -rw-r--r-- | R/pfm_degradation.R | 35 | ||||
| -rw-r--r-- | R/twa.R | 62 |
10 files changed, 414 insertions, 138 deletions
@@ -1,7 +1,7 @@ #' Groundwater ubiquity score based on Gustafson (1989) #' -#' The groundwater ubiquity score GUS is calculated according to -#' the following equation +#' The groundwater ubiquity score GUS is calculated according to the +#' equation #' \deqn{GUS = \log_{10} DT50_{soil} (4 - \log_{10} K_{oc})}{GUS = log10 DT50soil * (4 - log10 Koc)} #' #' @references Gustafson, David I. (1989) Groundwater ubiquity score: a simple diff --git a/R/PEC_soil.R b/R/PEC_soil.R index ac551a7..b1838d7 100644 --- a/R/PEC_soil.R +++ b/R/PEC_soil.R @@ -19,7 +19,7 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("destination", "study_typ #' Total soil and porewater PEC values for the scenarios as defined in the EFSA #' guidance (2017, p. 14/15) can easily be calculated. #' @note While time weighted average (TWA) concentrations given in the examples -#' from the EFSA guidance from 2015 (p. 80) are be reproduced, this is not +#' from the EFSA guidance from 2015 (p. 80) can be reproduced, this is not #' true for the TWA concentrations given for the same example in the EFSA guidance #' from 2017 (p. 92). #' @note According to the EFSA guidance (EFSA, 2017, p. 43), leaching should be @@ -68,7 +68,8 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("destination", "study_typ #' as Kom here #' @param t_avg Averaging times for time weighted average concentrations #' @param t_act Time series for actual concentrations -#' @param scenarios If this is 'default', the DT50 will be used without correction +#' @param scenarios If this is 'default', a soil bulk density of 1.5 kg/L will +#' be used. The DT50 will be used without correction #' and soil properties as specified in the REACH guidance (R.16, Table #' R.16-9) are used for porewater PEC calculations. If this is "EFSA_2015", #' the DT50 is taken to be a modelling half-life at 20°C and pF2 (for when diff --git a/R/PEC_sw_drainage_UK.R b/R/PEC_sw_drainage_UK.R index d773f40..40835b2 100644 --- a/R/PEC_sw_drainage_UK.R +++ b/R/PEC_sw_drainage_UK.R @@ -1,57 +1,105 @@ #' Calculate initial predicted environmental concentrations in surface water due to drainage using the UK method #' #' This implements the method specified in the UK data requirements handbook and was checked against the spreadsheet -#' published on the CRC website +#' published on the CRC website. Degradation between the end (30 April) and the start (1 October) of +#' the drainage period is taken into account if +#' `latest_application` is specified and the degradation parameters are given either as a `soil_DT50` or a `model`. #' -#' @param rate Application rate in g/ha +#' @param rate Application rate in g/ha or with a compatible unit specified +#' with the units package #' @param interception The fraction of the application rate that does not reach the soil -#' @param Koc The sorption coefficient normalised to organic carbon in L/kg +#' @param Koc The sorption coefficient normalised to organic carbon in L/kg or a unit specified +#' with the units package #' @param latest_application Latest application date, formatted as e.g. "01 July" -#' @param soil_DT50 Soil degradation half-life, if SFO kinetics are to be used +#' @param soil_DT50 Soil degradation half-life, if SFO kinetics are to be used, in +#' days or a time unit specified with the units package #' @param model The soil degradation model to be used. Either one of "FOMC", #' "DFOP", "HS", or "IORE", or an mkinmod object #' @param model_parms A named numeric vector containing the model parameters #' @return The predicted concentration in surface water in µg/L #' @references HSE's Chemicals Regulation Division (CRD) Active substance #' PECsw calculations (for UK specific authorisation requests) -#' \url{https://www.hse.gov.uk/pesticides/topics/pesticide-approvals/pesticides-registration/data-requirements-handbook/fate/active-substance-uk.htm} -#' accessed 2019-09-27 +#' \url{https://www.hse.gov.uk/pesticides/data-requirements-handbook/fate/pecsw-sed-via-drainflow.htm} +#' accessed 2026-02-13 #' -#' Drainage PECs Version 1.0 (2015) Spreadsheet published at -#' \url{https://www.hse.gov.uk/pesticides/topics/pesticide-approvals/pesticides-registration/data-requirements-handbook/fate/pec-tools-2015/PEC\%20sw-sed\%20(drainage).xlsx} -#' accessed 2019-09-27 +#' PECsw/sed spray drift and tier 1 drainflow calculator Version 2.1.1 (2025) Spreadsheet published at +#' \url{https://www.hse.gov.uk/pesticides/assets/docs/PEC%20sw-sed%20(spraydrift).xlsx)} +#' accessed 2026-02-13 #' @export #' @author Johannes Ranke #' @examples #' PEC_sw_drainage_UK(150, Koc = 100) -PEC_sw_drainage_UK <- function(rate, interception = 0, Koc, - latest_application = NULL, soil_DT50 = NULL, - model = NULL, model_parms = NULL) +#' PEC_sw_drainage_UK(60, interception = 0.5, Koc = 550, +#' latest_application = "01 July", soil_DT50 = 200) +PEC_sw_drainage_UK <- function(rate, + interception = 0, Koc, + latest_application = NULL, soil_DT50 = NULL, + model = NULL, model_parms = NULL) { - percentage_lost <- SSLRC_mobility_classification(Koc)[[2]] - amount_available <- rate * (1 - interception) # g/ha + # Set default units if not specified and convert to units used in the calculations + if (!inherits(rate, "units")) rate <- set_units(rate, "g/ha") + rate_g_ha <- as.numeric(set_units(rate, "g/ha")) + + if (!inherits(Koc, "units")) Koc <- set_units(Koc, "L/kg") + Koc_L_kg <- as.numeric(set_units(Koc, "L/kg")) + + if (!missing(soil_DT50)) { + if (!inherits(soil_DT50, "units")) { + soil_DT50_d <- soil_DT50 + } + soil_DT50_d <- as.numeric(set_units(soil_DT50, "d")) + } + + percentage_lost <- SSLRC_mobility_classification(Koc_L_kg)[[2]] + amount_available <- rate_g_ha * (1 - interception) # amount in g for 1 ha if (!missing(latest_application)) { lct <- Sys.getlocale("LC_TIME") tmp <- Sys.setlocale("LC_TIME", "C") - latest <- as.Date(paste(latest_application, "1999"), "%d %b %Y") + if (latest_application == "29 February") { # Use a leap year + ref_year <- 2000 + } else { ref_year <- 1999} # Use a non-leap year + latest <- as.Date(paste(latest_application, ref_year), "%d %b %Y") + if (is.na(latest)) stop("Please specify the latest application in the format '%d %b', e.g. '01 July'") tmp <- Sys.setlocale("LC_TIME", lct) - degradation_time <- as.numeric(difftime(as.Date("1999-10-01"), units = "days", latest)) - if (!missing(soil_DT50)) { - k = log(2)/soil_DT50 - as.Date(paste(latest_application, "1999"), "%d %B %Y") - - amount_available <- amount_available * exp(-k * degradation_time) - if (!missing(model)) stop("You already supplied a soil_DT50 value, implying SFO kinetics") + + drainage_date <- drainage_date_UK(latest) + degradation_time <- as.numeric(difftime(drainage_date, latest, units = "days")) + + if (degradation_time > 0) { + if (!missing(soil_DT50)) { + k = log(2)/soil_DT50_d + amount_available <- amount_available * exp(-k * degradation_time) + if (!missing(model)) stop("You already supplied a soil_DT50 value, implying SFO kinetics") + } + if (!missing(model)) { + fraction_left <- pfm_degradation(model, parms = model_parms, + times = degradation_time)[1, "parent"] + amount_available <- fraction_left * amount_available + } } - if (!missing(model)) { - fraction_left <- pfm_degradation(model, parms = model_parms, - times = degradation_time)[1, "parent"] - amount_available <- fraction_left * amount_available - } - } + } volume = 130000 # L/ha - PEC = 1e6 * (percentage_lost/100) * amount_available / volume + PEC = set_units(1e6 * (percentage_lost/100) * amount_available / volume, "\u00B5g/L") return(PEC) } + +#' @rdname PEC_sw_drainage_UK +#' @param application_date Application date +#' @export +#' @examples +#' drainage_date_UK("2023-07-10") +#' drainage_date_UK("2020-12-01") +#' drainage_date_UK(as.Date("2022-01-15")) +drainage_date_UK <- function(application_date) { + year <- substr(application_date, 1, 4) + drainage_end <- as.Date(paste0(year, "-04-30")) + drainage_start <- as.Date(paste0(year, "-10-01")) + if (application_date <= drainage_end | application_date >= drainage_start) { + drainage_date <- application_date + } else { + drainage_date <- drainage_start + } + return(drainage_date) +} diff --git a/R/PEC_sw_drift.R b/R/PEC_sw_drift.R index c1b24d3..621a0b6 100644 --- a/R/PEC_sw_drift.R +++ b/R/PEC_sw_drift.R @@ -1,27 +1,48 @@ +utils::globalVariables(c("A", "B", "C", "D", "H", "hinge", "z1", "z2", "distance", "pctg", "width")) #' Calculate predicted environmental concentrations in surface water due to drift #' #' This is a basic, vectorised form of a simple calculation of a contaminant #' concentration in surface water based on complete, instantaneous mixing #' with input via spray drift. #' +#' It is recommened to specify the arguments `rate`, `water_depth` and +#' `water_width` using [units::units] from the `units` package. +#' +#' Since pfm version 0.6.5, the function is vectorised with respect to rates, +#' applications, water depth, crop groups and distances +#' #' @inheritParams drift_percentages_rautmann +#' @importFrom testthat capture_output +#' @importFrom units as_units set_units #' @seealso [drift_parameters_focus], [drift_percentages_rautmann] -#' @param rate Application rate in units specified below +#' @param rate Application rate in units specified below, or with units defined via the +#' `units` package. +#' @param rate_units Defaults to g/ha. For backwards compatibility, only used +#' if the specified rate does not have [units::units]]. #' @param drift_percentages Percentage drift values for which to calculate PECsw. -#' Overrides 'drift_data' and 'distances' if not NULL. +#' Overrides 'drift_data', 'distances', 'applications', crop group and +#' formula arguments if not NULL. #' @param drift_data Source of drift percentage data. If 'JKI', the [drift_data_JKI] #' included in the package is used. If 'RF', the Rautmann drift data are calculated #' either in the original form or integrated over the width of the water body, depending #' on the 'formula' argument. #' @param crop_group_JKI When using the 'JKI' drift data, one of the German names -#' as used in [drift_parameters_focus]. Will only be used if drift_data is 'JKI'. +#' as used in [drift_data_JKI]. Will only be used if drift_data is 'JKI'. Available +#' crop groups are "Ackerbau", "Obstbau frueh", "Obstbau spaet", +#' "Weinbau frueh", "Weinbau spaet", "Hopfenbau", "Flaechenkulturen > 900 l/ha" and +#' "Gleisanlagen". #' @param water_depth Depth of the water body in cm -#' @param rate_units Defaults to g/ha #' @param PEC_units Requested units for the calculated PEC. Only µg/L currently supported #' @param water_width Width of the water body in cm #' @param side_angle The angle of the side of the water relative to the bottom which #' is assumed to be horizontal, in degrees. The SYNOPS model assumes 45 degrees here. -#' @return The predicted concentration in surface water +#' @importFrom tibble as_tibble +#' @importFrom dplyr bind_rows +#' @importFrom tidyr pivot_longer +#' @return A numeric vector with the predicted concentration in surface water. +#' In some cases, the vector is named with distances or drift percentages, for +#' backward compatibility with versions before the vectorisation of arguments +#' other than 'distances' was introduced in v0.6.5. #' @export #' @author Johannes Ranke #' @examples @@ -33,6 +54,10 @@ #' # This makes it possible to also use different distances #' PEC_sw_drift(100, distances = c(1, 3, 5, 6, 10, 20, 50, 100), drift_data = "RF") #' +#' # We can also specify distance units explicitly +#' library(units) +#' PEC_sw_drift(100, distances = set_units(c(1, 3, 5, 6, 10, 20, 50, 100), "m"), drift_data = "RF") +#' #' # or consider aerial application #' PEC_sw_drift(100, distances = c(1, 3, 5, 6, 10, 20, 50, 100), drift_data = "RF", #' crop_group_RF = "aerial") @@ -47,50 +72,120 @@ #' PEC_sw_drift(100, drift_data = "RF", formula = "FOCUS") #' PEC_sw_drift(100, drift_data = "RF", formula = "FOCUS", side_angle = 45) #' PEC_sw_drift(100, drift_data = "RF", formula = "FOCUS", side_angle = 45, water_width = 200) +#' +#' # The function is vectorised with respect to rates, applications, water depth, +#' # crop groups and distances +#' PEC_sw_drift( +#' rate = rep(100, 6), +#' applications = c(1, 2, rep(1, 4)), +#' water_depth = c(30, 30, 30, 60, 30, 30), +#' crop_group_JKI = c(rep("Ackerbau", 4), rep("Obstbau frueh", 2)), +#' distances = c(rep(5, 4), 10, 5)) +#' +#' # Try the same with the Rautmann formula +#' PEC_sw_drift( +#' rate = rep(100, 6), +#' applications = c(1, 2, rep(1, 4)), +#' water_depth = c(30, 30, 30, 60, 30, 30), +#' drift_data = "RF", +#' crop_group_RF = c(rep("arable", 4), rep("fruit, early", 2)), +#' distances = c(rep(5, 4), 10, 5)) +#' +#' # And with the FOCUS variant +#' PEC_sw_drift( +#' rate = rep(100, 6), +#' applications = c(1, 2, rep(1, 4)), +#' water_depth = c(30, 30, 30, 60, 30, 30), +#' drift_data = "RF", +#' formula = "FOCUS", +#' crop_group_RF = c(rep("arable", 4), rep("fruit, early", 2)), +#' distances = c(rep(5, 4), 10, 5)) PEC_sw_drift <- function(rate, applications = 1, - water_depth = 30, + water_depth = as_units("30 cm"), drift_percentages = NULL, drift_data = c("JKI", "RF"), - crop_group_JKI = c("Ackerbau", - "Obstbau frueh", "Obstbau spaet", "Weinbau frueh", "Weinbau spaet", - "Hopfenbau", "Flaechenkulturen > 900 l/ha", "Gleisanlagen"), - crop_group_RF = c("arable", "hops", "vines, late", "vines, early", - "fruit, late", "fruit, early", "aerial"), + crop_group_JKI = "Ackerbau", + crop_group_RF = "arable", distances = c(1, 5, 10, 20), formula = c("Rautmann", "FOCUS"), - water_width = 100, + water_width = as_units("100 cm"), side_angle = 90, rate_units = "g/ha", PEC_units = "\u00B5g/L") { + + # Check arguments and set default units if not specified rate_units <- match.arg(rate_units) PEC_units <- match.arg(PEC_units) + if (!inherits(rate, "units")) rate <- set_units(rate, rate_units, mode = "symbolic") + if (!inherits(distances, "units")) distances <- set_units(distances, "m") + if (!inherits(water_width, "units")) water_width <- set_units(water_width, "cm") + if (!inherits(water_depth, "units")) water_depth <- set_units(water_depth, "cm") drift_data <- match.arg(drift_data) - crop_group_JKI <- match.arg(crop_group_JKI) - crop_group_RF <- match.arg(crop_group_RF) - if (drift_data == "JKI" & crop_group_RF != "arable") { + + unmatched_crop_groups_JKI <- setdiff(crop_group_JKI, colnames(pfm::drift_data_JKI[[1]])) + if (length(unmatched_crop_groups_JKI) > 0) { + stop("Crop group(s) ", paste(unmatched_crop_groups_JKI, collapse = ", "), " not supported") + } + + unmatched_crop_groups_RF <- setdiff(crop_group_RF, unique(pfm::drift_parameters_focus$crop_group)) + if (length(unmatched_crop_groups_RF) > 0) { + stop("Crop group(s) ", paste(unmatched_crop_groups_RF, collapse = ", "), "not supported") + } + + if (drift_data == "JKI" & crop_group_RF[1] != "arable") { stop("Specifying crop_group_RF only makes sense if 'RF' is used for 'drift_data'") } - if (drift_data == "RF" & crop_group_JKI != "Ackerbau") { + if (drift_data == "RF" & crop_group_JKI[1] != "Ackerbau") { stop("Specifying crop_group_JKI only makes sense if 'JKI' is used for 'drift_data'") } formula <- match.arg(formula) + + # Check waterbody arguments and calculate mean water width (absolute and relative to water width) if (side_angle < 0 | side_angle > 90) stop("The side anglemust be between 0 and 90 degrees") - mean_water_width <- if (side_angle == 90) water_width + mean_water_width <- if (side_angle == 90) water_width # Mean water width over waterbody depth else water_width - (water_depth / tanpi(side_angle/180)) - water_volume <- 100 * mean_water_width * (water_depth/100) * 1000 # in L (for 1 ha) - PEC_sw_overspray <- rate * 1e6 / water_volume # in µg/L - dist_index <- as.character(distances) + if (as.numeric(mean_water_width) < 0) stop("Undefined geometry") + relative_mean_water_width <- mean_water_width / water_width # Always <= 1 + + # Check lengths of arguments advertised as vectorised for compatibility + arg_lengths <- sapply( + list(rate = rate, applications = applications, distances = distances, + water_depth = water_depth, crop_group_JKI = crop_group_JKI, + crop_group_RF = crop_group_RF), + length) + + arg_lengths_not_one <- arg_lengths[arg_lengths != 1] + if (length(unique(arg_lengths_not_one)) > 1) { + stop("The following argument lengths do not match:\n", + capture_output(print(arg_lengths_not_one))) + } + + # Base PEC sw drift for overspray + PEC_sw_overspray <- set_units(rate / (relative_mean_water_width * water_depth), PEC_units, mode = "symbolic") if (is.null(drift_percentages)) { - drift_percentages <- switch(drift_data, - JKI = pfm::drift_data_JKI[[applications]][dist_index, crop_group_JKI], - RF = drift_percentages_rautmann(distances, applications, + if (drift_data == "JKI") { + drift_data_JKI_long <- pfm::drift_data_JKI |> + lapply(as_tibble, rownames = "distance") |> + bind_rows(.id = "applications") |> + pivot_longer(3:10, names_to = "crop_group_JKI", values_to = "pctg") + + drift_percentages <- tibble( + applications = as.character(applications), + distance = as.character(distances), crop_group_JKI + ) |> + left_join(drift_data_JKI_long, by = c("applications", "distance", "crop_group_JKI")) |> + pull(pctg) + names(drift_percentages) <- paste(as.character(distances), "m") + } + if (drift_data == "RF") { + drift_percentages <- drift_percentages_rautmann(distances, applications, formula = formula, - crop_group_RF, widths = water_width/100) - ) - names(drift_percentages) <- paste(dist_index, "m") + crop_group_RF, widths = as.numeric(set_units(water_width, "m"))) + names(drift_percentages) <- paste(as.character(distances), "m") + } } else { names(drift_percentages) <- paste(drift_percentages, "%") } @@ -107,7 +202,11 @@ PEC_sw_drift <- function(rate, #' @param distances The distances in m for which to get PEC values #' @param widths The widths of the water bodies (only used in the FOCUS formula) #' @param applications Number of applications for selection of drift percentile -#' @param crop_group_RF One of the crop groups as used in [drift_parameters_focus] +#' @param crop_group_RF Crop group(s) as used in [drift_parameters_focus], i.e. +#' "arable", "hops", "vines, late", "vines, early", "fruit, late", "fruit, early" +#' or "aerial". +#' @importFrom tibble tibble +#' @importFrom dplyr if_else left_join mutate pull #' @seealso [drift_parameters_focus], [PEC_sw_drift] #' @references FOCUS (2014) Generic guidance for Surface Water Scenarios (version 1.4). #' FOrum for the Co-ordination of pesticde fate models and their USe. @@ -120,6 +219,16 @@ PEC_sw_drift <- function(rate, #' drift_percentages_rautmann(c(1, 3, 5)) #' drift_percentages_rautmann(c(1, 3, 5), formula = "FOCUS") #' +#' # Since pfm 0.6.5, the function can also take a vector of crop groups +#' drift_percentages_rautmann( +#' distances = c(1, 5, 5), +#' crop_group_RF = c("fruit, early", "fruit, early", "fruit, late")) +#' +#' # Two applications, all else equal +#' drift_data_JKI[[2]][as.character(c(1, 3, 5)), "Ackerbau"] +#' drift_percentages_rautmann(c(1, 3, 5), applications = 2) +#' drift_percentages_rautmann(c(1, 3, 5), formula = "FOCUS", app = 2) +#' #' # One application to early or late fruit crops #' drift_data_JKI[[1]][as.character(c(3, 5, 20, 50)), "Obstbau frueh"] #' drift_percentages_rautmann(c(3, 5, 20, 50), crop_group_RF = "fruit, early") @@ -140,41 +249,50 @@ PEC_sw_drift <- function(rate, #' main = "One application to fruit, early") #' abline(v = 11.4, lty = 2) drift_percentages_rautmann <- function(distances, applications = 1, - crop_group_RF = c("arable", "hops", "vines, late", "vines, early", "fruit, late", - "fruit, early", "aerial"), + crop_group_RF = "arable", formula = c("Rautmann", "FOCUS"), widths = 1 ) { - cg <- match.arg(crop_group_RF) - if (!applications %in% 1:8) stop("Only 1 to 8 applications are supported") + unmatched_crop_groups <- setdiff(crop_group_RF, unique(pfm::drift_parameters_focus$crop_group)) + if (!inherits(distances, "units")) distances <- set_units(distances, "m") + if (!inherits(widths, "units")) widths <- set_units(widths, "m") + if (length(unmatched_crop_groups) > 0) stop("Crop group(s) ", unmatched_crop_groups, " not supported") + if (!all(applications %in% 1:8)) stop("Only 1 to 8 applications are supported") formula <- match.arg(formula) + distances <- drop_units(distances) + widths <- drop_units(widths) - parms <- pfm::drift_parameters_focus[pfm::drift_parameters_focus$crop_group == cg & - pfm::drift_parameters_focus$n_apps == applications, c("A", "B", "C", "D", "hinge")] + # To avoid recycling of components with length != 1 but smaller than the longest argument, + # which would likely be unintended, we use tibble here + parms <- tibble(distance = distances, width = widths, n_apps = applications, crop_group = crop_group_RF) |> + left_join(pfm::drift_parameters_focus, by = c("n_apps", "crop_group")) if (formula[1] == "Rautmann") { - drift_percentages = with(as.list(parms), { - A <- ifelse(distances < hinge, A, C) - B <- ifelse(distances < hinge, B, D) - A * distances^B - }) + drift_percentages <- parms |> + mutate( + A = if_else(distance < hinge, A, C), + B = if_else(distance < hinge, B, D)) |> + mutate( + pctg = A * distances^B) |> + pull(pctg) } else { - drift_percentages = with(as.list(parms), { - z1 = distances - z2 = distances + widths - H = hinge - ifelse(z2 < hinge, + drift_percentages <- parms |> + mutate( + z1 = distance, + z2 = distance + width, + H = hinge) |> + mutate( + pctg = if_else(z2 < hinge, # farther edge closer than hinge distance - A/(widths * (B + 1)) * (z2^(B + 1) - z1^(B + 1)), - ifelse(z1 < hinge, + A/(width * (B + 1)) * (z2^(B + 1) - z1^(B + 1)), + if_else(z1 < hinge, # hinge distance in waterbody (between z1 and z2) - (A/(B + 1) * (H^(B + 1) - z1^(B + 1)) + C/(D + 1) * (z2^(D + 1) - H^(D + 1)))/widths, + (A/(B + 1) * (H^(B + 1) - z1^(B + 1)) + C/(D + 1) * (z2^(D + 1) - H^(D + 1)))/width, # z1 >= hinge, i.e. near edge farther than hinge distance - C/(widths * (D + 1)) * (z2^(D + 1) - z1^(D + 1)) - ) - ) - }) + C/(width * (D + 1)) * (z2^(D + 1) - z1^(D + 1))) + )) |> + pull(pctg) } return(drift_percentages) diff --git a/R/PEC_sw_exposit_runoff.R b/R/PEC_sw_exposit.R index 877dd92..45867c0 100644 --- a/R/PEC_sw_exposit_runoff.R +++ b/R/PEC_sw_exposit.R @@ -47,7 +47,13 @@ #' This is a reimplementation of the calculation described in the Exposit 3.02 spreadsheet file, #' in the worksheet "Konzept Runoff". #' -#' @param rate The application rate in g/ha +#' It is recommened to specify the arguments `rate`, `Koc`, `DT50`, `t_runoff`, `V_ditch` and `V_event` +#' using [units::units] from the `units` package. +#' +#' @importFrom units as_units set_units drop_units +#' @importFrom dplyr across mutate +#' @param rate Application rate in g/ha or with a compatible unit specified +#' with the units package #' @param interception The fraction intercepted by the crop #' @param Koc The sorption coefficient to soil organic carbon #' @param DT50 The soil half-life in days @@ -64,7 +70,7 @@ #' \describe{ #' \item{perc_runoff}{The runoff percentages for dissolved and bound substance} #' \item{runoff}{A matrix containing dissolved and bound input for the different distances} -#' \item{PEC_sw_runoff}{A matrix containing PEC values for dissolved and bound substance +#' \item{PEC_sw_runoff}{A dataframe containing PEC values for dissolved and bound substance #' for the different distances. If the rate was given in g/ha, the PECsw are in microg/L.} #' } #' @export @@ -74,43 +80,61 @@ #' @examples #' PEC_sw_exposit_runoff(500, Koc = 150) #' PEC_sw_exposit_runoff(600, Koc = 10000, DT50 = 195, exposit = "3.01a") -PEC_sw_exposit_runoff <- function(rate, interception = 0, Koc, DT50 = Inf, t_runoff = 3, +PEC_sw_exposit_runoff <- function(rate, interception = 0, Koc, + DT50 = set_units(Inf, "d"), + t_runoff = set_units(3, "days"), exposit_reduction_version = c("3.02", "3.01a", "3.01a2", "2.0"), - V_ditch = 30, V_event = 100, dilution = 2) + V_ditch = set_units(30, "m3"), V_event = set_units(100, "m3"), dilution = 2) { + # Set default units if not specified + if (!inherits(rate, "units")) rate <- set_units(rate, "g/ha") + if (!inherits(Koc, "units")) Koc <- set_units(Koc, "L/kg") + if (!inherits(DT50, "units")) DT50 <- set_units(DT50, "d") + if (!inherits(t_runoff, "units")) t_runoff <- set_units(t_runoff, "d") + if (!inherits(V_ditch, "units")) V_ditch <- set_units(V_ditch, "m3") + if (!inherits(V_event, "units")) V_event <- set_units(V_event, "m3") + k_deg <- log(2)/DT50 - input <- rate * (1 - interception) * 1 * exp(-k_deg * t_runoff) # assumes 1 ha treated area + + # The input is calculated for an area of 1 ha + input <- rate * as_units(1, "ha") * (1 - interception) * exp(as.numeric(-k_deg * t_runoff)) + input_units <- units(input) + input_numeric <- drop_units(input) if (length(Koc) > 1) stop("Only one compound at a time supported") exposit_reduction_version <- match.arg(exposit_reduction_version) - red_water <- pfm::perc_runoff_reduction_exposit[[exposit_reduction_version]]["dissolved"] / 100 - red_bound <- pfm::perc_runoff_reduction_exposit[[exposit_reduction_version]]["bound"] / 100 reduction_runoff <- pfm::perc_runoff_reduction_exposit[[exposit_reduction_version]] / 100 + transfer_runoff <- 1 - reduction_runoff - V_runoff <- V_event * (1 - reduction_runoff[["dissolved"]]) # m3 + V_runoff <- V_event * (1 - reduction_runoff[["dissolved"]]) V_flowing_ditch_runoff <- dilution * (V_ditch + V_runoff) - + f_runoff_exposit <- function(Koc) { - Koc_breaks <- c(pfm::perc_runoff_exposit$Koc_lower_bound, Inf) + Koc_breaks <- c(pfm::perc_runoff_exposit$Koc_lower_bound, set_units(Inf, "L/kg")) Koc_classes <- as.character(cut(Koc, Koc_breaks, labels = rownames(pfm::perc_runoff_exposit))) perc_runoff <- pfm::perc_runoff_exposit[Koc_classes, c("dissolved", "bound")] if (identical(Koc, 0)) perc_runoff <- c(dissolved = 0, bound = 0) return(unlist(perc_runoff) / 100) } f_runoff <- f_runoff_exposit(Koc) - runoff_dissolved <- input * f_runoff["dissolved"] * transfer_runoff["dissolved"] - runoff_bound <- input * f_runoff["bound"] * transfer_runoff["bound"] - runoff_input <- cbind(runoff_dissolved, runoff_bound) - runoff_input$total <- runoff_input$dissolved + runoff_input$bound - PEC_sw_runoff <- 1000 * runoff_input / V_flowing_ditch_runoff + runoff_dissolved <- input_numeric * f_runoff["dissolved"] * transfer_runoff[, "dissolved"] + runoff_bound <- input_numeric * f_runoff["bound"] * transfer_runoff[, "bound"] + runoff_input <- cbind(dissolved = runoff_dissolved, bound = runoff_bound, + total = runoff_dissolved + runoff_bound) + rownames(runoff_input) <- rownames(reduction_runoff) + units(runoff_input) <- input_units + + dn <- dimnames(runoff_input) + PEC_sw_runoff <- set_units(runoff_input / V_flowing_ditch_runoff, "\u00B5g/L") + dimnames(PEC_sw_runoff) <- dn result <- list( perc_runoff = 100 * f_runoff, - runoff = runoff_input, - PEC_sw_runoff = PEC_sw_runoff) + runoff = as.data.frame(runoff_input), + PEC_sw_runoff = as.data.frame(PEC_sw_runoff)) return(result) } @@ -123,6 +147,9 @@ PEC_sw_exposit_runoff <- function(rate, interception = 0, Koc, DT50 = Inf, t_run #' with modest to high mobility (groups 2, 3 and 4). In this implementation, #' the group is derived only from the Koc, if not given explicitly. For #' details, see the discussion of the function arguments below. + +#' It is recommened to specify the arguments `rate`, `Koc`, `DT50`, `t_drainage`, +#' `V_ditch` and `V_drainage` using [units::units] from the `units` package. #' #' @param rate The application rate in g/ha #' @param interception The fraction intercepted by the crop @@ -148,18 +175,31 @@ PEC_sw_exposit_runoff <- function(rate, interception = 0, Koc, DT50 = Inf, t_run #' @seealso \code{\link{perc_runoff_exposit}} for runoff loss percentages and \code{\link{perc_runoff_reduction_exposit}} for runoff reduction percentages used #' @examples #' PEC_sw_exposit_drainage(500, Koc = 150) -PEC_sw_exposit_drainage <- function(rate, interception = 0, Koc = NA, mobility = c(NA, "low", "high"), DT50 = Inf, t_drainage = 3, - V_ditch = 30, V_drainage = c(spring = 10, autumn = 100), dilution = 2) +PEC_sw_exposit_drainage <- function(rate, interception = 0, + Koc = NA, mobility = c(NA, "low", "high"), + DT50 = set_units(Inf, "d"), + t_drainage = set_units(3, "days"), + V_ditch = set_units(30, "m3"), + V_drainage = set_units(c(spring = 10, autumn = 100), "m3"), dilution = 2) { - # Rückstand zum Zeitpunkt des Niederschlagsereignisses (residue at the time of the drainage event) + # Set default units if not specified + if (!inherits(rate, "units")) rate <- set_units(rate, "g/ha") + if (!inherits(Koc, "units")) Koc <- set_units(Koc, "L/kg") + if (!inherits(DT50, "units")) DT50 <- set_units(DT50, "d") + if (!inherits(t_drainage, "units")) t_runoff <- set_units(t_drainage, "d") + if (!inherits(V_ditch, "units")) V_ditch <- set_units(V_ditch, "m3") + if (!inherits(V_drainage, "units")) V_event <- set_units(V_drainage, "m3") + k_deg <- log(2)/DT50 - residue <- rate * (1 - interception) * 1 * exp(-k_deg * t_drainage) # assumes 1 ha treated area + + # Total residue at the time of the drainage event, assumes 1 ha treated area + residue <- rate * as_units(1, "ha") * (1 - interception) * exp(as.numeric(-k_deg * t_drainage)) mobility <- match.arg(mobility) if (is.na(mobility)) { if (is.na(Koc)) stop("Koc is needed if the mobility is not specified") else { - if (Koc > 550) mobility = "low" + if (Koc > set_units(550, "L/kg")) mobility = "low" else mobility = "high" } } @@ -178,11 +218,11 @@ PEC_sw_exposit_drainage <- function(rate, interception = 0, Koc = NA, mobility = f_peak = c(spring = 0.125, autumn = 0.25) # Stoßbelastung (fraction drained at event) - PEC_sw_drainage <- 1000 * residue * f_drainage_total * f_peak / V_flowing_ditch_drainage + PEC_sw_drainage <- residue * f_drainage_total * f_peak / V_flowing_ditch_drainage result <- list( perc_drainage_total = 100 * f_drainage_total, perc_peak = 100 * f_peak, - PEC_sw_drainage = PEC_sw_drainage) + PEC_sw_drainage = set_units(PEC_sw_drainage, "\u00B5g/L")) return(result) } diff --git a/R/PEC_sw_sed.R b/R/PEC_sw_sed.R index 2865ab7..28e12d8 100644 --- a/R/PEC_sw_sed.R +++ b/R/PEC_sw_sed.R @@ -10,24 +10,30 @@ #' @param method The method used for the calculation #' @param sediment_depth Depth of the sediment layer #' @param water_depth Depth of the water body in cm -#' @param sediment_density The density of the sediment in L/kg (equivalent to +#' @param sediment_density The density of the sediment in kg/L (equivalent to #' g/cm3) #' @param PEC_sed_units The units of the estimated sediment PEC value #' @return The predicted concentration in sediment #' @export #' @author Johannes Ranke #' @examples +#' library(pfm) +#' library(units) #' PEC_sw_sed(PEC_sw_drift(100, distances = 1), percentage = 50) PEC_sw_sed <- function(PEC_sw, percentage = 100, method = "percentage", - sediment_depth = 5, water_depth = 30, - sediment_density = 1.3, + sediment_depth = set_units(5, "cm"), + water_depth = set_units(30, "cm"), + sediment_density = set_units(1.3, "kg/L"), PEC_sed_units = c("\u00B5g/kg", "mg/kg")) { - method = match.arg(method) - PEC_sed_units = match.arg(PEC_sed_units) + if (!inherits(PEC_sw, "units")) PEC_sw <- set_units(PEC_sw, "\u00B5g/L") + if (!inherits(sediment_depth, "units")) PEC_sw <- set_units(sediment_depth, "cm") + if (!inherits(water_depth, "units")) PEC_sw <- set_units(water_depth, "cm") + if (!inherits(sediment_density, "units")) PEC_sw <- set_units(sediment_density, "cm") + method <- match.arg(method) + PEC_sed_units <- match.arg(PEC_sed_units) if (method == "percentage") { - PEC_sed = PEC_sw * (percentage/100) * (water_depth / sediment_depth) * (1 / sediment_density) - if (PEC_sed_units == "mg/kg") PEC_sed <- PEC_sed / 1000 + PEC_sed <- PEC_sw * (percentage/100) * as.numeric((water_depth / sediment_depth)) * (1 / sediment_density) } - return(PEC_sed) + return(set_units(PEC_sed, PEC_sed_units, mode = "symbolic")) } diff --git a/R/TOXSWA_cwa.R b/R/TOXSWA_cwa.R index 310029e..ec9b5b7 100644 --- a/R/TOXSWA_cwa.R +++ b/R/TOXSWA_cwa.R @@ -147,10 +147,10 @@ plot.TOXSWA_cwa <- function(x, time_column = c("datetime", "t", "t_firstjan", "t #' and some associated statistics. like maximum moving window average #' concentrations, and dataframes holding the events exceeding specified #' thresholds. Usually, an instance of this class will be generated -#' by \code{\link{read.TOXSWA_cwa}}. +#' by [read.TOXSWA_cwa]. #' #' @export -#' @format An \code{\link{R6Class}} generator object. +#' @format An [R6::R6Class] generator object. #' @field filename Length one character vector holding the filename. #' @field basedir Length one character vector holding the directory where the file came from. #' @field zipfile If not null, giving the path to the zip file from which the file was read. @@ -275,9 +275,9 @@ TOXSWA_cwa <- R6::R6Class("TOXSWA_cwa", cwa_lines <- grep(cwa_string, outfile, value = TRUE) - cwa_all_segments <- read_fwf(paste(cwa_lines, collapse = "\n"), - fwf_empty(paste(tail(cwa_lines), collapse = "\n"))) - + cwa_all_segments <- read_fwf(I(paste(cwa_lines, collapse = "\n")), + fwf_empty(I(paste(tail(cwa_lines), collapse = "\n"))), + show_col_types = FALSE) available_segments = 1:(ncol(cwa_all_segments) - 3) if (segment == "last") segment = max(available_segments) diff --git a/R/endpoint.R b/R/endpoint.R index 08856f9..5415d84 100644 --- a/R/endpoint.R +++ b/R/endpoint.R @@ -41,7 +41,7 @@ endpoint <- function(chent, signif = 3) { if (!is(chent, "chent")) { - stop("Please supply a chent object as created using the package 'chents' available from jrwb.de") + stop("Please supply a chent object as created using the package 'chents' available from github") } ep_list <- chent$chyaml[[medium]][[type]] if (!is.na(lab_field[1])) { diff --git a/R/pfm_degradation.R b/R/pfm_degradation.R index 832a797..63a6e19 100644 --- a/R/pfm_degradation.R +++ b/R/pfm_degradation.R @@ -9,21 +9,48 @@ #' @param years For how many years should the degradation be predicted? #' @param step_days What step size in days should the output have? #' @param times The output times +#' @return A data frame containing the output times and the concentrations +#' assuming initial concentrations of 1 for the parent and zero for +#' metabolites, if any. #' @export #' @author Johannes Ranke #' @examples -#' head(pfm_degradation("SFO", DT50 = 10)) +#' # Simple example of an SFO decline curve +#' sfo_out <- pfm_degradation("SFO", DT50 = 10) +#' head(sfo_out) +#' +#' # Fictive example with a metabolite where we first generate an SFO-SFO model +#' sfo_sfo <- mkinmod( +#' parent = mkinsub("SFO", to = "metabolite"), +#' metabolite = mkinsub("SFO")) +#' +#' sfo_sfo_out <- pfm_degradation(sfo_sfo, +#' parms = c(k_parent = 0.1, f_parent_to_metabolite = 0.5, k_metabolite = 0.02)) +#' +#' plot( +#' sfo_sfo_out[, "time"], +#' sfo_sfo_out[, "parent"], type = "l", +#' xlab = "Time", ylab = "Relative concentration", +#' xlim = c(0, 100)) +#' lines( +#' sfo_sfo_out[, "time"], +#' sfo_sfo_out[, "metabolite"], lty = 2) +#' pfm_degradation <- function(model = "SFO", DT50 = 1000, parms = c(k_parent = log(2)/DT50), years = 1, step_days = 1, times = seq(0, years * 365, by = step_days)) { - if (model %in% c("SFO", "FOMC", "DFOP", "HS", "IORE")) { - model <- mkinmod(parent = list(type = model)) + if (!inherits(model, "mkinmod")) { + if (model[1] %in% c("SFO", "FOMC", "DFOP", "HS", "IORE")) { + model <- mkinmod(parent = list(type = model)) + } else { + stop("Please specify the model with a suitable name or an mkinmod object") + } } initial_state = c(1, rep(0, length(model$diffs) - 1)) names(initial_state) <- names(model$diffs) - time_course <- mkinpredict(model, odeparms = parms, + time_course <- mkinpredict(model, odeparms = parms, odeini = initial_state, outtimes = times, solution_type = ifelse(length(model$spec) == 1, "analytical", "deSolve")) @@ -3,9 +3,9 @@ #' @param x When numeric, this is the half-life to be used for an exponential #' decline. When a character string specifying a parent decline model is given #' e.g. \code{FOMC}, \code{parms} must contain the corresponding parameters. -#' If x is an \code{\link{mkinfit}} object, the decline is calculated from this +#' If x is an [mkinfit] object, the decline is calculated from this #' object. -#' @param ini The initial amount. If x is an \code{\link{mkinfit}} object, and +#' @param ini The initial amount. If x is an [mkinfit] object, and #' ini is 'model', the fitted initial concentrations are used. Otherwise, ini #' must be numeric. If it has length one, it is used for the parent and #' initial values of metabolites are zero, otherwise, it must give values for @@ -13,7 +13,7 @@ #' @param t_end End of the time series #' @param res Resolution of the time series #' @param ... Further arguments passed to methods -#' @return An object of class \code{one_box}, inheriting from \code{\link{ts}}. +#' @return An object of class \code{one_box}, inheriting from [ts]. #' @importFrom stats filter frequency time ts #' @export #' @examples @@ -21,13 +21,22 @@ #' pred_0 <- one_box(10) #' plot(pred_0) #' -#' # Use a fitted mkinfit model +#' # Use a custom mkin model #' require(mkin) +#' SFO_SFO <- mkinmod( +#' parent = mkinsub("SFO", to = "m1"), +#' m1 = mkinsub("SFO")) +#' c_0 = c(parent = 100, m1 = 0) +#' deg_parms = c(k_parent = 0.15, f_parent_to_m1 = 0.5, k_m1 = 0.01) +#' pred_sfo_sfo <- one_box(SFO_SFO, odeparms = deg_parms, ini = c_0) +#' plot(pred_sfo_sfo) +#' +#' # Use a fitted mkin model #' fit <- mkinfit("FOMC", FOCUS_2006_C, quiet = TRUE) #' pred_1 <- one_box(fit) #' plot(pred_1) #' -#' # Use a model with more than one observed variable +#' # Use a fitted model with a metabolite #' m_2 <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO")) #' fit_2 <- mkinfit(m_2, FOCUS_2006_D, quiet = TRUE) #' pred_2 <- one_box(fit_2, ini = "model") @@ -82,8 +91,32 @@ one_box.character <- function(x, ini = 1, parms, ..., } #' @rdname one_box +#' @param odeparms Will be passed to [mkinpredict] +#' @param solution_type Will be passed to [mkinpredict] #' @importFrom mkin mkinpredict #' @export +one_box.mkinmod <- function(x, + ini = c(1, rep(0, length(x$diffs) - 1)), + odeparms, + solution_type = "deSolve", + ..., + t_end = 100, res = 0.01) +{ + + # Set names for initial amounts to be compatible with the model + if (is.null(names(ini))) names(ini) <- names(x$diffs) + + t_out = seq(0, t_end, by = res) + + tmp <- mkinpredict(x, odeparms = odeparms, odeini = ini, + outtimes = t_out, solution_type = solution_type)[, -1, drop = FALSE] + result <- ts(tmp, 0, t_end, frequency = 1/res) + class(result) <- c("one_box", "ts") + return(result) +} + +#' @rdname one_box +#' @export one_box.mkinfit <- function(x, ini = "model", ..., t_end = 100, res = 0.01) { fit <- x if (ini[1] == "model") { @@ -108,7 +141,7 @@ one_box.mkinfit <- function(x, ini = "model", ..., t_end = 100, res = 0.01) { #' Plot time series of decline data #' -#' @param x The object of type \code{\link{one_box}} to be plotted +#' @param x The object of type [one_box] to be plotted #' @param xlim Limits for the x axis #' @param ylim Limits for the y axis #' @param xlab Label for the x axis @@ -119,7 +152,7 @@ one_box.mkinfit <- function(x, ini = "model", ..., t_end = 100, res = 0.01) { #' be shown if max_twa is not NULL. #' @param ... Further arguments passed to methods #' @importFrom stats plot.ts -#' @seealso \code{\link{sawtooth}} +#' @seealso [sawtooth] #' @export #' @examples #' dfop_pred <- one_box("DFOP", parms = c(k1 = 0.2, k2 = 0.02, g = 0.7)) @@ -131,6 +164,7 @@ one_box.mkinfit <- function(x, ini = "model", ..., t_end = 100, res = 0.01) { #' fit_2 <- mkinfit(m_2, FOCUS_2006_D, quiet = TRUE) #' pred_2 <- one_box(fit_2, ini = 1) #' pred_2_saw <- sawtooth(pred_2, 2, 7) +#' plot(pred_2_saw) #' plot(pred_2_saw, max_twa = 21, max_twa_var = "m1") plot.one_box <- function(x, xlim = range(time(x)), ylim = c(0, max(x)), @@ -140,7 +174,7 @@ plot.one_box <- function(x, obs_vars <- dimnames(x)[[2]] plot.ts(x, plot.type = "single", xlab = xlab, ylab = ylab, lty = 1:length(obs_vars), col = 1:length(obs_vars), - las = 1, xlim = xlim, ylim = ylim) + las = 1, xlim = xlim, ylim = ylim, ...) if (!is.null(max_twa)) { x_twa <- max_twa(x, window = max_twa) value <- x_twa$max[max_twa_var] @@ -148,7 +182,9 @@ plot.one_box <- function(x, x_twa$window_end[max_twa_var], value, col = "grey") text(x_twa$window_end[max_twa_var], value, paste("Maximum:", signif(value, 3)), pos = 4) # Plot a second time to cover the grey rectangle - matlines(time(x), as.matrix(x), lty = 1:length(obs_vars), col = 1:length(obs_vars)) + plot.ts(x, plot.type = "single", xlab = xlab, ylab = ylab, + lty = 1:length(obs_vars), col = 1:length(obs_vars), + las = 1, xlim = xlim, ylim = ylim, ...) } } @@ -156,7 +192,7 @@ plot.one_box <- function(x, #' #' If the application pattern is specified in \code{applications}, #' \code{n} and \code{i} are disregarded. -#' @param x A \code{\link{one_box}} object +#' @param x A [one_box] object #' @param n The number of applications. If \code{applications} is specified, \code{n} is ignored #' @param i The interval between applications. If \code{applications} is specified, \code{i} #' is ignored @@ -207,7 +243,7 @@ sawtooth <- function(x, n = 1, i = 365, #' the earliest possible time for the maximum in the time series returned #' is after one window has passed. #' -#' @param x An object of type \code{\link{one_box}} +#' @param x An object of type [one_box] #' @param window The size of the moving window #' @seealso \code{\link{max_twa}} #' @importFrom stats start end @@ -229,7 +265,7 @@ twa.one_box <- function(x, window = 21) resolution = 1/frequency(x) n_filter = window/resolution - result = filter(x, rep(1/n_filter, n_filter), method = "convolution", sides = 1) + result = stats::filter(x, rep(1/n_filter, n_filter), method = "convolution", sides = 1) class(result) = c("one_box", "ts") dimnames(result) <- dimnames(x) return(result) @@ -243,7 +279,7 @@ twa.one_box <- function(x, window = 21) #' \code{\link{plot.one_box}} using the window size for the argument #' \code{max_twa}. #' -#' The method working directly on fitted \code{\link{mkinfit}} objects uses the +#' The method working directly on fitted [mkinfit] objects uses the #' equations given in the PEC soil section of the FOCUS guidance and is restricted #' SFO, FOMC and DFOP models and to the parent compound #' @references FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence and |