diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/FOCUS_PELMO_data.R | 14 | ||||
| -rw-r--r-- | R/PELMO_runs.R | 163 |
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)])) +} |