From 12a31f4c130c551f82232d9ef7dfb608bd52c53f Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Tue, 27 Sep 2016 23:00:48 +0200 Subject: Reorganise repository using standard package layout --- R/FOCUS_GW_scenarios_2012.R | 11 ++ R/GUS.R | 98 ++++++++++ R/PEC_soil.R | 215 +++++++++++++++++++++ R/PEC_sw_drainage_UK.R | 66 +++++++ R/PEC_sw_drift.R | 65 +++++++ R/PEC_sw_sed.R | 50 +++++ R/SFO_actual_twa.R | 36 ++++ R/SSLRC_mobility_classification.R | 42 +++++ R/TOXSWA_cwa.R | 380 ++++++++++++++++++++++++++++++++++++++ R/drift_data_JKI.R | 44 +++++ R/endpoint.R | 109 +++++++++++ R/geomean.R | 38 ++++ R/pfm_degradation.R | 48 +++++ R/soil_scenario_data_EFSA_2015.R | 40 ++++ 14 files changed, 1242 insertions(+) create mode 100644 R/FOCUS_GW_scenarios_2012.R create mode 100644 R/GUS.R create mode 100644 R/PEC_soil.R create mode 100644 R/PEC_sw_drainage_UK.R create mode 100644 R/PEC_sw_drift.R create mode 100644 R/PEC_sw_sed.R create mode 100644 R/SFO_actual_twa.R create mode 100644 R/SSLRC_mobility_classification.R create mode 100644 R/TOXSWA_cwa.R create mode 100644 R/drift_data_JKI.R create mode 100644 R/endpoint.R create mode 100644 R/geomean.R create mode 100644 R/pfm_degradation.R create mode 100644 R/soil_scenario_data_EFSA_2015.R (limited to 'R') diff --git a/R/FOCUS_GW_scenarios_2012.R b/R/FOCUS_GW_scenarios_2012.R new file mode 100644 index 0000000..49cf7dd --- /dev/null +++ b/R/FOCUS_GW_scenarios_2012.R @@ -0,0 +1,11 @@ +#' A very small subset of the FOCUS Groundwater scenario defitions +#' +#' Currently, only a small subset of the soil definitions are provided. +#' +#' @name FOCUS_GW_scenarios_2012 +#' @references FOCUS (2012) Generic guidance for Tier 1 FOCUS ground water assessments. Version 2.1. +#' FOrum for the Co-ordination of pesticde fate models and their USe. +#' http://focus.jrc.ec.europa.eu/gw/docs/Generic_guidance_FOCV2_1.pdf +#' @examples +#' FOCUS_GW_scenarios_2012 +NULL diff --git a/R/GUS.R b/R/GUS.R new file mode 100644 index 0000000..4a6532d --- /dev/null +++ b/R/GUS.R @@ -0,0 +1,98 @@ +# Copyright (C) 2015 Johannes Ranke +# Contact: jranke@uni-bremen.de +# This file is part of the R package pfm + +# This program is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. + +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. + +# You should have received a copy of the GNU General Public License along with +# this program. If not, see + +#' Groundwater ubiquity score based on Gustafson (1989) +#' +#' The groundwater ubiquity score GUS is calculated according to +#' the following 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 +#' method for assessing pesticide leachability. \emph{Environmental +#' toxicology and chemistry} \bold{8}(4) 339–57. +#' @inheritParams endpoint +#' @param DT50 Half-life of the chemical in soil. Should be a field +#' half-life according to Gustafson (1989). However, leaching to the sub-soil +#' can not completely be excluded in field dissipation experiments and Gustafson +#' did not refer to any normalisation procedure, but says the field study should +#' be conducted under use conditions. +#' @param Koc The sorption constant normalised to organic carbon. Gustafson +#' does not mention the nonlinearity of the sorption constant commonly +#' found and usually described by Freundlich sorption, therefore it is +#' unclear at which reference concentration the Koc should be observed +#' (and if the reference concentration would be in soil or in porewater). +#' @param chent If a chent is given with appropriate information present in its +#' chyaml field, this information is used, with defaults specified below. +#' @param degradation_value Which of the available degradation values should +#' be used? +#' @param lab_field Should laboratory or field half-lives be used? This +#' defaults to lab in this implementation, in order to avoid +#' double-accounting for mobility. If comparability with the original GUS +#' values given by Gustafson (1989) is desired, non-normalised first-order +#' field half-lives obtained under actual use conditions should be used. +#' @param redox Aerobic or anaerobic degradation data +#' @param sorption_value Which of the available sorption values should be used? +#' Defaults to Kfoc as this is what is generally available from the European +#' pesticide peer review process. These values generally use a reference +#' concentration of 1 mg/L in porewater, that means they would be expected to +#' be Koc values at a concentration of 1 mg/L in the water phase. +#' @param degradation_aggregator Function for aggregating half-lives +#' @param sorption_aggregator Function for aggregation Koc values +#' @param ... Included in the generic to allow for further arguments later. Therefore +#' this also had to be added to the specific methods. +#' @return A list with the DT50 and Koc used as well as the resulting score +#' of class GUS_result +#' @author Johannes Ranke +#' @export +GUS <- function(...) UseMethod("GUS") + +#' @rdname GUS +#' @export +GUS.numeric <- function(DT50, Koc, ...) { + score <- log10(DT50) * (4 - log10(Koc)) + res <- list(DT50 = DT50, Koc = Koc, score = score) + class(res) <- "GUS_result" + return(res) +} + +#' @rdname GUS +#' @export +GUS.chent <- function(chent, + degradation_value = "DT50ref", + lab_field = "laboratory", + redox = "aerobic", + sorption_value = "Kfoc", + degradation_aggregator = geomean, + sorption_aggregator = geomean, + ...) +{ + DT50 = soil_DT50(chent, lab_field = lab_field, redox = redox, + value = degradation_value, + aggregator = degradation_aggregator, signif = 5) + Koc = soil_Kfoc(chent, value = sorption_value, + aggregator = sorption_aggregator, signif = 5) + GUS.numeric(DT50, Koc) +} + +#' @rdname GUS +#' @export +#' @param x An object of class GUS_result to be printed +#' @param digits The number of digits used in the print method +print.GUS_result = function(x, ..., digits = 1) { + cat("GUS: ", round(x$score, digits = 1), "\n") + cat("calculated from DT50 ", x$DT50, " and Koc ", x$Koc, "\n") +} diff --git a/R/PEC_soil.R b/R/PEC_soil.R new file mode 100644 index 0000000..0263e47 --- /dev/null +++ b/R/PEC_soil.R @@ -0,0 +1,215 @@ +# Copyright (C) 2015 Johannes Ranke +# Contact: jranke@uni-bremen.de +# This file is part of the R package pfm + +# This program is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. + +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. + +# You should have received a copy of the GNU General Public License along with +# this program. If not, see + +# Register global variables +if(getRversion() >= '2.15.1') utils::globalVariables(c("destination", "study_type", "TP_identifier", + "soil_scenario_data_EFSA_2015")) + +#' Calculate predicted environmental concentrations in soil +#' +#' This is a basic calculation of a contaminant concentration in bulk soil +#' based on complete, instantaneous mixing. If an interval is given, an +#' attempt is made at calculating a long term maximum concentration using +#' the concepts layed out for example in the PPR panel opinion (EFSA 2012). +#' +#' This assumes that the complete load to soil during the time specified by +#' 'interval' (typically 365 days) is dosed at once. As in the PPR panel +#' opinion cited below (PPR panel 2012), only temperature correction using the +#' Arrhenius equation is performed. +#' +#' Total soil and porewater PEC values for the scenarios as defined in the EFSA +#' guidance (2015, p. 13) can easily be calculated. +#' +#' @note If temperature information is available in the selected scenarios, as +#' e.g. in the EFSA scenarios, the DT50 for groundwater modelling +#' (destination 'PECgw') is taken from the chent object, otherwise the DT50 +#' with destination 'PECsoil'. +#' @importFrom methods is +#' @param rate Application rate in units specified below +#' @param rate_units Defaults to g/ha +#' @param interception The fraction of the application rate that does not reach the soil +#' @param mixing_depth Mixing depth in cm +#' @param interval Period of the deeper mixing, defaults to 365, which is a year if +#' rate units are in days +#' @param n_periods Number of periods to be considered for long term PEC calculations +#' @param PEC_units Requested units for the calculated PEC. Only mg/kg currently supported +#' @param PEC_pw_units Only mg/L currently supported +#' @param tillage_depth Periodic (see interval) deeper mixing in cm +#' @param chent An optional chent object holding substance specific information. Can +#' also be a name for the substance as a character string +#' @param DT50 If specified, overrides soil DT50 endpoints from a chent object +#' If DT50 is not specified here and not available from the chent object, zero +#' degradation is assumed +#' @param Koc If specified, overrides Koc endpoints from a chent object +#' @param Kom Calculated from Koc by default, but can explicitly be specified +#' as Kom here +#' @param t_avg Averaging times for time weighted average concentrations +#' @param scenarios If this is 'default', 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 +#' 'chents' is specified, the DegT50 with destination 'PECgw' will be used), +#' and corrected using an Arrhenius activation energy of 65.4 kJ/mol. Also +#' model and scenario adjustment factors from the EFSA guidance are used. +#' @param porewater Should equilibrium porewater concentrations be estimated +#' based on Kom and the organic carbon fraction of the soil instead of total +#' soil concentrations? Based on equation (7) given in the PPR panel opinion +#' (EFSA 2012, p. 24) and the scenarios specified in the EFSA guidance (2015, +#' p. 13). +#' @return The predicted concentration in soil +#' @references EFSA Panel on Plant Protection Products and their Residues (2012) +#' Scientific Opinion on the science behind the guidance for scenario +#' selection and scenario parameterisation for predicting environmental +#' concentrations of plant protection products in soil. \emph{EFSA Journal} +#' \bold{10}(2) 2562, doi:10.2903/j.efsa.2012.2562 +#' +#' EFSA (European Food Safety Authority) (2015) EFSA guidance document for +#' predicting environmental concentrations of active substances of plant +#' protection products and transformation products of these active substances +#' in soil. \emph{EFSA Journal} \bold{13}(4) 4093 +#' doi:10.2903/j.efsa.2015.4093 +#' @author Johannes Ranke +#' @export +#' @examples +#' PEC_soil(100, interception = 0.25) +#' +#' # This is example 1 starting at p. 79 of the EFSA guidance (2015) +#' PEC_soil(1000, interval = 365, DT50 = 250, t_avg = c(0, 21), +#' scenarios = "EFSA_2015") +#' PEC_soil(1000, interval = 365, DT50 = 250, t_av = c(0, 21), +#' Kom = 1000, scenarios = "EFSA_2015", porewater = TRUE) +#' +#' # The following is from example 4 starting at p. 85 of the EFSA guidance (2015) +#' # Metabolite M2 +#' # Calculate total and porewater soil concentrations for tier 1 scenarios +#' # Relative molar mass is 100/300, formation fraction is 0.7 * 1 +#' results_pfm <- PEC_soil(100/300 * 0.7 * 1 * 1000, interval = 365, DT50 = 250, t_avg = c(0, 21), +#' scenarios = "EFSA_2015") +#' results_pfm_pw <- PEC_soil(100/300 * 0.7 * 1000, interval = 365, DT50 = 250, t_av = c(0, 21), +#' Kom = 100, scenarios = "EFSA_2015", porewater = TRUE) + +PEC_soil <- function(rate, rate_units = "g/ha", interception = 0, + mixing_depth = 5, + PEC_units = "mg/kg", PEC_pw_units = "mg/L", + interval = NA, n_periods = Inf, + tillage_depth = 20, + chent = NA, + DT50 = NA, + Koc = NA, Kom = Koc / 1.724, + t_avg = 0, + scenarios = c("default", "EFSA_2015"), + porewater = FALSE) +{ + rate_to_soil = (1 - interception) * rate + rate_units = match.arg(rate_units) + PEC_units = match.arg(PEC_units) + scenarios = match.arg(scenarios) + sce <- switch(scenarios, + default = data.frame(rho = 1.5, T_arr = NA, theta_fc = 0.2, f_oc = 0.02, + f_sce = 1, f_mod = 1, row.names = "default"), + EFSA_2015 = if (porewater) soil_scenario_data_EFSA_2015[4:6, ] + else soil_scenario_data_EFSA_2015[1:3, ] + ) + n_sce = nrow(sce) + + soil_volume = 100 * 100 * (mixing_depth/100) # in m3 + soil_mass = soil_volume * sce$rho * 1000 # in kg + + # The following is C_T,ini from EFSA 2012, p. 22, but potentially with interception > 0 + PEC_soil_ini = rate_to_soil * 1000 / soil_mass # in mg/kg + + # Decide which DT50 to take, or set degradation to zero if no DT50 available + if (is.na(DT50) & is(chent, "chent")) { + if (all(is.na(sce$T_arr))) { # No temperature correction + DT50 <- subset(chent$soil_degradation_endpoints, destination == "PECsoil")$DT50 + } else { + DT50 <- subset(chent$soil_degradation_endpoints, destination == "PECgw")$DT50 + } + if (length(DT50) > 1) stop("More than one PECsoil DT50 in chent object") + if (length(DT50) == 0) DT50 <- Inf + } + k = log(2)/DT50 + + # Temperature correction of degradation (accumulation) + if (all(is.na(sce$T_arr))) { # No temperature correction + f_T = 1 + } else { + # Temperature correction as in EFSA 2012 p. 23 + f_T = ifelse(sce$T_arr == 0, + 0, + exp(- (65.4 / 0.008314) * (1/(sce$T_arr + 273.15) - 1/293.15))) + } + + # X is the fraction left after one period (EFSA guidance p. 23) + X = exp(- k * f_T * interval) + + # f_accu is the fraction left after n periods (X + X^2 + ...) + f_accu = 0 + if (!is.na(interval)) { + if (n_periods == Inf) { + f_accu = X/(1 - X) + } else { + for (i in 1:n_periods) { + f_accu = f_accu + X^i + } + } + } + + f_tillage = mixing_depth / tillage_depth + + PEC_background = f_accu * f_tillage * PEC_soil_ini + + PEC_soil = (1 + f_accu * f_tillage) * PEC_soil_ini + + # Get porewater PEC if requested + if (porewater) { + + # If Kom is not specified, try to get K(f)oc + if (is.na(Kom)) { + # If Koc not specified, try to get K(f)oc from chent + if (is.na(Koc) & is(chent, "chent")) { + Koc <- soil_Kfoc(chent) + } + Kom <- Koc / 1.724 + } + + if (is.na(Kom)) stop("No Kom information specified") + + PEC_soil = PEC_soil/((sce$theta_fc/sce$rho) + sce$f_om * Kom) + } + + # Scenario adjustment factors + PEC_soil_sce = PEC_soil * sce$f_sce + + # Model adjustment factors + PEC_soil_sce_mod = PEC_soil_sce * sce$f_mod + + result <- matrix(NA, ncol = n_sce, nrow = length(t_avg), + dimnames = list(t_avg = t_avg, scenario = rownames(sce))) + + result[1, ] <- PEC_soil_sce_mod + + for (i in seq_along(t_avg)) { + t_av_i <- t_avg[i] + if (t_av_i > 0) { + # Equation 10 from p. 24 (EFSA 2015) + result[i, ] <- PEC_soil_sce_mod/(t_av_i * f_T * k) * (1 - exp(- f_T * k * t_av_i)) + } + } + + return(result) +} diff --git a/R/PEC_sw_drainage_UK.R b/R/PEC_sw_drainage_UK.R new file mode 100644 index 0000000..e53f179 --- /dev/null +++ b/R/PEC_sw_drainage_UK.R @@ -0,0 +1,66 @@ +# Copyright (C) 2015 Johannes Ranke +# Contact: jranke@uni-bremen.de +# This file is part of the R package pfm + +# This program is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. + +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. + +# You should have received a copy of the GNU General Public License along with +# this program. If not, see + +#' 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 +#' +#' @param rate Application rate in g/ha +#' @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 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 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 +#' @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) +{ + percentage_lost <- SSLRC_mobility_classification(Koc)[[2]] + amount_available <- rate * (1 - interception) # g/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") + 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") + } + 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 + return(PEC) +} diff --git a/R/PEC_sw_drift.R b/R/PEC_sw_drift.R new file mode 100644 index 0000000..261706c --- /dev/null +++ b/R/PEC_sw_drift.R @@ -0,0 +1,65 @@ +# Copyright (C) 2015 Johannes Ranke +# Contact: jranke@uni-bremen.de +# This file is part of the R package pfm + +# This program is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. + +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. + +# You should have received a copy of the GNU General Public License along with +# this program. If not, see + +#' 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. +#' +#' @param rate Application rate in units specified below +#' @param applications Number of applications for selection of drift percentile +#' @param drift_percentages Percentage drift values for which to calculate PECsw. +#' 'drift_data' and 'distances' if not NULL. +#' @param drift_data Source of drift percentage data +#' @param crop Crop name (use German names for JKI data), defaults to "Ackerbau" +#' @param distances The distances in m for which to get PEC values +#' @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 +#' @return The predicted concentration in surface water +#' @export +#' @author Johannes Ranke +#' @examples +#' PEC_sw_drift(100) +PEC_sw_drift <- function(rate, + applications = 1, + water_depth = 30, + drift_percentages = NULL, + drift_data = "JKI", + crop = "Ackerbau", + distances = c(1, 5, 10, 20), + rate_units = "g/ha", + PEC_units = "\u00B5g/L") +{ + rate_units <- match.arg(rate_units) + PEC_units <- match.arg(PEC_units) + drift_data <- match.arg(drift_data) + water_volume <- 100 * 100 * (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 (is.null(drift_percentages)) { + drift_percentages <- pfm::drift_data_JKI[[applications]][dist_index, crop] + names(drift_percentages) <- paste(dist_index, "m") + } else { + names(drift_percentages) <- paste(drift_percentages, "%") + } + + PEC_sw_drift <- PEC_sw_overspray * drift_percentages / 100 + return(PEC_sw_drift) +} diff --git a/R/PEC_sw_sed.R b/R/PEC_sw_sed.R new file mode 100644 index 0000000..a4a3d6a --- /dev/null +++ b/R/PEC_sw_sed.R @@ -0,0 +1,50 @@ +# Copyright (C) 2015 Johannes Ranke +# Contact: jranke@uni-bremen.de +# This file is part of the R package pfm + +# This program is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. + +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. + +# You should have received a copy of the GNU General Public License along with +# this program. If not, see + +#' Calculate predicted environmental concentrations in sediment from surface +#' water concentrations +#' +#' The method 'percentage' is equivalent to what is used in the CRD spreadsheet +#' PEC calculator +#' +#' @param PEC_sw Numeric vector or matrix of surface water concentrations in µg/L for +#' which the corresponding sediment concentration is to be estimated +#' @param percentage The percentage in sediment, used for the percentage method +#' @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 +#' 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 +#' 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, + PEC_sed_units = c("\u00B5g/kg", "mg/kg")) +{ + 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 + } + return(PEC_sed) +} diff --git a/R/SFO_actual_twa.R b/R/SFO_actual_twa.R new file mode 100644 index 0000000..7facb6a --- /dev/null +++ b/R/SFO_actual_twa.R @@ -0,0 +1,36 @@ +# Copyright (C) 2015 Johannes Ranke +# Contact: jranke@uni-bremen.de +# This file is part of the R package pfm + +# This program is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. + +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. + +# You should have received a copy of the GNU General Public License along with +# this program. If not, see + +#' Actual and maximum moving window time average concentrations for SFO kinetics +#' +#' @param DT50 The half-life. +#' @param times The output times, and window sizes for time weighted average concentrations +#' @export +#' @author Johannes Ranke +#' @source FOCUS (2014) Generic Guidance for Estimating Persistence and Degradation +#' Kinetics from Environmental Fate Studies on Pesticides in EU Registratin, Version 1.1, +#' 18 December 2014, p. 251 +#' @examples +#' SFO_actual_twa(10) +SFO_actual_twa <- function(DT50 = 1000, times = c(0, 1, 2, 4, 7, 14, 21, 28, 42, 50, 100)) +{ + k = log(2)/DT50 + result <- data.frame(actual = 1 * exp(-k * times), + twa = (1 - exp(-k * times))/(k * times), + row.names = times) + return(result) +} diff --git a/R/SSLRC_mobility_classification.R b/R/SSLRC_mobility_classification.R new file mode 100644 index 0000000..deda5cf --- /dev/null +++ b/R/SSLRC_mobility_classification.R @@ -0,0 +1,42 @@ +# Copyright (C) 2015 Johannes Ranke +# Contact: jranke@uni-bremen.de +# This file is part of the R package pfm + +# This program is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. + +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. + +# You should have received a copy of the GNU General Public License along with +# this program. If not, see + +#' Determine the SSLRC mobility classification for a chemical substance from its Koc +#' +#' This implements the method specified in the UK data requirements handbook and was +#' checked against the spreadsheet published on the CRC website +#' +#' @param Koc The sorption coefficient normalised to organic carbon in L/kg +#' @return A list containing the classification and the percentage of the +#' compound transported per 10 mm drain water +#' @export +#' @author Johannes Ranke +#' @examples +#' SSLRC_mobility_classification(100) +SSLRC_mobility_classification <- function(Koc) +{ + if (!is.numeric(Koc) | length(Koc) != 1) stop("Please give a single number") + result <- list("Non mobile", 0.01) + if (Koc < 4000) result <- list("Slightly mobile", 0.02) + if (Koc < 1000) result <- list("Slightly mobile", 0.5) + if (Koc < 500) result <- list("Moderately mobile", 0.7) + if (Koc < 75) result <- list("Mobile", 1.9) + if (Koc < 15) result <- list("Very mobile", 1.9) + names(result) <- c("Mobility classification", + "Percentage drained per mm of drain water") + return(result) +} diff --git a/R/TOXSWA_cwa.R b/R/TOXSWA_cwa.R new file mode 100644 index 0000000..c5ddce9 --- /dev/null +++ b/R/TOXSWA_cwa.R @@ -0,0 +1,380 @@ +# Copyright (C) 2014,2015,2016 Johannes Ranke +# Contact: jranke@uni-bremen.de +# This file is part of the R package pfm + +# This program is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. + +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. + +# You should have received a copy of the GNU General Public License along with +# this program. If not, see + +#' Read TOXSWA surface water concentrations +#' +#' Read TOXSWA hourly concentrations of a chemical substance in a specific +#' segment of a TOXSWA surface water body. Per default, the data for the last +#' segment are imported. As TOXSWA 4 reports the values at the end of the hour +#' (ConLiqWatLayCur) in its summary file, we use this value as well instead +#' of the hourly averages (ConLiqWatLay). +#' +#' @param filename The filename of the cwa file (TOXSWA 2.x.y or similar) or the +#' out file (FOCUS TOXSWA 4, i.e. TOXSWA 4.4.2 or similar). +#' @param basedir The path to the directory where the cwa file resides. +#' @param zipfile Optional path to a zip file containing the cwa file. +#' @param segment The segment for which the data should be read. Either "last", or +#' the segment number. +#' @param total Set this to TRUE in order to read total concentrations as well. This is +#' only necessary for .out files as generated by TOXSWA 4.4.2 or similar, not for .cwa +#' files. For .cwa files, the total concentration is always read as well. +#' @param substance For TOXSWA 4 .out files, the default value "parent" leads +#' to reading concentrations of the parent compound. Alternatively, the substance +#' of interested can be selected by its code name. +#' @param windows Numeric vector of width of moving windows in days, for calculating +#' maximum time weighted average concentrations and areas under the curve. +#' @param thresholds Numeric vector of threshold concentrations in µg/L for +#' generating event statistics. +#' @importFrom readr read_fwf fwf_empty +#' @return An instance of an R6 object of class +#' \code{\link{TOXSWA_cwa}}. +#' @export +#' @author Johannes Ranke +#' @examples +#' H_sw_D4_pond <- read.TOXSWA_cwa("00001p_pa.cwa", +#' basedir = "SwashProjects/project_H_sw/TOXSWA", +#' zipfile = system.file("testdata/SwashProjects.zip", +#' package = "pfm")) +read.TOXSWA_cwa <- function(filename, basedir = ".", zipfile = NULL, + segment = "last", substance = "parent", + total = FALSE, + windows = NULL, thresholds = NULL) +{ + if (!missing(filename)) { + cwa <- TOXSWA_cwa$new(filename, basedir, zipfile, + substance = substance, total = total) + if (!is.null(windows[1])) cwa$moving_windows(windows) + if (!is.null(thresholds[1])) cwa$get_events(thresholds) + invisible(cwa) + } else { + message("You need to specify a filename for the cwa file to be read") + } +} + +#' Plot TOXSWA surface water concentrations +#' +#' Plot TOXSWA hourly concentrations of a chemical substance in a specific +#' segment of a TOXSWA surface water body. +#' +#' @import graphics +#' @param x The TOXSWA_cwa object to be plotted. +#' @param xlab,ylab Labels for x and y axis. +#' @param time_column What should be used for the time axis. If "t_firstjan" is chosen, +#' the time is given in days relative to the first of January in the first year. +#' @param add Should we add to an existing plot? +#' @param total Should the total concentration in water be plotted, including substance sorbed +#' to suspended matter? +#' @param LC_TIME Specification of the locale used to format dates +#' @param ... Further arguments passed to \code{plot} if we are not adding to an existing plot +#' @export +#' @author Johannes Ranke +#' @examples +#' H_sw_D4_pond <- read.TOXSWA_cwa("00001p_pa.cwa", +#' basedir = "SwashProjects/project_H_sw/TOXSWA", +#' zipfile = system.file("testdata/SwashProjects.zip", +#' package = "pfm")) +#' plot(H_sw_D4_pond) +plot.TOXSWA_cwa <- function(x, time_column = c("datetime", "t", "t_firstjan", "t_rel_to_max"), + xlab = "default", ylab = "default", + add = FALSE, + total = FALSE, LC_TIME = "C", ...) +{ + time_column = match.arg(time_column) + cwa_column = ifelse(total, "cwa_tot_mug_per_L", "cwa_mug_per_L") + lct <- Sys.getlocale("LC_TIME") + tmp <- Sys.setlocale("LC_TIME", LC_TIME) + if (xlab == "default") { + xlab = switch(time_column, + datetime = "Time", + t = "Time [days]", + t_firstjan = "Time since first of January [days]", + t_rel_to_max = "Time relative to maximum concentration [days]") + } + if (ylab == "default") { + ylab = paste( ifelse(total, "Total concentration", "Concentration"), "[\u03bcg/L]") + } + if (add) { + lines(x$cwas[c(time_column, cwa_column)], xlab = xlab, ylab = ylab, ...) + } else{ + if (time_column == "datetime") { + plot(x$cwas$datetime, x$cwas[[cwa_column]], type = "l", + xlab = xlab, ylab = ylab, xaxt = "n", ...) + seq_x <- seq(min(x$cwas$datetime), max(x$cwas$datetime), by = "quarter") + axis(1, at = seq_x, labels = format(seq_x, "%b %Y")) + } else { + plot(x$cwas[c(time_column, cwa_column)], type = "l", + xlab = xlab, ylab = ylab, ...) + } + } + tmp <- Sys.setlocale("LC_TIME", lct) +} + +#' R6 class for holding TOXSWA cwa concentration data and associated statistics +#' +#' An R6 class for holding TOXSWA cwa concentration data and some associated statistics. +#' Usually, an instance of this class will be generated by \code{\link{read.TOXSWA_cwa}}. +#' +#' @docType class +#' @importFrom R6 R6Class +#' @export +#' @format An \code{\link{R6Class}} generator object. +#' @field filename Length one character vector. +#' @field basedir Length one character vector. +#' @field segment Length one integer, specifying for which segment the cwa data were read. +#' @field cwas Dataframe holding the concentrations. +#' @field events List of dataframes holding the event statistics for each threshold. +#' @field windows Matrix of maximum time weighted average concentrations (TWAC_max) +#' and areas under the curve in µg/day * h (AUC_max_h) or µg/day * d (AUC_max_d) +#' for the requested moving window sizes in days. +#' @section Methods: +#' \describe{ +#' \item{\code{get_events(threshold, total = FALSE)}}{ +#' Populate a datataframe with event information for the specified threshold value +#' in µg/L. If \code{total = TRUE}, the total concentration including the amount +#' adsorbed to suspended matter will be used. The resulting dataframe is stored in the +#' \code{events} field of the object. +#' } +#' \item{\code{moving_windows(windows, total = FALSE)}}{ +#' Add to the \code{windows} field described above. +#' Again, if \code{total = TRUE}, the total concentration including the amount +#' adsorbed to suspended matter will be used. +#' } +#' } +#' @examples +#' H_sw_R1_stream <- read.TOXSWA_cwa("00003s_pa.cwa", +#' basedir = "SwashProjects/project_H_sw/TOXSWA", +#' zipfile = system.file("testdata/SwashProjects.zip", +#' package = "pfm")) +#' H_sw_R1_stream$get_events(c(2, 10)) +#' H_sw_R1_stream$moving_windows(c(7, 21)) +#' print(H_sw_R1_stream) +#' @keywords data + +TOXSWA_cwa <- R6Class("TOXSWA_cwa", + public = list( + filename = NULL, + basedir = NULL, + zipfile = NULL, + segment = NULL, + substance = NULL, + cwas = NULL, + windows = NULL, + events = list(), + initialize = function(filename, basedir, zipfile = NULL, + segment = "last", substance = "parent", total = FALSE) { + self$filename <- filename + self$basedir <- basedir + self$zipfile <- zipfile + if (!is.null(zipfile)) { + try(file_connection <- unz(zipfile, paste0(basedir, "/", filename))) + } else { + try(file_connection <- file(file.path(basedir, filename), "rt")) + } + + + if (grepl(".cwa$", filename)) { + # cwa file from FOCUS TOXSWA 3 (TOXSWA 2.x.y) + cwa_all_segments <- try( + read.table(file_connection, + sep = "", skip = 40, + encoding = "UTF-8", + colClasses = c("character", "numeric", + "integer", rep("numeric", 5)), + col.names = c("datetime", "t", "segment", + "xcd", "cwa_tot", "cwa", "Xss", "Xmp"))) + + if (is.null(zipfile)) close(file_connection) # only needed for files + + if (!inherits(cwa_all_segments, "try-error")) { + available_segments = 1:max(cwa_all_segments$segment) + if (segment == "last") segment = max(available_segments) + if (!segment %in% available_segments) stop("Invalid segment specified") + self$segment <- segment + cwa <- subset(cwa_all_segments, segment == self$segment, + c("datetime", "t", "segment", "cwa", "cwa_tot")) + lct <- Sys.getlocale("LC_TIME"); Sys.setlocale("LC_TIME", "C") + cwa$datetime <- strptime(cwa$datetime, "%d-%b-%Y-%H:%M") + Sys.setlocale("LC_TIME", lct) + startyear = format(cwa$datetime[1], "%Y") + firstjan <- strptime(paste0(startyear, "-01-01"), "%Y-%m-%d") + cwa$t_firstjan <- as.numeric(difftime(cwa$datetime, + firstjan, units = "days")) + + t_max = cwa[which.max(cwa$cwa), "t"] + cwa$t_rel_to_max = cwa$t - t_max + cwa$cwa_mug_per_L <- cwa$cwa * 1000 + cwa$cwa_tot_mug_per_L <- cwa$cwa_tot * 1000 + self$cwas <- cwa[c("datetime", "t", "t_firstjan", + "t_rel_to_max", + "cwa_mug_per_L", + "cwa_tot_mug_per_L")] + } else { + stop("Could not read ", filename) + } + } else { + # out file from FOCUS TOXSWA 4 (TOXSWA 4.4.2 or similar) + outfile <- try(readLines(file_connection)) + + close(file_connection) # only needed for files + + if (inherits(outfile, "try-error")) { + stop("Could not read ", filename) + } else { + # Get the substance name(s) + sub_lines <- grep(".*0.000.*ConLiqWatLayCur_", outfile[1:50], value = TRUE) + substances <- gsub(".*ConLiqWatLayCur_(.*?) *[0-9].*", "\\1", sub_lines) + if (!substance %in% c("parent", substances)) { + stop("No data for substance ", substance, " present in the .out file.") + } + + # Generate field name for the concentrations at end of hour for the + # substance of interest + if (substance == "parent") { + cwa_string = paste0("ConLiqWatLayCur_", substances[1]) + } else { + cwa_string = paste0("ConLiqWatLayCur_", substance) + } + + 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"))) + + # Append time "-00h00" to datetime in first row, as this is not (always?) present + # in the line ConLiqWatLayCur + if (nchar(cwa_all_segments[1, "X2"]) == 11) { + cwa_all_segments[1, "X2"] = paste0(cwa_all_segments[1, "X2"], "-00h00") + } + + available_segments = 1:(ncol(cwa_all_segments) - 3) + if (segment == "last") segment = max(available_segments) + if (!segment %in% available_segments) stop("Invalid segment specified") + self$segment <- segment + cwa <- data.frame( + datetime = as.character(cwa_all_segments$X2), + t = cwa_all_segments$X1, + cwa = cwa_all_segments[[3 + segment]] + ) + if (total) { + cwa_tot_lines <- outfile[grep("ConSysWatLay_", outfile)] # hourly total conc. + cwa_tot_all_segments <- read.table(text = cwa_lines) + cwa$cwa_tot = cwa_tot_all_segments[[3 + segment]] + } + lct <- Sys.getlocale("LC_TIME"); Sys.setlocale("LC_TIME", "C") + cwa$datetime <- strptime(cwa$datetime, "%d-%b-%Y-%Hh%M") + Sys.setlocale("LC_TIME", lct) + + startyear = format(cwa$datetime[1], "%Y") + firstjan <- strptime(paste0(startyear, "-01-01"), "%Y-%m-%d") + cwa$t_firstjan <- as.numeric(difftime(cwa$datetime, + firstjan, units = "days")) + t_max = cwa[which.max(cwa$cwa), "t"] + cwa$t_rel_to_max = cwa$t - t_max + cwa$cwa_mug_per_L <- cwa$cwa * 1000 + self$cwas <- cwa[c("datetime", "t", "t_firstjan", + "t_rel_to_max", + "cwa_mug_per_L")] + if (total) { + self$cwas$cwa_tot_mug_per_L <- cwa$cwa_tot * 1000 + } + } + } + }, + moving_windows = function(windows, total = FALSE) { + window_names = paste(windows, "days") + n = length(window_names) + self$windows <- data.frame(window = window_names, + max_TWAC = numeric(n), + max_AUC_h = numeric(n), + max_AUC_d = numeric(n)) + if (missing(windows)) { + stop("You need to specify at least one moving window size in days") + } + cwa_column = ifelse(total, "cwa_tot_mug_per_L", "cwa_mug_per_L") + for (i in seq_along(windows)) { + window_size = windows[i] + filter_size = window_size * 24 + max_TWAC = max(filter(self$cwas[cwa_column], rep(1/filter_size, filter_size), + "convolution"), na.rm = TRUE) + max_AUC_h = max_TWAC * filter_size + max_AUC_d = max_TWAC * window_size + + self$windows[i, -1] = c(max_TWAC, max_AUC_h, max_AUC_d) + } + invisible(self) + }, + get_events = function(thresholds, total = FALSE) { + if (missing(thresholds)) { + stop("You need to specify at least one threshold concentration in \u03bcg/L") + } + for (threshold in thresholds) { + events = data.frame(t_start = numeric(), cwa_max = numeric(), + duration = numeric(), pre_interval = numeric(), + AUC_h = numeric(), AUC_d = numeric()) + cwa_column = ifelse(total, "cwa_tot_mug_per_L", "cwa_mug_per_L") + + event_end = 0 + event = FALSE + event_max = 0 + event_nr = 0 + n_rows = nrow(self$cwas) + for (i in 1:n_rows) { + cwa_cur = self$cwas[i, cwa_column] + if (event == FALSE) { + if (cwa_cur > threshold) { + event_start = self$cwas[i, "t"] + pre_interval = event_start - event_end + i_start = i + event = TRUE + event_max = cwa_cur + } + } else { + if (cwa_cur > event_max) event_max = cwa_cur + if (cwa_cur < threshold || i == n_rows) { + event_nr = event_nr + 1 + i_end = i + if (i == n_rows) i_end = i + event_end = self$cwas[i_end, "t"] + event_length = event_end - event_start + event_cwas <- self$cwas[i_start:(i_end - 1), cwa_column] + event_AUC_h = sum(event_cwas) + event_AUC_d = event_AUC_h / 24 + events[event_nr, ] = c(event_start, event_max, event_length, + pre_interval, event_AUC_h, event_AUC_d) + event = FALSE + } + } + } + + self$events[[as.character(threshold)]] <- events + } + invisible(self) + }, + print = function() { + cat(" data from file", self$filename, "segment", self$segment, "\n") + print(head(self$cwas)) + cat("Moving window analysis\n") + print(self$windows) + for (threshold in names(self$events)) { + cat("Event statistics for threshold", threshold, "\n") + if (nrow(self$events[[threshold]]) == 0) cat("No events found\n") + else print(self$events[[threshold]]) + } + } + ) +) +# vim: set ts=2 sw=2 expandtab: diff --git a/R/drift_data_JKI.R b/R/drift_data_JKI.R new file mode 100644 index 0000000..888d3b5 --- /dev/null +++ b/R/drift_data_JKI.R @@ -0,0 +1,44 @@ +#' Deposition from spray drift expressed as percent of the applied dose as +#' published by the JKI +#' +#' Deposition from spray drift expressed as percent of the applied dose as +#' published by the German Julius-Kühn Institute (JKI). +#' +#' The data were extracted from the spreadsheet cited below using the R code +#' given in the example section. The spreadsheet is not included in the package +#' as its licence is not clear. +#' +#' +#' @name drift_data_JKI +#' @docType data +#' @format A list currently containing matrices with spray drift percentage +#' data for field crops (Ackerbau), and Pome/stone fruit, early and late +#' (Obstbau frueh, spaet). +#' @source JKI (2010) Spreadsheet 'Tabelle der Abdrifteckwerte.xls', retrieved +#' from +#' http://www.jki.bund.de/no_cache/de/startseite/institute/anwendungstechnik/abdrift-eckwerte.html +#' on 2015-06-11 +#' @keywords datasets +#' @examples +#' +#' \dontrun{ +#' # This is the code that was used to extract the data +#' library(readxl) +#' abdrift_path <- "inst/extdata/Tabelle der Abdrifteckwerte.xls" +#' JKI_crops <- c("Ackerbau", "Obstbau frueh", "Obstbau spaet") +#' names(JKI_crops) <- c("Field crops", "Pome/stone fruit, early", "Pome/stone fruit, late") +#' drift_data_JKI <- list() +#' +#' for (n in 1:8) { +#' drift_data_raw <- read_excel(abdrift_path, sheet = n + 1, skip = 2) +#' drift_data <- as.matrix(drift_data_raw[1:9, 2:4]) +#' dimnames(drift_data) <- list(distance = as.integer(drift_data_raw[1:9, 1]), +#' crop = JKI_crops) +#' drift_data_JKI[[n]] <- drift_data +#' } +#' save(drift_data_JKI, file = "data/drift_data_JKI.RData") +#' } +#' +#' # And this is the resulting data +#' drift_data_JKI +NULL diff --git a/R/endpoint.R b/R/endpoint.R new file mode 100644 index 0000000..56a314b --- /dev/null +++ b/R/endpoint.R @@ -0,0 +1,109 @@ +#' Retrieve endpoint information from the chyaml field of a chent object +#' +#' R6 class objects of class chent represent chemical entities +#' and can hold a list of information loaded from a chemical yaml file in their +#' chyaml field. Such information is extracted and optionally aggregated by +#' this function. +#' +#' The functions \code{soil_*} are functions to extract soil specific endpoints. +#' For the Freundlich exponent, the capital letter \code{N} is used in order to +#' facilitate dealing with such data in R. In pesticide fate modelling, this +#' exponent is often called 1/n. +#' +#' @export +#' @param chent The chent object to get the information from +#' @param medium The medium for which information is sought +#' @param type The information type +#' @param lab_field If not NA, do we want laboratory or field endpoints +#' @param redox If not NA, are we looking for aerobic or anaerobic data +#' @param value The name of the value we want. The list given in the +#' usage section is not exclusive +#' @param aggregator The aggregator function. Can be mean, +#' \code{\link{geomean}}, or identity, for example. +#' @param raw Should the number(s) be returned as stored in the chent +#' object (could be a character value) to retain original information +#' about precision? +#' @param signif How many significant digits do we want +#' @return The result from applying the aggregator function to +#' the values converted to a numeric vector, rounded to the +#' given number of significant digits, or, if raw = TRUE, +#' the values as a character value, retaining any implicit +#' information on precision that may be present. +#' +endpoint <- function(chent, + medium = "soil", + type = c("degradation", "sorption"), + lab_field = c(NA, "laboratory", "field"), + redox = c(NA, "aerobic", "anaerobic"), + value = c("DT50ref", "Kfoc", "N"), + aggregator = geomean, + raw = FALSE, + signif = 3) +{ + if (!is(chent, "chent")) { + stop("Please supply a chent object as created using the package 'chents' available from jrwb.de") + } + ep_list <- chent$chyaml[[medium]][[type]] + if (!is.na(lab_field[1])) { + ep_list <- ep_list[[lab_field]] + } + if (!is.na(redox[1])) { + ep_list <- ep_list[[redox]] + } + values <- ep_list$data[[value]] + if (raw) return(values) + else return(signif(aggregator(as.numeric(values)), signif)) +} + +#' @inheritParams endpoint +#' @rdname endpoint +#' @export +soil_DT50 <- function(chent, aggregator = geomean, signif = 3, + lab_field = "laboratory", value = "DT50ref", + redox = "aerobic", raw = FALSE) { + ep <- endpoint(chent, medium = "soil", type = "degradation", + lab_field = "laboratory", redox = redox, + signif = signif, + value = value, aggregator = aggregator, raw = raw) + return(ep) +} + +#' @inheritParams endpoint +#' @rdname endpoint +#' @export +soil_Kfoc <- function(chent, aggregator = geomean, signif = 3, + value = "Kfoc", raw = FALSE) { + ep <- endpoint(chent, medium = "soil", type = "sorption", + signif = signif, + value = value, aggregator = aggregator, raw = raw) + return(ep) +} + +#' @inheritParams endpoint +#' @rdname endpoint +#' @export +soil_N <- function(chent, aggregator = mean, signif = 3, raw = FALSE) { + ep <- endpoint(chent, medium = "soil", type = "sorption", + signif = signif, + value = "N", aggregator = aggregator, raw = raw) + return(ep) +} + +#' @inheritParams endpoint +#' @rdname endpoint +#' @param values The values to be returned +#' @param aggregators A named vector of aggregator functions to be used +#' @export +soil_sorption <- function(chent, values = c("Kfoc", "N"), + aggregators = c(Kfoc = geomean, Koc = geomean, N = mean), + signif = c(Kfoc = 3, N = 3), + raw = FALSE) { + res <- sapply(values, + function(x) { + endpoint(chent, medium = "soil", type = "sorption", + signif = signif[[x]], + value = x, aggregator = aggregators[[x]], raw = raw) + } + ) + return(res) +} diff --git a/R/geomean.R b/R/geomean.R new file mode 100644 index 0000000..626829b --- /dev/null +++ b/R/geomean.R @@ -0,0 +1,38 @@ +# Copyright (C) 2015 Johannes Ranke +# Contact: jranke@uni-bremen.de +# This file is part of the R package pfm + +# This program is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. + +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. + +# You should have received a copy of the GNU General Public License along with +# this program. If not, see + +#' Calculate the geometric mean +#' +#' Based on some posts in a thread on Stackoverflow +#' \url{http://stackoverflow.com/questions/2602583/geometric-mean-is-there-a-built-in} +#' This function checks for negative values, removes NA values per default and +#' returns 0 if at least one element of the vector is 0. +#' +#' @param x Vector of numbers +#' @param na.rm Should NA values be omitted? +#' @return The geometric mean +#' @export +#' @author Johannes Ranke +#' @examples +#' geomean(c(1, 3, 9)) +#' geomean(c(1, 3, NA, 9)) +#' \dontrun{geomean(c(1, -3, 9)) # returns an error} +geomean = function(x, na.rm = TRUE){ + if (any(is.na(x)) & na.rm == FALSE) stop("Removal of NA values was not requested") + if (any(x < 0, na.rm = na.rm)) stop("Only defined for positive numbers") + exp(mean(log(x), na.rm = na.rm)) +} diff --git a/R/pfm_degradation.R b/R/pfm_degradation.R new file mode 100644 index 0000000..6c8610e --- /dev/null +++ b/R/pfm_degradation.R @@ -0,0 +1,48 @@ +# Copyright (C) 2015 Johannes Ranke +# Contact: jranke@uni-bremen.de +# This file is part of the R package pfm + +# This program is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. + +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. + +# You should have received a copy of the GNU General Public License along with +# this program. If not, see + +#' Calculate a time course of relative concentrations based on an mkinmod model +#' +#' @import mkin +#' @param model The degradation model to be used. Either a parent only model like +#' 'SFO' or 'FOMC', or an mkinmod object +#' @param DT50 The half-life. This is only used when simple exponential decline +#' is calculated (SFO model). +#' @param parms The parameters used for the degradation model +#' @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 +#' @export +#' @author Johannes Ranke +#' @examples +#' head(pfm_degradation("SFO", DT50 = 10)) +pfm_degradation <- function(model = "SFO", DT50 = 1000, parms = c(k_parent_sink = 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)) + } + initial_state = c(1, rep(0, length(model$diffs) - 1)) + names(initial_state) <- names(model$diffs) + time_course <- mkinpredict(model, odeparms = parms, + odeini = initial_state, + outtimes = times, + solution_type = ifelse(length(model$spec) == 1, + "analytical", "deSolve")) + invisible(time_course) +} diff --git a/R/soil_scenario_data_EFSA_2015.R b/R/soil_scenario_data_EFSA_2015.R new file mode 100644 index 0000000..fb096bc --- /dev/null +++ b/R/soil_scenario_data_EFSA_2015.R @@ -0,0 +1,40 @@ +#' Properties of the predefined scenarios from the EFSA guidance from 2015 +#' +#' Properties of the predefined scenarios used at Tier 1, Tier 2A and Tier 3A for the +#' concentration in soil as given in the EFSA guidance (2015, p. 13/14). Also, the +#' scenario and model adjustment factors from p. 15 and p. 17 are included. +#' +#' @name soil_scenario_data_EFSA_2015 +#' @docType data +#' @format A data frame with one row for each scenario. Row names are the scenario codes, +#' e.g. CTN for the Northern scenario for the total concentration in soil. Columns are +#' mostly self-explanatory. \code{rho} is the dry bulk density of the top soil. +#' @source EFSA (European Food Safety Authority) (2015) +#' EFSA guidance document for predicting environmental concentrations +#' of active substances of plant protection products and transformation products of these +#' active substances in soil. \emph{EFSA Journal} \bold{13}(4) 4093 +#' doi:10.2903/j.efsa.2015.4093 +#' @keywords datasets +#' @examples +#' \dontrun{ +#' # This is the code that was used to define the data +#' soil_scenario_data_EFSA_2015 <- data.frame( +#' Zone = rep(c("North", "Central", "South"), 2), +#' Country = c("Estonia", "Germany", "France", "Denmark", "Czech Republik", "Spain"), +#' T_arit = c(4.7, 8.0, 11.0, 8.2, 9.1, 12.8), +#' T_arr = c(7.0, 10.1, 12.3, 9.8, 11.2, 14.7), +#' Texture = c("Coarse", "Coarse", "Medium fine", "Medium", "Medium", "Medium"), +#' f_om = c(0.118, 0.086, 0.048, 0.023, 0.018, 0.011), +#' theta_fc = c(0.244, 0.244, 0.385, 0.347, 0.347, 0.347), +#' rho = c(0.95, 1.05, 1.22, 1.39, 1.43, 1.51), +#' f_sce = c(3, 2, 2, 2, 1.5, 1.5), +#' f_mod = c(2, 2, 2, 4, 4, 4), +#' stringsAsFactors = FALSE, +#' row.names = c("CTN", "CTC", "CTS", "CLN", "CLC", "CLS") +#' ) +#' save(soil_scenario_data_EFSA_2015, file = '../data/soil_scenario_data_EFSA_2015.RData') +#' } +#' +#' # And this is the resulting dataframe +#' soil_scenario_data_EFSA_2015 +NULL -- cgit v1.2.1