aboutsummaryrefslogtreecommitdiff
path: root/R/PELMO_runs.R
diff options
context:
space:
mode:
Diffstat (limited to 'R/PELMO_runs.R')
-rw-r--r--R/PELMO_runs.R163
1 files changed, 141 insertions, 22 deletions
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)]))
+}

Contact - Imprint