R/edr.wrapper.R

In ashiklom/PEcAnRTM: PEcAn functions used for radiative transfer modeling

#' ED radiative transfer module (EDR) wrapper function
#'
#' This function provides a convenient way to call the ED
#' radiative transfer module (EDR, which simulates full spectral return of an
#' ED patch for a given point in time) directly from R.
#'
#' @param img_path Path to Singularity container (usually a `.simg` file)
#' @param ed2in_path Path to ED2IN file.
#' @param spectra_list List of spectral data matrices. Names must exactly match
#' the PFTs given in `trait.values`. Each item must be a matrix of
#' reflectance (labeled 'R' or 'reflectance') and transmittance (labeled 'T' or
#' 'transmittance') values. If the matrix is not of class `spectra` (see
#' [spectra()]), it must also have a column of wavelength values (labeled
#' 'wl'). Such a matrix is returned by default by all versions of PROSPECT in
#' this package.
#' @param trait.values Named, hierarchical list of trait values for generating
#' config.xml file. Names must be PFTs, and must exactly match names of
#' `spectra_list`.
#' @param soil_reflect_path Path to soil reflectance file (defaults to spectra
#' in package `extdata`)
#' @param wood_reflect_path Path to wood reflectance file (defaults to spectra
#' in package `extdata`)
#' @param par.wl Vector of wavelengths defining PAR region
#' @param nir.wl Vector of wavelengths defining NIR region
#' @param edr_exe_path If `img_path` is `NULL`, a full path to the EDR
#' executable. Ignored otherwise.
#' @param output.path Directory in which to execute the run. Defaults to
#' `dirname(ed2in_path)`.
#' @param settings PEcAn settings list. Default is `list(model = list(revision
#' = "git", config.header = NULL))`.
#' @param singularity_args Additional arguments to be passed to `singularity run` (before)
#' @param clean Logical. If `TRUE`, remove all files generated by this function
#' (e.g. cloned history file, ED2IN, output HDF files).
#' @param stderr Logical. If `TRUE` (default), internalize `system2` results as
#' R character vector. `TRUE` is recommended because it allows EDR to check its
#' execution and to run more quietly.
#' @param verbose_error Logical. If `TRUE` (default), spit out the full ED
#' if EDR execution fails.
#' @param ... Additional arguments to `system2`
#'
#' @author Alexey Shiklomanov
#' @export
EDR <- function(img_path,
                ed2in_path,
                spectra_list,
                trait.values,
                soil_reflect_path = system.file("extdata", "soil_reflect_par.dat", package = "PEcAnRTM"),
                wood_reflect_path = system.file("extdata", "wood_reflect_par.dat", package = "PEcAnRTM"),
                par.wl = 400:2499,
                nir.wl = 2500,
                edr_exe_path = NULL,
                output.path = dirname(normalizePath(ed2in_path, mustWork = TRUE)),
                settings = list(model = list(revision = "git", config.header = NULL)),
                singularity_args = list(),
                clean = FALSE,
                stderr = TRUE,
                verbose_error = TRUE,
                ...) {

  ed2in_path <- normalizePath(ed2in_path, mustWork = TRUE)

  # Write ED2 config.xml file
  defaults <- list()
  xml <- PEcAn.ED2::write.config.xml.ED2(
    defaults = defaults,
    settings = settings,
    trait.values = trait.values
  )

  new.config.path <- file.path(output.path, "config.xml")
  PREFIX_XML <- '<?xml version="1.0"?>\n<!DOCTYPE config SYSTEM "ed.dtd">\n'
  XML::saveXML(xml, file = new.config.path, indent=TRUE, prefix = PREFIX_XML)

  # Generate input files
  files_list <- file.path(
    output.path,
    c("lengths.dat",
      "reflect_par.dat", "reflect_nir.dat",
      "trans_par.dat", "trans_nir.dat")
  )
  names(files_list) <- c("lengths", "reflect_par", "reflect_nir",
                         "trans_par", "trans_nir")
  file.create(files_list)

  write_dat <- function(value, file) {
    cat(value, file = file, sep = " ", append = TRUE)
    cat("\n", file = file, append = TRUE)
  }

  par.nir.lengths <- c(length(par.wl), length(nir.wl))
  write_dat(par.nir.lengths, files_list["lengths"])

  if (is_spectra(spectra_list[[1]])) {
    wl <- wavelengths(spectra_list[[1]])
  } else {
    wl <- spectra_list[[1]][, "wl"]
  }
  par.ind <- which(wl %in% par.wl)
  nir.ind <- which(wl %in% nir.wl)

  # Set up soil and wood reflectance files
  if (!is.na(soil_reflect_path)) {
    file.copy(
      soil_reflect_path,
      file.path(output.path, "soil_reflect_par.dat"),
      overwrite = TRUE
    )
  }

  if (!is.na(wood_reflect_path)) {
    file.copy(
      wood_reflect_path,
      file.path(output.path, "wood_reflect_par.dat"),
      overwrite = TRUE
    )
  }

  # Multi-PFT settings
  if (length(spectra_list) != length(trait.values)) {
    stop("Spectral data and trait.values do not have same length. ",
         "Spectral data length: ", length(spectra_list),
         "trait.values length: ", length(trait.values))
  }
  pft_names <- names(trait.values)
  if (any(!names(spectra_list) %in% pft_names)) {
    stop("Spectral data and trait.values do not have same PFT names. ",
         "Spectral data names: ", names(spectra_list),
         "trait.values names: ", pft_names)
  }
  data(pftmapping, package = "PEcAn.ED2")
  npft <- length(pft_names)
  write_dat(npft, files_list["lengths"])
  pft_numbers <- numeric(npft)
  for (i in seq_len(npft)) {
    pft_numbers[i] <- pftmapping[pftmapping$PEcAn == pft_names[i], "ED"]
    write_dat(pft_numbers[i], files_list["lengths"])
    spectra_list_pft <- spectra_list[[pft_names[i]]]
    rind <- which(colnames(spectra_list_pft) %in% c("R", "reflectance"))
    tind <- which(colnames(spectra_list_pft) %in% c("T", "transmittance"))
    write_dat(spectra_list_pft[par.ind, rind], files_list["reflect_par"])
    write_dat(spectra_list_pft[nir.ind, rind], files_list["reflect_nir"])
    write_dat(spectra_list_pft[par.ind, tind], files_list["trans_par"])
    write_dat(spectra_list_pft[nir.ind, tind], files_list["trans_nir"])
  }

  oldwd <- getwd()
  setwd(output.path)
  on.exit(setwd(oldwd), add = TRUE)

  if (!is.null(img_path)) {
    ex <- PEcAn.ED2::run_ed_singularity(
      img_path, ed2in_path,
      app = "EDR",
      singularity_args = singularity_args,
      stderr = stderr,
      stdout = stderr,
      ...
    )
  } else {
    ex <- system2(
      edr_exe_path,
      args = c("-f", ed2in_path),
      stderr = stderr,
      stdout = stderr,
      ...
    )
  }

  # Call EDR -- NOTE that this requires that the ED2IN
  if (stderr && any(grepl("fatal error|fortran runtime error", ex, ignore.case = TRUE))) {
    if (verbose_error) {
      print(ex)
    }
    stop("Error executing EDR.")
  }

  # Analyze output
  albedo <- get.EDR.output(output.path)
  # Optionally, clean up all generated files
  if(clean) {
    delete.files <- file.remove(files_list)
    # NOTE that currently, not all files are deleted (e.g. history file, copied ED2IN)
    if (!delete.files) {
      warning("Error in deleting files.")
    }
  }
  return(albedo)
} # EDR

