Setting up PELMO runs, execution and evaluation

It all works!

2 files changed, 148 insertions, 29 deletions

diff --git a/R/FOCUS_PELMO_data.R b/R/FOCUS_PELMO_data.R
index 28888e7..89fe5f6 100644
--- a/R/FOCUS_PELMO_data.R
+++ b/R/FOCUS_PELMO_data.R
@@ -2,10 +2,10 @@
 #'
 #' A named character vector with the crop names used in the PELMO 5.5.3 GUI.
 #' For the names, three letter codes were constructed by generally taking the
-#' first three letters in lower case. Only when there is an expression in 
-#' parentheses, the first letter of this expression becomes the first letter 
+#' first three letters in lower case. Only when there is an expression in
+#' parentheses, the first letter of this expression becomes the first letter
 #' in the three letter code, i.e. 'Peas (animals)' has the code \code{ape}.
-#' 
+#'
 #' @docType data
 #' @export
 #' @examples
@@ -37,13 +37,13 @@ FOCUS_PELMO_crops = c(
   tom = "Tomatoes",
   vin = "Vines")
 
-#' FOCUS PELMO crop acronyms used for naming sceneario files 
+#' FOCUS PELMO crop acronyms used for naming sceneario files
 #'
 #' A named character vector with the crop acronyms used in FOCUS PELMO 5.5.3
 #' for naming the .sze files located in the FOCUS directory. The crop acronyms
 #' in the file names are sometimes capitalized, sometimes not. The scenario
 #' files used for Beans (field) and Beans (vegetable) are the same.
-#' 
+#'
 #' @docType data
 #' @export
 #' @examples
