#' 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
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.