#' Add soil layers (by interpolation) if not already present and store in
#' soil data input file
#'
#' @section Notes: Splitting existing soil layers into two new layers will
#' "exhaust" the values of fields \var{\dQuote{EvapCoeff}},
#' \var{\dQuote{TranspCoeff}}, and \var{\dQuote{Imperm}}, i.e., sum remains
#' constant; and will "interpolate" the values of the fields
#' \var{\dQuote{Matricd}}, \var{\dQuote{GravelContent}}, \var{\dQuote{Sand}},
#' \var{\dQuote{Clay}}, \var{\dQuote{SoilTemp}}.
#'
#' @export
calc_RequestedSoilLayers <- function(SFSW2_prj_meta,
SFSW2_prj_inputs, runIDs_adjust, keep_old_depth = TRUE, verbose = FALSE) {
requested_soil_layers <-
SFSW2_prj_meta[["opt_input"]][["requested_soil_layers"]]
if (verbose) {
t1 <- Sys.time()
temp_call <- shQuote(match.call()[1])
print(paste0("rSFSW2's ", temp_call, ": started at ", t1))
on.exit({
print(paste0("rSFSW2's ", temp_call, ": ended after ",
round(difftime(Sys.time(), t1, units = "secs"), 2), " s"))
cat("\n")}, add = TRUE)
}
# Column name pattern of soil layers, e.g., `depth_L1`
cn_depth <- "depth_L"
# How to add different soil variables
# values will be exhausted:
sl_vars_sub <- c("EvapCoeff", "TranspCoeff", "Imperm")
# Requested layers
requested_soil_layers <- as.integer(round(requested_soil_layers))
stopifnot(requested_soil_layers > 0, diff(requested_soil_layers) > 0)
# Available layers
ids_depth <- strsplit(names(SFSW2_prj_inputs[["sw_input_soils_use"]])
[SFSW2_prj_inputs[["sw_input_soils_use"]]], "_", fixed = TRUE)
stopifnot(length(ids_depth) > 0)
var_layers <- unique(sapply(ids_depth, function(x)
paste0(x[-length(x)], collapse = "_")))
ids_depth2 <- unique(sapply(ids_depth, function(x) x[length(x)]))
use_layers <- paste0("depth_", ids_depth2)
layers_depth <- round(as.matrix(SFSW2_prj_inputs[["sw_input_soillayers"]]
[runIDs_adjust, use_layers, drop = FALSE]))
i_nodata <- apply(is.na(layers_depth), 1, all)
if (any(i_nodata)) {
layers_depth <- layers_depth[!i_nodata, ]
runIDs_adjust_ws <- runIDs_adjust[!i_nodata]
} else {
runIDs_adjust_ws <- runIDs_adjust
}
i_nodata <- apply(is.na(layers_depth), 2, all)
if (any(i_nodata))
layers_depth <- layers_depth[, !i_nodata]
ids_layers <- seq_len(dim(layers_depth)[2])
avail_sl_ids <- apply(layers_depth, 1, paste0, collapse = "x")
# Loop through runs with same layer profile and adjust
layer_sets <- unique(avail_sl_ids)
if (length(layer_sets) > 0) {
has_changed <- FALSE
sw_input_soils_data <- lapply(var_layers, function(x)
as.matrix(SFSW2_prj_inputs[["sw_input_soils"]][
runIDs_adjust_ws,
grep(x, names(SFSW2_prj_inputs[["sw_input_soils"]]))[ids_layers],
drop = FALSE]))
sw_input_soils_data2 <- NULL
for (ils in seq_along(layer_sets)) {
il_set <- avail_sl_ids == layer_sets[ils]
if (sum(il_set, na.rm = TRUE) == 0) next
# Identify which requested layers to add
ldset <- stats::na.exclude(layers_depth[which(il_set)[1], ])
req_sd_toadd <- setdiff(requested_soil_layers, ldset)
if (isTRUE(keep_old_depth)) {
req_sd_toadd <- req_sd_toadd[req_sd_toadd < max(ldset)]
}
if (length(req_sd_toadd) == 0) next
# Add identified layers
sw_input_soils_data2 <- lapply(seq_along(var_layers),
function(iv) sw_input_soils_data[[iv]][il_set, , drop = FALSE])
for (lnew in req_sd_toadd) {
ilnew <- findInterval(lnew, ldset)
il_weight <- calc_weights_from_depths(ilnew, lnew, ldset)
sw_input_soils_data2 <- lapply(seq_along(var_layers), function(iv)
add_layer_to_soil(sw_input_soils_data2[[iv]], il = ilnew,
w = il_weight, method = if (var_layers[iv] %in% sl_vars_sub)
"exhaust" else "interpolate"))
ldset <- sort(c(ldset, lnew))
}
# Update soil datafiles
lyrs <- seq_along(ldset)
irows <- runIDs_adjust_ws[il_set]
for (iv in seq_along(var_layers)) {
icol <- grep(var_layers[iv],
names(SFSW2_prj_inputs[["sw_input_soils_use"]]))[lyrs]
SFSW2_prj_inputs[["sw_input_soils"]][irows, icol] <-
round(sw_input_soils_data2[[iv]][, lyrs],
if (var_layers[iv] %in% sl_vars_sub) 4L else 2L)
SFSW2_prj_inputs[["sw_input_soils_use"]][icol] <- TRUE
}
icol <- paste0(cn_depth, lyrs)
SFSW2_prj_inputs[["sw_input_soillayers"]][irows, icol] <-
matrix(ldset, nrow = sum(il_set), ncol = length(ldset), byrow = TRUE)
has_changed <- TRUE
}
if (has_changed) {
#write data to disk
utils::write.csv(SFSW2_prj_inputs[["sw_input_soillayers"]],
file = SFSW2_prj_meta[["fnames_in"]][["fslayers"]], row.names = FALSE)
utils::write.csv(reconstitute_inputfile(
SFSW2_prj_inputs[["sw_input_soils_use"]],
SFSW2_prj_inputs[["sw_input_soils"]]),
file = SFSW2_prj_meta[["fnames_in"]][["fsoils"]], row.names = FALSE)
unlink(SFSW2_prj_meta[["fnames_in"]][["fpreprocin"]])
print(paste("'InterpolateSoilDatafileToRequestedSoilLayers':",
"don't forget to adjust lookup tables with per-layer values if",
"applicable for this project"))
}
SFSW2_prj_meta[["opt_input"]][["requested_soil_layers"]] <-
requested_soil_layers
}
list(SFSW2_prj_meta = SFSW2_prj_meta, SFSW2_prj_inputs = SFSW2_prj_inputs)
}
#' Calculate potential bare-soil evaporation coefficients
#'
#' Soil texture influence based on re-analysis of data from Wythers et al. 1999.
#' Default of \code{depth_max_bs_evap} = 15 cm from Torres et al. 2010.
#'
#' @references Torres EA, Calera A (2010) Bare soil evaporation under high
#' evaporation demand: a proposed modification to the FAO-56 model.
#' Hydrological Sciences Journal- Journal des Sciences Hydrologiques, 55,
#' 303-315.
#'
#' @references Wythers K.R., Lauenroth W.K., Paruelo J.M. (1999) Bare-Soil
#' Evaporation Under Semiarid Field Conditions. Soil Science Society of
#' America Journal, 63, 1341-1349.
#'
#' @param layers_depth A numeric vector, matrix, or data.frame. Values describe
#' the lower soil layer depths in units of centimeters.
#' @param sand A numeric vector, matrix, or data.frame. Values are sand contents
#' in units of mass-percentage / 100.
#' @param clay A numeric vector, matrix, or data.frame. Values are clay contents
#' in units of mass-percentage / 100.
#' @param depth_max_bs_evap_cm A numeric value. The maximal soil depth in
#' centimeters from which bare-soil evaporation is potentially drawing
#' moisture.
#'
#' @section Notes: Rows of soil input arguments \code{layers_depth},
#' \code{sand}, and \code{clay} correspond to sites and columns to soil
#' layers. If \code{sand} and/or \code{clay} are vectors, then they are
#' converted to 1-row matrices. If \code{layers_depth} is a vector, then it is
#' converted to a matrix with as many sites/rows as \code{sand} and
#' \code{clay} have. That is the code assumes identical soil layer depths for
#' each site. All soil input arguments must have a the same number of sites
#' and of soil layers, i.e., identical matrix dimensions.
#' @section Warning: Influence of gravel is not accounted for.
#'
#' @return A numeric matrix with potential bare-soil evaporation coefficients
#' where rows correspond to sites and columns to soil layers.
#'
#' @export
calc_BareSoilEvapCoefs <- function(layers_depth, sand, clay,
depth_max_bs_evap_cm = 15) {
#--- If inputs are not site x layers, then convert them into 1 site x layers
if (is.null(dim(sand))) {
sand <- matrix(sand, nrow = 1, ncol = length(sand))
}
if (is.null(dim(clay))) {
clay <- matrix(clay, nrow = 1, ncol = length(clay))
}
if (is.null(dim(layers_depth))) {
layers_depth <- matrix(layers_depth, nrow = dim(sand)[1],
ncol = length(layers_depth), byrow = TRUE)
}
#--- Test inputs
# - sand and clay have identical number of sites and layers
# - all soil inputs have identical number of sites and at least as many
# layers as depths
# - soil layer depths are numeric and positive -- or NA, if all deeper
# layers are NA
# - sand and clay are numeric and values between 0 and 1 -- or NA, if all
# deeper layers are NA as well
# - the sum of sand and clay is less or equal to 1
sand_and_clay <- sand + clay
stopifnot(
identical(dim(sand), dim(clay)),
identical(dim(sand)[1], dim(layers_depth)[1]),
dim(sand)[2] >= dim(layers_depth)[2],
is.numeric(layers_depth),
layers_depth > 0 | has_NAs_pooled_at_depth(layers_depth),
is.numeric(unlist(sand)),
sand >= 0 & sand <= 1 | has_NAs_pooled_at_depth(sand),
is.numeric(unlist(clay)),
clay >= 0 & clay <= 1 | has_NAs_pooled_at_depth(clay),
sand_and_clay <= 1 | has_NAs_pooled_at_depth(sand_and_clay),
is.finite(depth_max_bs_evap_cm) & depth_max_bs_evap_cm >= 0)
#--- Calculate
depth_min_bs_evap <- min(layers_depth[, 1], na.rm = TRUE)
if (depth_min_bs_evap > depth_max_bs_evap_cm) {
# all sites have first layer with coeff = 1
res <- array(1, dim = dim(sand))
res[, -1] <- 0
return(res)
}
lyrs_max_bs_evap <- t(apply(layers_depth, 1, function(x) {
xdm <- depth_max_bs_evap_cm - x
i0 <- abs(xdm) < SFSW2_glovars[["tol"]]
ld <- if (any(i0, na.rm = TRUE)) {
which(i0)
} else {
temp <- which(xdm < 0)
if (length(temp) > 0) temp[1] else length(x)
}
c(diff(c(0, x))[seq_len(ld)], rep(0L, length(x) - ld))
}))
ldepth_max_bs_evap <- rowSums(lyrs_max_bs_evap)
sand_mean <- rowSums(lyrs_max_bs_evap * sand, na.rm = TRUE) /
ldepth_max_bs_evap
clay_mean <- rowSums(lyrs_max_bs_evap * clay, na.rm = TRUE) /
ldepth_max_bs_evap
# equation from re-analysis
temp_depth <- 4.1984 + 0.6695 * sand_mean ^ 2 + 168.7603 * clay_mean ^ 2
depth_bs_evap <- pmin(pmax(temp_depth, depth_min_bs_evap, na.rm = TRUE),
depth_max_bs_evap_cm, na.rm = TRUE)
lyrs_bs_evap0 <- t(apply(depth_bs_evap - layers_depth, 1, function(x) {
i0 <- abs(x) < SFSW2_glovars[["tol"]]
ld <- if (any(i0, na.rm = TRUE)) {
which(i0)
} else {
temp <- which(x < 0)
if (length(temp) > 0) temp[1] else sum(!is.na(x))
}
ld0 <- max(0, ld - 1)
c(rep(TRUE, ld0), rep(FALSE, length(x) - ld0))
}))
# function made up to match previous cummulative distributions
temp_coeff <- 1 - exp(- 5 * layers_depth / depth_bs_evap)
temp_coeff[!lyrs_bs_evap0 | is.na(temp_coeff)] <- 1
coeff_bs_evap <- round(t(apply(cbind(0, temp_coeff), 1, diff)), 4)
coeff_bs_evap / rowSums(coeff_bs_evap, na.rm = TRUE)
}
#' Calculate bare-soil evaporation coefficients based on soil texture and store
#' in soil input file
get_BareSoilEvapCoefs <- function(SFSW2_prj_meta, SFSW2_prj_inputs,
runIDs_adjust, resume = TRUE, verbose = FALSE) {
if (verbose) {
t1 <- Sys.time()
temp_call <- shQuote(match.call()[1])
print(paste0("rSFSW2's ", temp_call, ": started at ", t1))
on.exit({
print(paste0("rSFSW2's ", temp_call, ": ended after ",
round(difftime(Sys.time(), t1, units = "secs"), 2), " s"))
cat("\n")}, add = TRUE)
}
icol_bsE <- grep("EvapCoeff", names(SFSW2_prj_inputs[["sw_input_soils_use"]]))
icol_sand <- grep("Sand_L", names(SFSW2_prj_inputs[["sw_input_soils_use"]]))
icol_clay <- grep("Clay_L", names(SFSW2_prj_inputs[["sw_input_soils_use"]]))
use_layers <- which(SFSW2_prj_inputs[["sw_input_soils_use"]][icol_sand] &
SFSW2_prj_inputs[["sw_input_soils_use"]][icol_clay])
stopifnot(length(use_layers) > 0)
do_calc <- TRUE
if (resume) {
temp <- icol_bsE[use_layers]
icols <- temp[SFSW2_prj_inputs[["sw_input_soils_use"]][temp]]
if (length(icols) > 0L) {
x <- SFSW2_prj_inputs[["sw_input_soils"]][runIDs_adjust, icols,
drop = FALSE]
do_calc <- anyNA(x) || !all(rowSums(x, na.rm = TRUE) > 0)
rm(x)
}
}
if (do_calc) {
icols <- grep("depth_L",
names(SFSW2_prj_inputs[["sw_input_soillayers"]]))[use_layers]
temp <- SFSW2_prj_inputs[["sw_input_soillayers"]][runIDs_adjust, icols,
drop = FALSE]
layers_depth <- as.matrix(temp)
sand <- SFSW2_prj_inputs[["sw_input_soils"]][runIDs_adjust, icol_sand,
drop = FALSE]
clay <- SFSW2_prj_inputs[["sw_input_soils"]][runIDs_adjust, icol_clay,
drop = FALSE]
coeff_bs_evap <- calc_BareSoilEvapCoefs(
layers_depth, sand, clay, depth_max_bs_evap_cm =
SFSW2_prj_meta[["opt_sim"]][["depth_max_bs_evap_cm"]])
#add data to sw_input_soils and set the use flags
icol <- seq_len(sum(apply(coeff_bs_evap, 2, function(x)
any(x > SFSW2_glovars[["tol"]]))))
icols_bsE_used <- icol_bsE[icol]
icols_bse_notused <- icol_bsE[-icol]
SFSW2_prj_inputs[["sw_input_soils_use"]][icols_bsE_used] <- TRUE
SFSW2_prj_inputs[["sw_input_soils"]][runIDs_adjust, icols_bsE_used] <-
round(coeff_bs_evap[, icol], 4)
SFSW2_prj_inputs[["sw_input_soils_use"]][icols_bse_notused] <- FALSE
SFSW2_prj_inputs[["sw_input_soils"]][runIDs_adjust, icols_bse_notused] <- 0
stopifnot(!is.na(SFSW2_prj_inputs[["sw_input_soils"]][runIDs_adjust,
icols_bsE_used]))
#write data to disk
utils::write.csv(reconstitute_inputfile(
SFSW2_prj_inputs[["sw_input_soils_use"]],
SFSW2_prj_inputs[["sw_input_soils"]]),
file = SFSW2_prj_meta[["fnames_in"]][["fsoils"]], row.names = FALSE)
unlink(SFSW2_prj_meta[["fnames_in"]][["fpreprocin"]])
}
SFSW2_prj_inputs
}
#' Look-up input values from spreadsheet tables
#' @export
do_prior_TableLookups <- function(SFSW2_prj_meta, SFSW2_prj_inputs,
resume = TRUE, verbose = FALSE) {
if (verbose) {
t1 <- Sys.time()
temp_call <- shQuote(match.call()[1])
print(paste0("rSFSW2's ", temp_call, ": started at ", t1))
on.exit({
print(paste0("rSFSW2's ", temp_call, ": ended after ",
round(difftime(Sys.time(), t1, units = "secs"), 2), " s"))
cat("\n")}, add = TRUE)
}
do_prior_lookup <- list(
LookupEvapCoefs = list(),
LookupTranspRegions = list(),
LookupSnowDensity = list()
)
done_prior <- rep(FALSE, length(do_prior_lookup))
names(done_prior) <- names(do_prior_lookup)
temp <- SFSW2_prj_inputs[["create_treatments"]] %in% names(do_prior_lookup)
if (any(temp)) {
do_prior_lookup[["LookupEvapCoefs"]] <- list(
flag = "LookupEvapCoefs",
pattern = "EvapCoeff",
tr_input = SFSW2_prj_inputs[["tr_input_EvapCoeff"]],
sw_input_use = "sw_input_soils_use",
sw_input = "sw_input_soils",
nvars = SFSW2_glovars[["slyrs_maxN"]],
do_fill = FALSE,
datafile = SFSW2_prj_meta[["fnames_in"]][["fsoils"]])
do_prior_lookup[["LookupTranspRegions"]] <- list(
flag = "LookupTranspRegions",
pattern = "TranspRegion",
tr_input = SFSW2_prj_inputs[["tr_input_TranspRegions"]],
sw_input_use = "sw_input_soils_use",
sw_input = "sw_input_soils",
nvars = SFSW2_glovars[["slyrs_maxN"]],
do_fill = FALSE,
datafile = SFSW2_prj_meta[["fnames_in"]][["fsoils"]])
do_prior_lookup[["LookupSnowDensity"]] <- list(
flag = "LookupSnowDensity",
pattern = "(snowd)|(SnowD_Hemisphere)",
tr_input = SFSW2_prj_inputs[["tr_input_SnowD"]],
sw_input_use = "sw_input_cloud_use",
sw_input = "sw_input_cloud",
nvars = 12 + 1,
do_fill = TRUE,
fill_pattern = "snowd",
fill_value = 76, # 76 kg/m3 = median of medians over 6 sites in
# Colorado and Wyoming: Judson, A. & Doesken, N. (2000) Density of
# Freshly Fallen Snow in the Central Rocky Mountains. Bulletin of the
# American Meteorological Society, 81, 1577-1587.
datafile = SFSW2_prj_meta[["fnames_in"]][["fclimnorm"]])
for (pc in do_prior_lookup) {
if (any(SFSW2_prj_inputs[["create_treatments"]] == pc$flag)) {
# lookup values per category for each simulation run and copy values
# to datafile
temp1 <- SFSW2_prj_inputs[["sw_input_experimentals_use"]][pc$flag]
temp <- unique(SFSW2_prj_inputs[["sw_input_experimentals"]][, pc$flag])
temp2 <- length(temp) == 1L
# Lookup prior to do_OneSite() only if option is off in
# sw_input_experimentals or constant
if (!temp1 || (temp1 && temp2)) {
if (resume) {
# Determine whether lookup already carried out and stored to file
sw_input_use <- SFSW2_prj_inputs[[pc$sw_input_use]]
icols <- grep(pc$pattern, names(sw_input_use))
icols <- icols[sw_input_use[icols]]
temp <- SFSW2_prj_inputs[[pc$sw_input]][, icols, drop = FALSE]
if (all(!apply(is.na(temp), 2, all))) {
# if no layer has only NAs for which the _use flag is on, then
# consider as completed
done_prior[pc$flag] <- TRUE
next
}
}
if (verbose)
print(paste(Sys.time(), ": performing", shQuote(pc$flag)))
temp <- SFSW2_prj_inputs[["sw_input_experimentals_use"]][pc$flag]
trtype <- if (temp) {
unique(SFSW2_prj_inputs[["sw_input_experimentals"]][, pc$flag])
} else {
SFSW2_prj_inputs[["sw_input_treatments"]][, pc$flag]
}
if (any(is.na(trtype)))
stop("ERROR: ", pc$flag, " column cannot have any NAs.")
if (!all(unique(trtype) %in% rownames(pc$tr_input)))
stop("ERROR: ", pc$flag, " column values do not match up with ",
"trfile. ", pc$flag, " row names.")
tempdat <- try(get.LookupFromTable(
pattern = pc$pattern,
trtype = trtype,
tr_input = pc$tr_input,
sw_input_use = SFSW2_prj_inputs[[pc$sw_input_use]],
sw_input = SFSW2_prj_inputs[[pc$sw_input]],
nvars = pc$nvars))
done_prior[pc$flag] <- !inherits(tempdat, "try-error")
if (done_prior[pc$flag]) {
if (!is.null(pc$do_fill) && pc$do_fill) {
tempdat <- fill_empty(tempdat, pattern = pc$fill_pattern,
fill = pc$fill_value)
}
SFSW2_prj_inputs[[pc$sw_input_use]] <- tempdat$sw_input_use
SFSW2_prj_inputs[[pc$sw_input]] <- tempdat$sw_input
#write data to datafile
utils::write.csv(reconstitute_inputfile(
SFSW2_prj_inputs[[pc$sw_input_use]],
SFSW2_prj_inputs[[pc$sw_input]]), file = pc$datafile,
row.names = FALSE)
unlink(SFSW2_prj_meta[["fnames_in"]][["fpreprocin"]])
}
} else {
done_prior[pc$flag] <- FALSE
}
}
}
}
SFSW2_prj_inputs[["done_prior"]] <- done_prior
SFSW2_prj_inputs
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.