From 92d97ba99d7b90d95a67796cb5e68f28f752b70b Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 17 Dec 2014 21:39:42 +0100 Subject: Initial commit: R6 class for TOXSWA cwa files --- pkg/R/TOXSWA_cwa.R | 267 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 267 insertions(+) create mode 100644 pkg/R/TOXSWA_cwa.R (limited to 'pkg/R/TOXSWA_cwa.R') diff --git a/pkg/R/TOXSWA_cwa.R b/pkg/R/TOXSWA_cwa.R new file mode 100644 index 0000000..113dc42 --- /dev/null +++ b/pkg/R/TOXSWA_cwa.R @@ -0,0 +1,267 @@ +# Copyright (C) 2014 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. +#' +#' @param filename The filename of the cwa file. +#' @param basedir The full path to the directory where the cwa file resides. +#' @param segment The segment for which the data should be read. Either "last", or +#' the segment number. +#' @param windows Numeric vector of width of moving windows for calculating +#' time weighted average concentrations and areas under the curve. +#' @param thresholds Numeric vector of threshold concentrations in µg/L for +#' generating event statistics. +#' @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", +#' system.file("testdata/SwashProjects/project_H_sw/TOXSWA", +#' package = "pfm")) +read.TOXSWA_cwa <- function(filename, basedir = ".", segment = "last", + windows = c(7, 21), thresholds = NULL) +{ + if (!missing(filename)) { + cwa <- TOXSWA_cwa$new(filename, basedir) + 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 segment of a TOXSWA surface water body. +#' +#' @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", +#' system.file("testdata/SwashProjects/project_H_sw/TOXSWA", +#' package = "pfm")) +#' plot(H_sw_D4_pond) +plot.TOXSWA_cwa <- function(x, time_column = c("datetime", "t", "t_firstjan"), + 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]") + } + 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{ + 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", +#' system.file("testdata/SwashProjects/project_H_sw/TOXSWA", +#' package = "pfm")) +#' H_sw_R1_stream$get_events(c(2, 10)) +#' @keywords data + +TOXSWA_cwa <- R6Class("TOXSWA_cwa", + public = list( + filename = "character", + basedir = "character", + segment = "integer", + cwas = "data.frame", + windows = "matrix", + events = "list", + initialize = function(filename, basedir, segment = "last") { + self$filename <- filename + self$basedir <- basedir + cwa_all_segments <- try(read.table(file.path(basedir, filename), + sep = "", skip = 40, + colClasses = c("character", "numeric", + "integer", rep("numeric", 5)), + col.names = c("datetime", "t", "segment", + "xcd", "cwa_tot", "cwa", + "Xss", "Xmp"))) + 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")] + self$events <- list() + } else { + stop("Could not read ", filename) + } + }, + 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 + } + } 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: -- cgit v1.2.1