Read cwa data from TOXSWA 4.4.2 .out files

2 files changed, 93 insertions, 40 deletions

diff --git a/pkg/DESCRIPTION b/pkg/DESCRIPTION
index 7e6d064..1cc6817 100644
--- a/pkg/DESCRIPTION
+++ b/pkg/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: pfm
 Type: Package
 Title: Utilities for Pesticide Fate Modelling
-Version: 0.3-1
-Date: 2015-12-27
+Version: 0.3-2
+Date: 2016-02-04
 Authors@R: person("Johannes Ranke", email = "jranke@uni-bremen.de", 
                   role = c("aut", "cre", "cph"))
 Description: Utilities for simple calculations of predicted environmental
diff --git a/pkg/R/TOXSWA_cwa.R b/pkg/R/TOXSWA_cwa.R
index df7f627..46b1995 100644
--- a/pkg/R/TOXSWA_cwa.R
+++ b/pkg/R/TOXSWA_cwa.R
@@ -21,11 +21,15 @@
 #' 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 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 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
@@ -40,11 +44,11 @@
 #'                                  zipfile = system.file("testdata/SwashProjects.zip",
 #'                                                        package = "pfm"))
 read.TOXSWA_cwa <- function(filename, basedir = ".", zipfile = NULL, 
-                            segment = "last",
+                            segment = "last", total = FALSE,
                             windows = NULL, thresholds = NULL)
 {
   if (!missing(filename)) {
-    cwa <- TOXSWA_cwa$new(filename, basedir, zipfile)
+    cwa <- TOXSWA_cwa$new(filename, basedir, zipfile, total = total)
     if (!is.null(windows[1])) cwa$moving_windows(windows)
     if (!is.null(thresholds[1])) cwa$get_events(thresholds)
     invisible(cwa)
@@ -160,8 +164,8 @@ TOXSWA_cwa <- R6Class("TOXSWA_cwa",
     segment = NULL,
     cwas = NULL,
     windows = NULL,
-    events = NULL,
-    initialize = function(filename, basedir, zipfile = NULL, segment = "last") {
+    events = list(),
+    initialize = function(filename, basedir, zipfile = NULL, segment = "last", total = FALSE) {
       self$filename <- filename
       self$basedir <- basedir
       self$zipfile <- zipfile
@@ -170,42 +174,91 @@ TOXSWA_cwa <- R6Class("TOXSWA_cwa",
       } else {
         try(file_connection <- file(file.path(basedir, filename), "rt"))
       }
-      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"))
+      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
 
-        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()
+        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 {
-        stop("Could not read ", filename)
+        # out file from FOCUS TOXSWA 4 (TOXSWA 4.4.2 or similar)
+        outfile <- try(readLines(file_connection))
+
+        if (is.null(zipfile)) {
+          close(file_connection) # only needed for files
+        } 
+
+        if (inherits(outfile, "try-error")) {
+          stop("Could not read ", filename)
+        } else {
+          cwa_lines <- outfile[grep("ConLiqWatLay_", outfile)] # hourly concentrations
+          cwa_all_segments <- read.table(text = cwa_lines) 
+
+          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$V2),
+            t = cwa_all_segments$V1,
+            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) {