aboutsummaryrefslogtreecommitdiff
path: root/R/PELMO_runs.R
blob: c0d126668e7772ae85ee299fd6115d97c6adfb09 (plain) (blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
#' 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{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 matrices 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 existing run directories be overwritten?
#' @references PELMO.installeR \url{https://pkgdown.jrwb.de/PELMO.installeR}
#'
#' Wine \url{https://winehq.org}
#'
#' PELMO test results \url{http://esdac.jrc.ec.europa.eu/public_path/projects_data/focus/gw/models/pelmo/test-results/test_results_FOCUS_PELMO_5_5_3.doc}
#' @export
#' @examples
#' # At the moment I can not run the examples, as my wine installation is not working
#' \dontrun{
#' # Reproduce the official test results for annual application of Pesticide D
#' # to winter cereals at the day before emergence
#' runs_1 <- list(
#'   list(psm = 'Pesticide_D',
#'        win = c("Cha", "Ham", "Jok", "Kre", "Oke", "Pia", "Por", "Sev", "Thi")),
#'   list(psm = 'Pesticide_D_1_day_pre_em_every_third_year',
#'        pot = c("Cha", "Ham", "Jok", "Kre", "Oke", "Pia", "Por", "Sev", "Thi")))
#' time_1 <- system.time(
#'   PECgw_1 <- PELMO_runs(runs_1, psm_dir = system.file("testdata", package = "pfm"),
#'     cores = 6, overwrite = TRUE)
#' )
#' print(PECgw_1)
#' # We get exactly the same PECgw values (on Linux, calling PELMO using Wine).
#' print(time_1)
#'
#' # Demonstrate some results with metabolites.
#' runs_2 <- list(list(psm = 'Pesticide_D_1_May_every_other_year_mets',
#'                     win = c("Cha", "Ham", "Kre")))
#' PECgw_2 <- PELMO_runs(runs_2, psm_dir = system.file("testdata", package = "pfm"),
#'   cores = 3, overwrite = TRUE)
#' print(PECgw_2)
#' }
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)
  }

  lapply(run_list, setup_run)

  if (execute) {
    run_PELMO(runs, version = version, PELMO_base = PELMO_base, cores = cores)
  }

  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
#'
#' @param psm_file The path to the .psm file
#' @param location_code The location code
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
#' 
#' @param annual The annual flux as obtained by \code{get_flux}
#' @param interval The interval in years
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
#' 
#' @param chem_file The full path to 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