diff options
| -rw-r--r-- | ChangeLog | 17 | ||||
| -rw-r--r-- | DESCRIPTION | 7 | ||||
| -rw-r--r-- | NAMESPACE | 15 | ||||
| -rw-r--r-- | R/twa.R | 220 | ||||
| -rw-r--r-- | man/max_twa.Rd | 19 | ||||
| -rw-r--r-- | man/one_box.Rd | 45 | ||||
| -rw-r--r-- | man/plot.one_box.Rd | 35 | ||||
| -rw-r--r-- | man/sawtooth.Rd | 39 | ||||
| -rw-r--r-- | man/twa.Rd | 29 |
9 files changed, 423 insertions, 3 deletions
@@ -1,3 +1,20 @@ +commit bba2cf3a70849ba86f37520d3e909cf1c706f416 +Author: Johannes Ranke <jranke@uni-bremen.de> +Date: 2016-12-22 11:06:57 +0100 + + Fix reading in times from .out files + + The code from the previous commit was broken. Also, the time + zone for the times that are read is now wet to 'UTC', in order to + avoid setting different time zones due to daylight savings, which + introduces artificial one-hour offsets on changeover days. + +commit 0af7c7b8c34067fc4756929925239c329b28ed32 +Author: Johannes Ranke <jranke@uni-bremen.de> +Date: 2016-12-14 18:26:14 +0100 + + Changelog update and roxygen run + commit 5a04ad3061c1484b45703e44149f49ec97cfbf15 Author: Johannes Ranke <jranke@uni-bremen.de> Date: 2016-12-14 16:52:14 +0100 diff --git a/DESCRIPTION b/DESCRIPTION index 7207ca4..f92e6da 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: pfm Type: Package Title: Utilities for Pesticide Fate Modelling -Version: 0.3-9 -Date: 2016-12-22 +Version: 0.4-1 +Date: 2017-01-18 Authors@R: person("Johannes Ranke", email = "jranke@uni-bremen.de", role = c("aut", "cre", "cph")) Description: Utilities for simple calculations of predicted environmental @@ -17,7 +17,8 @@ Imports: methods Suggests: testthat, - chents + chents, + magrittr License: GPL LazyLoad: yes LazyData: yes @@ -2,8 +2,13 @@ S3method(GUS,chent) S3method(GUS,numeric) +S3method(max_twa,one_box) +S3method(one_box,mkinfit) +S3method(one_box,numeric) S3method(plot,TOXSWA_cwa) +S3method(plot,one_box) S3method(print,GUS_result) +S3method(twa,one_box) export(GUS) export(PEC_soil) export(PEC_sw_drainage_UK) @@ -14,15 +19,25 @@ export(SSLRC_mobility_classification) export(TOXSWA_cwa) export(endpoint) export(geomean) +export(max_twa) +export(one_box) export(pfm_degradation) export(read.TOXSWA_cwa) +export(sawtooth) export(soil_DT50) export(soil_Kfoc) export(soil_N) export(soil_sorption) +export(twa) import(graphics) import(mkin) importFrom(R6,R6Class) importFrom(methods,is) +importFrom(mkin,mkinpredict) importFrom(readr,fwf_empty) importFrom(readr,read_fwf) +importFrom(stats,filter) +importFrom(stats,frequency) +importFrom(stats,plot.ts) +importFrom(stats,time) +importFrom(stats,ts) @@ -0,0 +1,220 @@ +# Copyright (C) 2017 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 <http://www.gnu.org/licenses/> + +#' Create a time series of decline data +#' +#' The time series starts with the amount specified for the first application. +#' This does not create objects of type \code{\link{ts}}. +#' +#' @param x When numeric, this is the half-life to be used for an exponential +#' decline. If x is an mkinfit object, the decline is calculated from this object +#' @param t_end End of the time series +#' @param res Resolution of the time series +#' @param ... Further arguments passed to methods +#' @importFrom stats filter frequency time ts +#' @export +#' @examples +#' # Only use a half-life +#' pred_0 <- one_box(10) +#' plot(pred_0) +#' +#' # Use a fitted mkinfit model +#' require(mkin) +#' fit <- mkinfit("FOMC", FOCUS_2006_C) +#' pred_1 <- one_box(fit) +#' plot(pred_1) +#' +#' # Use a model with more than one observed variable +#' m_2 <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO")) +#' fit_2 <- mkinfit(m_2, FOCUS_2006_D) +#' pred_2 <- one_box(fit_2) +#' plot(pred_2) +one_box <- function(x, + t_end = 100, res = 0.01, ...) +{ + UseMethod("one_box") +} + +#' @rdname one_box +#' @export +one_box.numeric <- function(x, + t_end = 100, res = 0.01, ...) +{ + half_life = x + k = log(2)/half_life + t_out <- seq(0, t_end, by = res) + raw <- matrix(exp( - k * t_out), ncol = 1) + dimnames(raw) <- list(NULL, "parent") + result <- ts(raw, 0, t_end, frequency = 1/res) + class(result) <- c("one_box", "ts") + return(result) +} + +#' @rdname one_box +#' @importFrom mkin mkinpredict +#' @export +one_box.mkinfit <- function(x, t_end = 100, res = 0.01, ...) { + fit <- x + t_out = seq(0, t_end, by = res) + odeini <- c(1, rep(0, length(fit$mkinmod$spec) - 1)) + names(odeini) <- names(fit$mkinmod$spec) + if (length(fit$mkinmod$spec) == 1) solution_type = "analytical" + else solution_type = "deSolve" + + tmp <- mkinpredict(fit$mkinmod, odeparms = fit$bparms.ode, odeini = odeini, + outtimes = t_out, solution_type = solution_type)[-1] + result <- ts(tmp, 0, t_end, frequency = 1/res) + class(result) <- c("one_box", "ts") + return(result) +} + +#' Plot time series of decline data +#' +#' @param x The object of type \code{\link{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 +#' @param ylab Label for the y axis +#' @param max_twa If a numeric value is given, the maximum time weighted +#' average concentration(s) is/are shown in the graph. +#' @param max_twa_var Variable for which the maximum time weighted average should +#' be shown if max_twa is not NULL. +#' @param ... Further arguments passed to methods +#' @importFrom stats plot.ts +#' @export +#' @examples +#' plot(sawtooth(one_box(10), 3, 7), max_twa = 21) +plot.one_box <- function(x, + xlim = range(time(x)), ylim = c(0, max(x)), + xlab = "Time", ylab = "Fraction of initial", + max_twa = NULL, max_twa_var = dimnames(x)[[2]][1], ...) +{ + 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) + if (!is.null(max_twa)) { + x_twa <- max_twa(x, window = max_twa) + value <- x_twa$max[max_twa_var] + rect(x_twa$window_start[max_twa_var], 0, + 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)) + } +} + +#' Create decline time series for multiple applications +#' +#' If the number of application cycles \code{n} is greater than 1, the +#' application pattern specified in \code{applications} is repeated \code{n} +#' times, with an interval \code{i}. +#' @param x A \code{\link{one_box}} object +#' @param n The number of applications. If \code{applications} is specified, \code{n} is not used +#' @param i The interval between applications. If \code{applications} is specified, \code{i} +#' is not used +#' @param applications A data frame holding the application times in the first column and +#' the corresponding amounts applied in the second column for each application cycle. +#' If \code{n} is one, the application pattern specified here is used only once. +#' @export +#' @examples +#' applications = data.frame(time = seq(0, 14, by = 7), amount = c(1, 2, 3)) +#' pred <- one_box(10) +#' plot(sawtooth(pred, applications = applications)) +#' +#' m_2 <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO")) +#' fit_2 <- mkinfit(m_2, FOCUS_2006_D) +#' pred_2 <- one_box(fit_2) +#' pred_2_saw <- sawtooth(pred_2, 2, 7) +#' plot(pred_2_saw, max_twa = 21, max_twa_var = "m1") +#' +#' max_twa(pred_2_saw) +sawtooth <- function(x, n = 1, i = 365, + applications = data.frame(time = seq(0, 0 + n * i, length.out = n), + amount = 1)) +{ + n_obs = ncol(as.matrix(x)) + t_end = max(time(x)) + freq = frequency(x) + empty <- ts(matrix(0, nrow = t_end * freq, ncol = n_obs), 0, t_end, freq) + result <- empty + for (i_app in 1:nrow(applications)) { + t_app <- applications[i_app, "time"] + amount_app <- applications[i_app, "amount"] + if (t_app == 0) { + result <- result + x * amount_app + } else { + lag_phase <- as.matrix(empty)[1:(t_app * freq), , drop = FALSE] + app_phase <- amount_app * as.matrix(x)[1:((t_end - t_app) * freq + 1), , drop = FALSE] + app_ts <- ts(rbind(lag_phase, app_phase), 0, t_end, frequency = freq) + result <- result + app_ts + } + } + class(result) = c("one_box", "ts") + dimnames(result) <- dimnames(x) + return(result) +} + +#' Calculate a time weighted average concentration +#' +#' The moving average is built only using the values in the past, so +#' 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 window The size of the moving window +#' @seealso max_twa +#' @export +#' @examples +#' pred <- sawtooth(one_box(10), +#' applications = data.frame(time = c(0, 7), amount = c(1, 1))) +#' max_twa(pred) +twa <- function(x, window = 21) UseMethod("twa") + +#' @rdname twa +#' @export +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) + class(result) = c("one_box", "ts") + dimnames(result) <- dimnames(x) + return(result) +} + +#' The maximum time weighted average concentration for a moving window +#' +#' @seealso twa +#' @inheritParams twa +#' @export +max_twa <- function(x, window = 21) UseMethod("max_twa") + +#' @export +max_twa.one_box <- function(x, window = 21) +{ + freq = frequency(x) + + twa_ts <- twa(x, window = window) + window_end <- apply(twa_ts, 2, which.max) / freq + result <- list() + result$max <- apply(twa_ts, 2, max, na.rm = TRUE) + result$window_start <- window_end - window + result$window_end <- window_end + return(result) +} diff --git a/man/max_twa.Rd b/man/max_twa.Rd new file mode 100644 index 0000000..13e7bbf --- /dev/null +++ b/man/max_twa.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/twa.R +\name{max_twa} +\alias{max_twa} +\title{The maximum time weighted average concentration for a moving window} +\usage{ +max_twa(x, window = 21) +} +\arguments{ +\item{x}{An object of type \code{\link{one_box}}} + +\item{window}{The size of the moving window} +} +\description{ +The maximum time weighted average concentration for a moving window +} +\seealso{ +twa +} diff --git a/man/one_box.Rd b/man/one_box.Rd new file mode 100644 index 0000000..86440ce --- /dev/null +++ b/man/one_box.Rd @@ -0,0 +1,45 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/twa.R +\name{one_box} +\alias{one_box} +\alias{one_box.numeric} +\alias{one_box.mkinfit} +\title{Create a time series of decline data} +\usage{ +one_box(x, t_end = 100, res = 0.01, ...) + +\method{one_box}{numeric}(x, t_end = 100, res = 0.01, ...) + +\method{one_box}{mkinfit}(x, t_end = 100, res = 0.01, ...) +} +\arguments{ +\item{x}{When numeric, this is the half-life to be used for an exponential +decline. If x is an mkinfit object, the decline is calculated from this object} + +\item{t_end}{End of the time series} + +\item{res}{Resolution of the time series} + +\item{...}{Further arguments passed to methods} +} +\description{ +The time series starts with the amount specified for the first application. +This does not create objects of type \code{\link{ts}}. +} +\examples{ +# Only use a half-life +pred_0 <- one_box(10) +plot(pred_0) + +# Use a fitted mkinfit model +require(mkin) +fit <- mkinfit("FOMC", FOCUS_2006_C) +pred_1 <- one_box(fit) +plot(pred_1) + +# Use a model with more than one observed variable +m_2 <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO")) +fit_2 <- mkinfit(m_2, FOCUS_2006_D) +pred_2 <- one_box(fit_2) +plot(pred_2) +} diff --git a/man/plot.one_box.Rd b/man/plot.one_box.Rd new file mode 100644 index 0000000..7d30bd2 --- /dev/null +++ b/man/plot.one_box.Rd @@ -0,0 +1,35 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/twa.R +\name{plot.one_box} +\alias{plot.one_box} +\title{Plot time series of decline data} +\usage{ +\method{plot}{one_box}(x, xlim = range(time(x)), ylim = c(0, max(x)), + xlab = "Time", ylab = "Fraction of initial", max_twa = NULL, + max_twa_var = dimnames(x)[[2]][1], ...) +} +\arguments{ +\item{x}{The object of type \code{\link{one_box}} to be plotted} + +\item{xlim}{Limits for the x axis} + +\item{ylim}{Limits for the y axis} + +\item{xlab}{Label for the x axis} + +\item{ylab}{Label for the y axis} + +\item{max_twa}{If a numeric value is given, the maximum time weighted +average concentration(s) is/are shown in the graph.} + +\item{max_twa_var}{Variable for which the maximum time weighted average should +be shown if max_twa is not NULL.} + +\item{...}{Further arguments passed to methods} +} +\description{ +Plot time series of decline data +} +\examples{ +plot(sawtooth(one_box(10), 3, 7), max_twa = 21) +} diff --git a/man/sawtooth.Rd b/man/sawtooth.Rd new file mode 100644 index 0000000..4576c8b --- /dev/null +++ b/man/sawtooth.Rd @@ -0,0 +1,39 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/twa.R +\name{sawtooth} +\alias{sawtooth} +\title{Create decline time series for multiple applications} +\usage{ +sawtooth(x, n = 1, i = 365, applications = data.frame(time = seq(0, 0 + n + * i, length.out = n), amount = 1)) +} +\arguments{ +\item{x}{A \code{\link{one_box}} object} + +\item{n}{The number of applications. If \code{applications} is specified, \code{n} is not used} + +\item{i}{The interval between applications. If \code{applications} is specified, \code{i} +is not used} + +\item{applications}{A data frame holding the application times in the first column and +the corresponding amounts applied in the second column for each application cycle. +If \code{n} is one, the application pattern specified here is used only once.} +} +\description{ +If the number of application cycles \code{n} is greater than 1, the +application pattern specified in \code{applications} is repeated \code{n} +times, with an interval \code{i}. +} +\examples{ +applications = data.frame(time = seq(0, 14, by = 7), amount = c(1, 2, 3)) +pred <- one_box(10) +plot(sawtooth(pred, applications = applications)) + +m_2 <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO")) +fit_2 <- mkinfit(m_2, FOCUS_2006_D) +pred_2 <- one_box(fit_2) +pred_2_saw <- sawtooth(pred_2, 2, 7) +plot(pred_2_saw, max_twa = 21, max_twa_var = "m1") + +max_twa(pred_2_saw) +} diff --git a/man/twa.Rd b/man/twa.Rd new file mode 100644 index 0000000..a6afa87 --- /dev/null +++ b/man/twa.Rd @@ -0,0 +1,29 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/twa.R +\name{twa} +\alias{twa} +\alias{twa.one_box} +\title{Calculate a time weighted average concentration} +\usage{ +twa(x, window = 21) + +\method{twa}{one_box}(x, window = 21) +} +\arguments{ +\item{x}{An object of type \code{\link{one_box}}} + +\item{window}{The size of the moving window} +} +\description{ +The moving average is built only using the values in the past, so +the earliest possible time for the maximum in the time series returned +is after one window has passed. +} +\examples{ +pred <- sawtooth(one_box(10), + applications = data.frame(time = c(0, 7), amount = c(1, 1))) +max_twa(pred) +} +\seealso{ +max_twa +} |