diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/PEC_soil.R | 3 | ||||
| -rw-r--r-- | R/PEC_sw_drainage_UK.R | 106 | ||||
| -rw-r--r-- | R/PEC_sw_drift.R | 184 | ||||
| -rw-r--r-- | R/PEC_sw_exposit.R | 38 | ||||
| -rw-r--r-- | R/TOXSWA_cwa.R | 4 | ||||
| -rw-r--r-- | R/endpoint.R | 2 | ||||
| -rw-r--r-- | R/twa.R | 25 |
7 files changed, 265 insertions, 97 deletions
diff --git a/R/PEC_soil.R b/R/PEC_soil.R index ac551a7..d6356b9 100644 --- a/R/PEC_soil.R +++ b/R/PEC_soil.R @@ -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 7028101..05f90dd 100644 --- a/R/PEC_sw_drift.R +++ b/R/PEC_sw_drift.R @@ -1,3 +1,4 @@ +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 @@ -7,7 +8,11 @@ #' 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, or with units defined via the @@ -15,19 +20,29 @@ #' @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 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 @@ -53,16 +68,41 @@ #' 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 = 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 = as_units("100 cm"), @@ -70,38 +110,77 @@ PEC_sw_drift <- function(rate, 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) - # Set default units if not specified if (!inherits(rate, "units")) rate <- set_units(rate, rate_units, mode = "symbolic") 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 over waterbody depth else water_width - (water_depth / tanpi(side_angle/180)) 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") - dist_index <- as.character(distances) 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(distances, "m") + } + if (drift_data == "RF") { + drift_percentages <- drift_percentages_rautmann(distances, applications, formula = formula, crop_group_RF, widths = as.numeric(set_units(water_width, "m"))) - ) - names(drift_percentages) <- paste(dist_index, "m") + names(drift_percentages) <- paste(distances, "m") + } } else { names(drift_percentages) <- paste(drift_percentages, "%") } @@ -118,7 +197,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. @@ -131,6 +214,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") @@ -151,41 +244,46 @@ 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 (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) - 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.R b/R/PEC_sw_exposit.R index fe23284..45867c0 100644 --- a/R/PEC_sw_exposit.R +++ b/R/PEC_sw_exposit.R @@ -52,7 +52,8 @@ #' #' @importFrom units as_units set_units drop_units #' @importFrom dplyr across mutate -#' @param rate The 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 intercepted by the crop #' @param Koc The sorption coefficient to soil organic carbon #' @param DT50 The soil half-life in days @@ -69,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 @@ -79,7 +80,8 @@ #' @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 = set_units(Inf, "d"), +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 = set_units(30, "m3"), V_event = set_units(100, "m3"), dilution = 2) @@ -145,6 +147,9 @@ PEC_sw_exposit_runoff <- function(rate, interception = 0, Koc, DT50 = set_units( #' 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 @@ -170,18 +175,31 @@ PEC_sw_exposit_runoff <- function(rate, interception = 0, Koc, DT50 = set_units( #' @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" } } @@ -200,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/TOXSWA_cwa.R b/R/TOXSWA_cwa.R index 310029e..b2f7619 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. 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])) { @@ -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 @@ -108,7 +108,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 +119,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 +131,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 +141,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 +149,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 +159,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 +210,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 +232,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 +246,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 |