#' Set up runs for FOCUS PELMO
#'
#' Per default, the runs are not only set up but also executed with FOCUS
#' PELMO, the results are processed 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 the package installation directory of \code{PELMO.installeR}
#' and run using \code{wine}.
#'
#' As a side effect, an R data file (period_pfm.rda) is generated in each
#' run directory, holding the results for all FOCUS periods, equivalent to
#' the period.plm file generated by the FOCUS PELMO GUI.
#'
#' @return If evaluate is TRUE, a list of lists of data frames holding the
#' PEC data.
#' @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,
#' as used in \code{\link{FOCUS_PELMO_crops}},
#' that hold character vectors of three letter scenario acronyms
#' as used in \code{\link{FOCUS_GW_scenarios_2012}}.
#' @param psm_dir The directory where the psm files are located
#' @param version The FOCUS PELMO version
#' @param PELMO_base Where the FOCUS PELMO installation is located
#' @param execute Should PELMO be executed directly?
#' @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),
evaluate = TRUE, overwrite = FALSE)
{
if (PELMO_base[1] == "auto") {
PELMO_base = file.path(system.file(package = "PELMO.installeR"),
paste0("FOCUSPELMO.", gsub("\\.", "", version)))
}
PELMO_exe = list.files(PELMO_base, "^pelmo.*.exe", full.names = TRUE)
if (length(PELMO_exe) == 0) {
stop("Could not find PELMO executable. Did you run PELMO.installeR::install_PELMO()?")
}
run_list <- create_run_list(runs, psm_dir = psm_dir, check_psm_files = TRUE)
setup_run <- function(x) {
psm <- x$psm
crop <- x$crop
scenario <- x$scenario
location_code <- FOCUS_PELMO_location_codes[scenario]
run_dir <- file.path(PELMO_base, "FOCUS", PELMO_path(psm, crop, scenario))
if (dir.exists(run_dir)) {
if (overwrite) {
unlink(run_dir, recursive = TRUE)
} else {
stop("Run directory for ", crop, " in ", scenario, "\n", run_dir, "\nalready exists")
}
}
# Create the run directory
dir.create(run_dir, recursive = TRUE)
# Copy the psm file
psm_file <- file.path(psm_dir, paste0(psm, ".psm"))
file.copy(psm_file, run_dir)
# Copy Haude factors
file.copy(file.path(PELMO_base, "FOCUS", "HAUDE.DAT"), file.path(run_dir, "Haude.dat"))
# Copy scenario file
sze_file_upper <- paste0(location_code, "_",
toupper(FOCUS_PELMO_crop_sze_names[crop]), ".sze")
sze_file_lower <- paste0(location_code, "_",
FOCUS_PELMO_crop_sze_names[crop], ".sze")
if (file.exists(file.path(PELMO_base, "FOCUS", sze_file_lower))) {
sze_file <- sze_file_lower
} else {
if (file.exists(file.path(PELMO_base, "FOCUS", sze_file_upper))) {
sze_file <- sze_file_upper
} else {
stop("Could not find szenario file for ", crop, " in ", scenario)
}
}
file.copy(file.path(PELMO_base, "FOCUS", sze_file), file.path(run_dir, sze_file_lower))
# Generate PELMO input file
input_file <- file(file.path(run_dir, "pelmo.inp"), encoding = "latin1", open = "w+")
on.exit(close(input_file))
add <- function(x) cat(paste0(x, "\r\n"), file = input_file, append = TRUE)
interval <- get_interval(psm_file, location_code)
n_years <- switch(as.character(interval),
"1" = 26,
"2" = 46,
"3" = 66)
# First line with number of climate years
add(paste0(formatC(n_years, width = 3, flag = "0"),
" 01 01 31 12 Version ", substr(version, 1, 1)))
# Second line with scenario file (why did we copy it, if this line refers to the FOCUS dir...
add(paste0("..\\..\\..\\", sze_file_lower))
# Third line with psm file
add(basename(psm_file))
# Lines with climate files
for (year in 1:n_years) {
add(paste0("..\\..\\..\\", location_code, "_", formatC(year, width = 2, flag = "0"), ".cli"))
}
# Output control section
add("00000015 00.")
add(" PRSN TSER 000.0 1.0 00000001")
add(" TETD TSER 000.0 1.0 00000001")
add(" INFL TSER 100.0 1.0 00000001")
add(" RUNF TSER 000.0 1.0 00000001")
add(" THET TSER 000.0 1.0 00000001")
add(" THET TSER 030.0 1.0 00000001")
add(" TEMP TSER 000.0 1.0 00000001")
add(" TEMP TSER 030.0 1.0 00000001")
add(" TPAP TSER 000.0 1.0E05 00000001")
add(" TDKF TSER 000.0 1.0E05 00000001")
add(" TUPF TSER 000.0 1.0E05 00000001")
add(" TPST TSER 005.0 1.0E06 00000001")
add(" PFLX TSER 100.0 1.0E05 00000001")
add(" RFLX TSER 000.0 1.0E05 00000001")
add(" LEAC TSER 100.0 1.0E09 00000001")
# Copy pelmo executable
file.copy(PELMO_exe, run_dir)
# In addition to the files copied by the FOCUS PELMO GUI, we
# need the error file in the run directory, as we start
# 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) {
pfm_PECgw <- evaluate_PELMO(runs, version = version, PELMO_base = PELMO_base)
return(pfm_PECgw)
} else {
invisible(NULL)
}
}
#' @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))
{
if (PELMO_base[1] == "auto") {
PELMO_base = file.path(system.file(package = "PELMO.installeR"),
paste0("FOCUSPELMO.", gsub("\\.", "", version)))
}
pelmo_exe = list.files(PELMO_base, "^pelmo.*.exe")
# Create list of runs to traverse using mclapply
run_list <- create_run_list(runs, check_psm_files = FALSE)
# In order to run in parallel, we need to make a directory tree starting
# from the base of the PELMO installation, as pelmo401.exe seems to look
# for ..\..\..\summary.PLM and does not start if this is in use by another
# instance
run_exe <- function(x) {
psm <- x$psm
crop <- x$crop
scenario <- x$scenario
run_dir <- file.path(PELMO_base, "FOCUS", PELMO_path(psm, crop, scenario))
exe_dir <- file.path(PELMO_base, paste(sample(c(0:9, letters, LETTERS),
size = 8, replace = TRUE),
collapse = ""))
dir.create(exe_dir)
# Fresh PELMO base directory for every run with random name
run_dir_exe <- file.path(exe_dir, "FOCUS", PELMO_path(psm, crop, scenario))
dir.create(run_dir_exe, recursive = TRUE)
# Copy the contents of the run directory
run_files <- list.files(run_dir, full.names = TRUE)
file.copy(run_files, run_dir_exe)
# Copy FOCUS files
location_code <- FOCUS_PELMO_location_codes[x$scenario]
location_files <- list.files(file.path(PELMO_base, "FOCUS"), paste0("^", location_code),
full.names = TRUE)
focus_dir_exe <- file.path(exe_dir, "FOCUS")
file.copy(location_files, focus_dir_exe)
file.copy(file.path(PELMO_base, "FOCUS", "FOCUSCHECK.DAT"), focus_dir_exe)
# We need to go the directory to simplify calling pelmo with wine
setwd(run_dir_exe)
psm_file <- paste0(psm, ".psm")
message("Starting ", pelmo_exe, " with ", psm_file)
system(paste("wine", pelmo_exe, psm_file), ignore.stdout = TRUE)
# Copy the results to the original run directory
plm_files <- list.files(run_dir_exe, ".plm$", full.names = TRUE)
PLM_files <- list.files(run_dir_exe, ".PLM$", full.names = TRUE)
file.copy(c(plm_files, PLM_files), run_dir)
# Clean up
unlink(exe_dir, recursive = TRUE)
}
mclapply(run_list, run_exe, mc.cores = cores)
}
#' Create a path of run directories as the PELMO GUI does
#'
#' @export
#' @importFrom utils data
#' @param psm The psm identifier
#' @param crop The PELMO crop acronym
#' @param scenario The scenario
PELMO_path <- function(psm, crop, scenario = NA) {
if (crop %in% names(FOCUS_PELMO_crops)) {
crop <- FOCUS_PELMO_crops[crop]
}
if (!crop %in% FOCUS_PELMO_crops) stop(crop, " is not in FOCUS_PELMO_crops.")
psm_dir <- paste0(psm, ".run")
crop_dir <- paste0(gsub(" ", "_-_", crop), ".run")
# Deal with 'irrigation'. It only affects the naming of the scenario directory,
# but it does not appear to affect the PELMO run. This can be seen when
# comparing runs set up for Beans (field) and Beans (vegetable) for the Porto
# scenario
irrigation_string <- ""
if (is.na(scenario)) {
return(file.path(psm_dir, crop_dir))
} else {
# 'Irrigation' is only possible in the GUI for five scenarios
if (scenario %in% c("Cha", "Pia", "Por", "Sev", "Thi")) {
crop_acronyms <- names(FOCUS_PELMO_crops)
names(crop_acronyms) <- FOCUS_PELMO_crops
crop_acronym <- crop_acronyms[[crop]]
# Some crops are not 'irrigated' according to the GUI
if (crop_acronym %in% c("win", "fbe", "woi", "ape", "spr")) {
irrigation_string <- "--_-_no_-_irrigation"
} else {
irrigation_string <- "--_-_irrigated"
}
}
scenario_dir <- paste0(
FOCUS_GW_scenarios_2012$names[scenario], "_-_(", FOCUS_PELMO_location_codes[scenario], ")",
irrigation_string, ".run")
return(file.path(psm_dir, crop_dir, scenario_dir))
}
}
#' Create a list of runs that we can traverse
#'
#' @inheritParams PELMO_runs
#' @param check_psm_files Should we check if the psm file exists
create_run_list <- function(runs, psm_dir = ".", check_psm_files = FALSE) {
i <- 0
run_list <- list()
for (run in runs) {
psm <- run$psm
if (check_psm_files) {
psm_file <- file.path(psm_dir, paste0(psm, ".psm"))
if (file.access(psm_file) == -1) {
stop(psm_file, " is not readable")
}
}
crops <- setdiff(names(run), "psm")
for (crop in crops) {
crop_acronyms <- names(FOCUS_PELMO_crops)
names(crop_acronyms) <- FOCUS_PELMO_crops
if (!crop %in% crop_acronyms) {
if (crop %in% FOCUS_PELMO_crops) {
crop <- crop_acronyms[crop]
} else {
stop("Invalid crop specification ", crop)
}
}
for (scenario in run[[crop]]) {
i <- i + 1
run_list[[i]] <- list(psm = psm, crop = crop, scenario = scenario)
}
}
}
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)))
}
pfm_PECgw <- list()
for (run in runs) {
psm <- run$psm
pfm_PECgw[[psm]] <- list()
crops <- setdiff(names(run), "psm")
# Get acronyms of simulated compounds
example_run_dir <- file.path(PELMO_base, "FOCUS", PELMO_path(psm, crops[1], run[[crops[1]]][1]))
example_echo_file <- readLines(file.path(example_run_dir, "ECHO.PLM"), encoding = "latin1")
parm_lines <- grep("\\*\\*\\* PARAMETERS OF", example_echo_file, value = TRUE)
acronyms <- gsub(".*\\((.*)\\).*", "\\1", parm_lines)
met_codes <- gsub(".*METABOLITE (..).*", "\\1", parm_lines)
met_codes[1] <- NA
names(met_codes) <- acronyms
# Loop over runs to get results
for (crop in crops) {
scenarios <- run[[crop]]
pfm_PECgw[[psm]][[crop]] <- matrix(nrow = length(scenarios), ncol = length(acronyms),
dimnames = list(scenarios, acronyms))
for (scenario in scenarios) {
run_dir <- file.path(PELMO_base, "FOCUS", PELMO_path(psm, crop, scenario))
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
PECgw <- focus_80th(periods$conc)
results_pfm[[acronym]]$focus <- PECgw
pfm_PECgw[[psm]][[crop]][scenario, acronym] <- round(PECgw, 3)
}
save(results_pfm, file = file.path(run_dir, "period_pfm.rda"))
}
}
}
return(pfm_PECgw)
}
#' 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)]))
}
