diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/PEC_soil.R | 3 | ||||
| -rw-r--r-- | R/PEC_sw_drainage_UK.R | 36 | ||||
| -rw-r--r-- | R/PEC_sw_drift.R | 2 | ||||
| -rw-r--r-- | R/PEC_sw_exposit.R | 33 | ||||
| -rw-r--r-- | R/TOXSWA_cwa.R | 2 | ||||
| -rw-r--r-- | R/endpoint.R | 2 | ||||
| -rw-r--r-- | R/twa.R | 25 |
7 files changed, 67 insertions, 36 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..4b3111e 100644 --- a/R/PEC_sw_drainage_UK.R +++ b/R/PEC_sw_drainage_UK.R @@ -1,7 +1,8 @@ #' 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 before the start 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 interception The fraction of the application rate that does not reach the soil @@ -24,6 +25,9 @@ #' @author Johannes Ranke #' @examples #' PEC_sw_drainage_UK(150, Koc = 100) +#' 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) @@ -34,20 +38,26 @@ PEC_sw_drainage_UK <- function(rate, interception = 0, Koc, 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") { + ref_year <- 2000 + } else { ref_year <- 1999} + 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") + degradation_time <- as.numeric(difftime(as.Date(paste0(ref_year,"-10-01")), units = "days", latest)) + if (degradation_time > 0) { + 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") - } - if (!missing(model)) { - fraction_left <- pfm_degradation(model, parms = model_parms, - times = degradation_time)[1, "parent"] - amount_available <- fraction_left * amount_available + 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 + } } } diff --git a/R/PEC_sw_drift.R b/R/PEC_sw_drift.R index 7028101..5c7fff4 100644 --- a/R/PEC_sw_drift.R +++ b/R/PEC_sw_drift.R @@ -21,7 +21,7 @@ #' 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'. #' @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 diff --git a/R/PEC_sw_exposit.R b/R/PEC_sw_exposit.R index fe23284..282a1ac 100644 --- a/R/PEC_sw_exposit.R +++ b/R/PEC_sw_exposit.R @@ -79,7 +79,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 +146,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 +174,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 +217,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..b132beb 100644 --- a/R/TOXSWA_cwa.R +++ b/R/TOXSWA_cwa.R @@ -150,7 +150,7 @@ plot.TOXSWA_cwa <- function(x, time_column = c("datetime", "t", "t_firstjan", "t #' by \code{\link{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 |