#' Read EDR output
#'
#' @param path Path to directory containing `albedo_par/nir.dat` files
#' @export
get.EDR.output <- function(path = getwd()) {
  alb.par <- as.matrix(read.table(file.path(path, "albedo_par.dat")))[1, ]
  alb.nir <- as.matrix(read.table(file.path(path, "albedo_nir.dat")))[1, ]
  albedo <- c(alb.par, alb.nir)
  return(albedo)
}

#' Preprocess history file for EDR
#'
#' Locate history file based on path and prefix, copy to specified output
#' directory, and rename to correct time.
#' @param history.path Path to directory containing history file.
#' @param history.prefix String describing the history file prefix in
#' `history.path`. Default = 'history'
#' @param datetime POSIX date and time for run
#' @export
EDR.preprocess.history <- function(history.path, output.path, datetime, history.prefix = "history") {
  # Check inputs
  stopifnot(
    is.character(history.path),
    is.character(history.prefix)
  )

  # Extract date and time
  day          <- strftime(datetime, "%d", tz = "UTC")
  month        <- strftime(datetime, "%m", tz = "UTC")
  year         <- strftime(datetime, "%Y", tz = "UTC")
  time.history <- strftime(datetime, "%H%M%S", tz = "UTC")
  # Locate history file
  history.search <- sprintf("%1$s-S-%2$s-%3$s-%4$s",
                            history.prefix, year, month, day)
  history.name <- list.files(history.path, history.search)
  history.full.path <- file.path(history.path, history.name)
  if (length(history.name) > 1) {
    stop("Multiple history files matched")
  }
  if (length(history.name) == 0) {
    stop("No history files found")
  }
  # Copy and rename history file
  history.new.name <- gsub('([[:digit:]]{6})', time.history, history.name)
  history.new.path <- file.path(output.path, history.new.name)
  history.copy <- file.copy(history.full.path, history.new.path, overwrite = FALSE)
  if (!history.copy) {
    warning("Could not copy history with overwrite=FALSE. Attempting with overwrite=TRUE")
    history.copy <- file.copy(history.full.path, history.new.path, overwrite = TRUE)
    if (!history.copy) {
      stop("Unable to copy history file, even with overwrite=TRUE. Check permissions on both input and output directories.")
    }
  }
  history.full.prefix <- file.path(output.path, history.prefix)
  return(history.full.prefix)
} # EDR.preprocess.history