@@ -85,6 +85,6 @@ FOCUS_PELMO_crop_sze_names = c(
 #' @examples
 #' print(FOCUS_PELMO_location_codes)
 FOCUS_PELMO_location_codes = c(
-  Cha = "C", Ham = "H", Jok = "J", 
-  Kre = "K", Oke = "N", Pia = "P", 
+  Cha = "C", Ham = "H", Jok = "J",
+  Kre = "K", Oke = "N", Pia = "P",
   Por = "O", Sev = "S", Thi = "T")
diff --git a/R/PELMO_runs.R b/R/PELMO_runs.R
index f599110..5d401ea 100644
--- a/R/PELMO_runs.R
+++ b/R/PELMO_runs.R
@@ -1,7 +1,10 @@
 #' Set up runs for FOCUS PELMO
 #'
 #' Per default, the runs are also executed with FOCUS PELMO, and the results are processed
-#' and returned.
+#' and returned. Currently, only FOCUS PELMO as installed on Linux (or other Unix systems)
+#' using the \code{\link{install_PELMO}} from the \code{PELMO.installeR} package
+#' maintained on github is supported. In such installations, FOCUS PELMO is installed into
+#' package installation directory of \code{PELMO.installeR} and run with \code{wine}.
 #'
 #' @param runs A list of lists. Each inner lists has an element named 'psm'
 #'   that holds the psm string, and elements named using three letter crop acronyms,
@@ -15,6 +18,9 @@
 #' @param cores The number of cores to execute PELMO runs in parallel
 #' @param evaluate Should the results be returned?
 #' @param overwrite Should an existing run directories be overwritten
+#' @references PELMO.installeR \url{https://jranke.github.io/PELMO.installeR}
+#' 
+#' Wine \url{https://winehq.org}
 #' @export
 PELMO_runs <- function(runs, psm_dir = ".", version = "5.5.3", PELMO_base = "auto",
                        execute = TRUE, cores = getOption("mc.cores", 2L),
@@ -80,26 +86,9 @@ PELMO_runs <- function(runs, psm_dir = ".", version = "5.5.3", PELMO_base = "aut
 
     add <- function(x) cat(paste0(x, "\r\n"), file = input_file, append = TRUE)
 
-    # How many years do we calculate (26, 46 or 66)?
-    psm <- readLines(psm_file, encoding = "latin1")
-    number_of_apps_lines <- grep("number of application location", psm)
-    absolute_apps_line <- grep(location_code,
-                               psm[number_of_apps_lines])
-    period <- 1 # application every year
-    if (length(absolute_apps_line) == 1) {
-      apps_root <- number_of_apps_lines[absolute_apps_line]
-    } else {
-      apps_root <- number_of_apps_lines[1]
-    }
-
-    number_of_apps <- as.integer(substr(psm[apps_root], 1, 3))
-    last_app_line <- psm[apps_root + number_of_apps]
-    last_app_year <- as.integer(gsub("^.{2,3}  ..  (..)   .*", "\\1",
-                                     last_app_line))
-    if (last_app_year > 26) period <- 2
-    if (last_app_year > 46) period <- 3
+    interval <- get_interval(psm_file, location_code)
 
-    n_years <- switch(as.character(period),
+    n_years <- switch(as.character(interval),
       "1" = 26,
       "2" = 46,
       "3" = 66)
@@ -145,18 +134,23 @@ PELMO_runs <- function(runs, psm_dir = ".", version = "5.5.3", PELMO_base = "aut
     # the exe file from this directory
     file.copy(file.path(PELMO_base, "lf90.eer"), run_dir)
   }
+  run_list[[7]]
+  setup_run(run_list[[7]])
 
   lapply(run_list, setup_run)
 
   if (execute) {
     run_PELMO(runs, version = version, PELMO_base = PELMO_base)
   }
+
+  if (evaluate) {
+    evaluate_PELMO(runs, version = version, PELMO_base = PELMO_base)
+  }
 }
 
-#' Run PELMO
-#'
 #' @inheritParams PELMO_runs
 #' @importFrom parallel mclapply
+#' @rdname PELMO_runs
 #' @export
 run_PELMO <- function(runs, version = "5.5.3", PELMO_base = "auto",
                       cores = getOption("mc.cores", 2L))
@@ -295,3 +289,128 @@ create_run_list <- function(runs, psm_dir = ".", check_psm_files = FALSE) {
   }
   return(run_list)
 }
+
+#' @rdname PELMO_runs
+#' @inheritParams PELMO_runs
+#' @importFrom parallel mclapply
+#' @export
+evaluate_PELMO <- function(runs, version = "5.5.3", PELMO_base = "auto")
+{
+
+  if (PELMO_base[1] == "auto") {
+    PELMO_base = file.path(system.file(package = "PELMO.installeR"),
+                                       paste0("FOCUSPELMO.", gsub("\\.", "", version)))
+  }
+
+  for (run in runs) {
+    psm <- run$psm
+    crops <- setdiff(names(run), "psm")
+    for (crop in crops) {
+      for (scenario in run[[crop]]) {
+        run_dir <- file.path(PELMO_base, "FOCUS", PELMO_path(psm, crop, scenario))
+        echo_file <- readLines(file.path(run_dir, "ECHO.PLM"), encoding = "latin1")
+        parm_lines <- grep("\\*\\*\\* PARAMETERS OF", echo_file, value = TRUE)
+        acronyms <- gsub(".*\\((.*)\\).*", "\\1", parm_lines)
+        met_codes <- gsub(".*METABOLITE (..).*", "\\1", parm_lines)
+        met_codes[1] <- NA
+        names(met_codes) <- acronyms
+
+        psm_file <- file.path(run_dir, paste0(psm, ".psm"))
+        location_code <- FOCUS_PELMO_location_codes[scenario]
+        interval <- get_interval(psm_file, location_code)
+
+        # Get percolate for each period
+        wasser <- readLines(file.path(run_dir, "WASSER.PLM"))
+        percolate_lines <- grep("RECHARGE BELOW ROOT ZONE", wasser, value = TRUE)
+        percolate <- as.numeric(substr(percolate_lines, 40, 46))
+        percolate_period <- sum_periods(percolate, interval)
+
+        # Set up results that should match period.plm generated by the PELMO GUI
+        results_pfm <- list()
+        for (acronym in acronyms) {
+          results_pfm[[acronym]] <- list()
+          periods <- data.frame(
+            period = 1:20,
+            flux = NA,
+            percolate = 10 * percolate_period,
+            conc = NA)
+
+          if (is.na(met_codes[acronym])) {
+            chem_file <- file.path(run_dir, "CHEM.PLM")
+          } else {
+            chem_file <- file.path(run_dir, paste0("CHEM_", met_codes[acronym], ".PLM"))
+          }
+
+          annual_flux <- get_flux(chem_file)
+          periods$flux <- sum_periods(annual_flux, interval) * 1000
+
+          periods$conc <- 100 * periods$flux / periods$percolate
+
+          results_pfm[[acronym]]$periods <- periods
+          results_pfm[[acronym]]$focus <- focus_80th(periods$conc)
+        }
+        save(results_pfm, file = file.path(run_dir, "period_pfm.rda"))
+      }
+    }
+  }
+}
+
+#' Get the application interval in years from a psm file
+get_interval <- function(psm_file, location_code) {
+  # How many years do we calculate (26, 46 or 66)?
+  psm <- readLines(psm_file, encoding = "latin1")
+  number_of_apps_lines <- grep("number of application location", psm)
+  absolute_apps_line <- grep(location_code,
+                             psm[number_of_apps_lines])
+  interval <- 1 # application every year
+  if (length(absolute_apps_line) == 1) {
+    apps_root <- number_of_apps_lines[absolute_apps_line]
+  } else {
+    apps_root <- number_of_apps_lines[1]
+  }
+
+  number_of_apps <- as.integer(substr(psm[apps_root], 1, 3))
+  last_app_line <- psm[apps_root + number_of_apps]
+  last_app_year <- as.integer(gsub("^.{2,3}  ..  (..)   .*", "\\1",
+                                   last_app_line))
+  if (last_app_year > 26) interval <- 2
+  if (last_app_year > 46) interval <- 3
+  return(interval)
+}
+
+#' Sum up values according to FOCUS periods
+sum_periods <- function(annual, interval) {
+  n_years <- switch(as.character(interval),
+    "1" = 26,
+    "2" = 46,
+    "3" = 66)
+  period_start_years <- seq(from = 7, to = n_years, by = interval)
+  sapply(1:20, function(x) {
+    years_i <- period_start_years[x] + (0 : (interval - 1))
+    sum(annual[years_i])
+  })
+}
+
+#' Get the flux of a chemical out of the FOCUS layer from a CHEM*.PLM file
+get_flux <- function(chem_file) {
+  chem <- readLines(chem_file)
+  lowest_focus_comp_lines <- grep("^    .     21      ", chem, value = TRUE)
+  lowest_focus_comp <- read.table(text = lowest_focus_comp_lines)
+  return(lowest_focus_comp$V9)
+}
+
+#' Calculate the 80th percentile according to FOCUS guidance
+#'
+#' This is nowadays defined as the mean of the 16th and the 17th
+#' highest value. Previously, the 17th highest values was used (FOCUS 2014, p.
+#' 18). NaN values need to be set to zero in order to reproduce the
+#' values obtained by PELMO.
+#' @param c_period A numeric vector of values to calculate the percentile from
+#' @param old Should the old calculation method be used (the 17th highest value)?
+#' @export
+focus_80th <- function(c_period, old = FALSE) {
+  c_period <- ifelse(is.na(c_period), 0, c_period)
+   c_period_sorted <- sort(c_period)
+  if (old) return(c_period_sorted[17])
+  else return(mean(c_period_sorted[c(16, 17)]))
+}