#------ Deprecated functions ------
#' Deprecated pedotransfer functions to convert between soil moisture
#' (volumetric water content, \var{VWC}) and soil water potential (\var{SWP})
#'
#' @param sand A numeric value or vector. Sand content of the soil layer(s) as
#' fractional value in \code{[0,1]}.
#' @param clay A numeric value or vector. Clay content of the soil layer(s) as
#' fractional value in \code{[0,1]}.
#' @param thetas Soon obsolete ... (see `feature_swrc`)
#' @param psis Soon obsolete ... (see `feature_swrc`)
#' @param b Soon obsolete ... (see `feature_swrc`)
#' @param MPa_toBar Soon obsolete ... (see `feature_swrc`)
#' @param bar_conversion Soon obsolete ... (see `feature_swrc`)
#' @param bar_toMPa Soon obsolete ... (see `feature_swrc`)
#' @param ... Additional arguments.
#'
#' @references Cosby, B. J., G. M. Hornberger, R. B. Clapp, and T. R. Ginn.
#' 1984. A statistical exploration of the relationships of soil moisture
#' characteristics to the physical properties of soils. Water Resources Research
#' 20:682-690.
#'
#' @seealso The use of these functions is deprecated;
#' please use `ptf_estimate()` and `swrc_conversion()` instead.
#'
#' @name pedotransfer
#' @md
NULL
#' @rdname pedotransfer
#' @section Note: either \code{swp} or \code{sand}/\code{clay} needs be a
#' single value
pdf_to_vwc <- function(swp, sand, clay, thetas, psis, b, MPa_toBar = -10,
bar_conversion = 1024) {
thetas * (psis / (swp * MPa_toBar * bar_conversion)) ^ (1 / b) / 100
}
#' @rdname pedotransfer
#' @section Note: either \code{vwc} or \code{sand}/\code{clay} needs be a
#' single value
pdf_to_swp <- function(vwc, sand, clay, thetas, psis, b, bar_toMPa = -0.1,
bar_conversion = 1024) {
psis / ((vwc * 100 / thetas) ^ b * bar_conversion) * bar_toMPa
}
pedotransfer <- function(x, sand, clay, pdf) {
stopifnot(
length(sand) > 0,
length(sand) == length(clay)
)
sand <- rSW2utils::finite01(sand, NA, NA)
clay <- rSW2utils::finite01(clay, NA, NA)
if (any(complete.cases(sand, clay))) {
thetas <- -14.2 * sand - 3.7 * clay + 50.5
psis <- 10 ^ (-1.58 * sand - 0.63 * clay + 2.17)
b <- -0.3 * sand + 15.7 * clay + 3.10
if (any(b <= 0, na.rm = TRUE))
stop("Pedotransfer for soil texture with b <= 0 is not possible.")
np_x <- NROW(x) * NCOL(x)
if (NROW(x) == 1 || NCOL(x) == 1) {
# cases 1-4
if (np_x == 1 || length(sand) == 1) {
# cases 1-3
res <- pdf(x, sand, clay, thetas, psis, b)
} else {
# case 4; Note: case 3 could also be calculated with the code for
# case 4, but is much slower, unless x is a data.frame with one column
temp <- lapply(
x,
pdf,
sand = sand,
clay = clay,
thetas = thetas,
psis = psis,
b = b
)
res <- matrix(unlist(temp), nrow = np_x, byrow = TRUE)
}
} else {
# cases 5-6
dx <- dim(x)
if (length(sand) == 1) {
# case 5
res <- vapply(seq_len(dx[2]), function(d) {
pdf(x[, d], sand, clay, thetas, psis, b)
}, rep(1, dx[1]), USE.NAMES = FALSE)
} else {
# case 6
stopifnot(dx[2] == length(sand))
res <- vapply(seq_len(dx[2]), function(d) {
pdf(x[, d], sand[d], clay[d], thetas[d], psis[d], b[d])
}, rep(1, dx[1]), USE.NAMES = FALSE)
}
}
} else {
res <- x
res[] <- NA
}
# if SWP then in units of MPa [-Inf, 0]; if VWC then in units of m3/m3 [0, 1]
res
}
#' Calculate volumetric water content from soil water potential and soil texture
#' @rdname pedotransfer
#' @param swp A numeric value, vector, or 2-dimensional object
#' (matrix or data.frame). The soil water potential (of the soil matrix) in
#' units of \var{MPa}, i.e., the soil without the volume of rock and gravel.
#'
#' @return Volumetric water content in units of m^3 (of water) / m^3 (of soil)
#' \code{[0, 1]}. There are six use cases:\enumerate{
#' \item 1) \itemize{
#' \item Input: \code{SWP} [single value]; \code{sand} and \code{clay}
#' [single values]
#' \item Output: \code{VWC} [single value]}
#' \item 2) \itemize{
#' \item Input: \code{SWP} [single value]; \code{sand} and \code{clay}
#' [vectors of length d]
#' \item Output: \code{VWC} [vector of length d]}
#' \item 3) \itemize{
#' \item Input: \code{SWP} [vector of length l]; \code{sand} and
#' \code{clay} infraction [single values]
#' \item Output: \code{VWC} [vector of length l]}
#' \item 4) \itemize{
#' \item Input: \code{SWP} [vector of length l]; \code{sand} and
#' \code{clay} [vectors of length d]
#' \item Output: \code{VWC} [l x d matrix] where \code{SWP} is
#' repeated for each column}
#' \item 5) \itemize{
#' \item Input: \code{SWP} [l x d matrix]; \code{sand} and \code{clay}
#' [single values]
#' \item Output: \code{VWC} [l x d matrix]}
#' \item 6) \itemize{
#' \item Input: \code{SWP} [l x d matrix]; \code{sand} and \code{clay}
#' [vectors of length d]
#' \item Output: \code{VWC} [l x d matrix], \code{sand} and \code{clay}
#' vectors are repeated for each row}
#' }
#' @export
SWPtoVWC <- function(swp, sand, clay, ...) {
.Deprecated("swrc_swp_to_vwc")
swrc_swp_to_vwc(
swp_MPa = swp,
sand = sand,
clay = clay,
...
)
}
SWPtoVWC_old <- function(swp, sand, clay) {
pedotransfer(swp, sand, clay, pdf = pdf_to_vwc)
}
#' Calculate soil water potential from volumetric water content and soil texture
#' @rdname pedotransfer
#' @param vwc A numeric value, vector, or 2-dimensional object
#' (matrix or data.frame). The matric soil moisture, i.e., reduced by the
#' volume of rock and gravel.
#'
#' @return Soil water potential in units of \var{MPa} \code{[-Inf, 0]}.
#' There are six use cases: \enumerate{
#' \item 1) \itemize{
#' \item Input: \code{VWC} [single value]; \code{sand} and \code{clay}
#' [single values]
#' \item Output: \code{SWP} [single value]}
#' \item 2) \itemize{
#' \item Input: \code{VWC} [single value]; \code{sand} and \code{clay}
#' [vectors of length d]
#' \item Output: \code{SWP} [vector of length d]}
#' \item 3) \itemize{
#' \item Input: \code{VWC} [vector of length l]; \code{sand} and
#' \code{clay} in fraction [single values]
#' \item Output: \code{SWP} [vector of length l]}
#' \item 4) \itemize{
#' \item Input: \code{VWC} [vector of length l]; \code{sand} and
#' \code{clay} [vectors of length d]
#' \item Output: \code{SWP} [l x d matrix] where \code{VWC} is repeated for
#' each column}
#' \item 5) \itemize{
#' \item Input: \code{VWC} [l x d matrix]; \code{sand} and \code{clay}
#' [single values]
#' \item Output: \code{SWP} [l x d matrix]}
#' \item 6) \itemize{
#' \item Input: \code{VWC} [l x d matrix]; \code{sand} and \code{clay}
#' [vectors of length d]
#' \item Output: \code{SWP} [l x d matrix], \code{sand} and \code{clay}
#' vectors are repeated for each row}
#' }
#' @export
VWCtoSWP <- function(vwc, sand, clay, ...) {
.Deprecated("swrc_vwc_to_swp")
swrc_vwc_to_swp(
vwcBulk = vwc,
sand = sand,
clay = clay,
...
)
}
VWCtoSWP_old <- function(vwc, sand, clay) {
pedotransfer(vwc, sand, clay, pdf = pdf_to_swp)
}
#------ SWRC parameters & pedotransfer functions ------
# MAINTENANCE:
# Notes for implementing a new PTF "YYY" in `rSOILWAT2`
# * new `ptf_YYY_for_XXX()`
# * new `ptf_YYY_availability(verbose = interactive(), ...)`
# * update `ptfs_implemented_by_rSW2()`
# * update `rSW2_check_SWRC_vs_PTF()`
# * update `rSW2_SWRC_PTF_estimate_parameters()`
# * update `check_ptf_availability()`
# * update examples and unit tests to utilize new functions
#
# Notes for implementing a new SWRC "XXX" and/or PTF "YYY" in `SOILWAT2`
# 1) SOILWAT2: see notes in SOILWAT2/src/SW_Site.h
# 2) rSOILWAT2: everything should automatically work with new XXX/YYY
# * update examples and unit tests to utilize new functions
#' Functionality for Soil Water Retention Curves (`SWRC`)
#'
#' @description
#' `SWRCs` convert between soil water content and soil water potential
#' using a set of parameters, see [swrc_swp_to_vwc()] and [swrc_vwc_to_swp()].
#'
#' The `SWRC` parameters may be estimated from soil properties with suitable
#' pedotransfer functions `PTFs`, see [ptf_estimate()].
#'
#' The `SWRC` parameters can be checked for consistency with [check_swrcp()].
#'
#'
#' @param sand A numeric value or vector.
#' Sand content of the matric soil component
#' (< 2 mm fraction; units of `[g/g]`) of each soil layer.
#' @param clay A numeric value or vector.
#' Clay content of the matric soil component
#' (< 2 mm fraction; units of `[g/g]`) of each soil layer.
#' @param fcoarse A numeric value or vector.
#' Coarse fragments, e.g., gravel, (> 2 mm; units of `[m3/m3]`)
#' relative to the whole soil of each soil layer.
#' `fcoarse` is required, for instance, to translate between
#' values relative to the matric soil component (< 2 mm fraction) and
#' relative to the whole soil (matric soil plus coarse fragments).
#' @param bdensity A numeric value or vector.
#' Density of the whole soil
#' (matric soil plus coarse fragments; units `[g/cm3]`).
#' @param layer_width A numeric value or vector.
#' Depth interval, width, of each soil layer (units of `cm`).
#' `layer_width` is required to translate between
#' soil water content of a soil layer and volumetric water content.
#' @param swrc_name An character string or vector.
#' The selected `SWRC` name
#' (one of [swrc_names()], with default `"Campbell1974"`).
#' @param ptf_name An character string or vector.
#' The selected `PTF` name
#' (one of [ptf_names()], with default `"Cosby1984AndOthers"`).
#' @param swrcp A numeric vector or matrix.
#' The parameters of a selected `SWRC`;
#' each row represents one `SWRC`, e.g., one per soil layer.
#' @param swrc A named list.
#' Contains all necessary elements of a `SWRC`,
#' i.e., `name` (short for `swrc_name`) and `swrcp`,
#' or all necessary elements to estimate parameters of a `SWRC` given
#' soil parameters, i.e., `swrc_name` and `ptf_name`.
#' @param fail A logical value.
#' Issue a warning (`FALSE`, default) or throw an error (`TRUE`)
#' if request fails.
#' @param verbose A logical value. If `TRUE`, then display
#' `SOILWAT2` internal warnings and other messages.
#' @param ... Additional function arguments passed on or ignored.
#'
#' @section Details:
#' [swrc_names()] lists implemented `SWRCs`;
#' [ptf_names()] lists implemented `PTFs`.
#'
#' @inherit ptf_Rosetta_for_vanGenuchten1980 references
#' @inherit ptf_neuroFX2021_for_FXW references
#' @references
#' Cosby, B. J., G. M. Hornberger, R. B. Clapp, & T. R. Ginn. 1984.
#' A statistical exploration of the relationships of soil moisture
#' characteristics to the physical properties of soils.
#' Water Resources Research, 20:682-690, \doi{10.1029/WR020i006p00682}
#'
#' @seealso
#' [swrc_names()],
#' [ptf_names()],
#' [check_ptf_availability()],
#' [ptf_estimate()],
#' [check_swrcp()],
#' [swrc_swp_to_vwc()],
#' [swrc_vwc_to_swp()]
#'
#' @name SWRCs
#' @md
NULL
#' List Soil Water Retention Curves `SWRCs`
#'
#' @return An integer vector with names of implemented `SWRCs`
#'
#' @details Notes:
#' The integer values may change with new versions of `SOILWAT2.`
#'
#' @seealso [`SWRCs`], [ptf_names()], [check_ptf_availability()]
#'
#' @md
#' @export
swrc_names <- function() {
rSW2_glovars[["kSOILWAT2"]][["SWRC_types"]]
}
#' List Pedotransfer Functions `PTFs`
#'
#' @return An named integer vector with names of implemented `PTFs`
#'
#' @details Notes:
#' The integer values may change with new versions of `SOILWAT2.`
#'
#' @seealso [`SWRCs`], [swrc_names()], [check_ptf_availability()]
#'
#' @md
#' @export
ptf_names <- function() {
tmp1 <- ptfs_implemented_by_SW2(names_only = FALSE)
tmp2 <- ptfs_implemented_by_rSW2()
c(
tmp1,
stats::setNames(max(tmp1) + seq_along(tmp2), tmp2)
)
}
#' List `PTFs` implemented by `rSOILWAT2`
#' @md
ptfs_implemented_by_rSW2 <- function() {
c(
# `Rosetta3` estimates parameters of `vanGenuchten1980` SWRC
"Rosetta3",
# `neuroFX2021` estimates parameters of `FXW` SWRC
"neuroFX2021"
)
}
#' List `PTFs` implemented by `SOILWAT2`
#'
#' @param names_only A logical value, see `return` value.
#'
#' @return An named integer vector (if not `names_only`)
#' with or a character vector (if `names_only`) names of implemented `PTFs`.
#'
#' @md
ptfs_implemented_by_SW2 <- function(names_only = FALSE) {
res <- rSW2_glovars[["kSOILWAT2"]][["PTF_types"]]
if (isTRUE(names_only)) names(res) else res
}
#' Standardize a `SWRC` name
#'
#' `"Campbell1974"` is the backwards compatible default.
#'
#' @md
#' @noRd
std_swrc <- function(swrc_name) {
if (missing(swrc_name) || is.null(swrc_name) || all(is.na(swrc_name))) {
"Campbell1974"
} else {
as.character(swrc_name)
}
}
#' Standardize a `PTF` name
#'
#' `"Cosby1984AndOthers"` is the backwards compatible default.
#'
#' @md
#' @noRd
std_ptf <- function(ptf_name) {
if (missing(ptf_name) || is.null(ptf_name) || all(is.na(ptf_name))) {
"Cosby1984AndOthers"
} else {
as.character(ptf_name)
}
}
#' Translate a `SWRC` name to its `SOILWAT2` internal integer code
#'
#' @return An integer value. `NA` if `swrc_name` is not implemented.
#'
#' @md
#' @noRd
encode_name2swrc <- function(swrc_name, fail = TRUE) {
res <- as.integer(unname(swrc_names()[std_swrc(swrc_name)]))
if (isTRUE(fail) && anyNA(res)) {
stop(
"Requested SWRC ",
shQuote(swrc_name[is.na(res)]),
" is not available."
)
}
res
}
#' Translate a `PTF` name to its `SOILWAT2` internal integer code
#'
#' @return An integer value. `NA` if `ptf_name` is not implemented.
#'
#' @md
#' @noRd
encode_name2ptf <- function(ptf_name, fail = TRUE) {
res <- as.integer(unname(ptfs_implemented_by_SW2()[std_ptf(ptf_name)]))
if (isTRUE(fail) && anyNA(res)) {
stop("Requested PTF ", shQuote(ptf_name[is.na(res)]), " is not available.")
}
res
}
#' Matching pairs of implemented `SWRCs` and `PTFs`
#'
#' @inheritParams SWRCs
#'
#' @return A `data.frame` with two columns `SWRC` and `PTF` where each
#' row contains a matching pair of `SWRC` and `PTF` that are implemented.
#'
#' @examples
#' # Data frame of SWRC-PTF combinations
#' df_swrc_ptfs <- rSOILWAT2::list_matched_swrcs_ptfs()
#'
#' # List of SWRC-PTF combinations
#' list_swrcs_ptfs <- unname(as.list(as.data.frame(t(df_swrc_ptfs))))
#'
#' # Available SWRC-PTF combinations
#' has_ptf <- check_ptf_availability(df_swrc_ptfs[, "PTF"])
#' df_swrc_ptfs[has_ptf, , drop = FALSE]
#' list_swrcs_ptfs[has_ptf]
#'
#' @md
#' @export
list_matched_swrcs_ptfs <- function(swrc_name = names(swrc_names())) {
res <- expand.grid(
SWRC = std_swrc(swrc_name),
PTF = names(ptf_names()),
stringsAsFactors = FALSE,
KEEP.OUT.ATTRS = FALSE
)
ids <- check_SWRC_vs_PTF(res[, "SWRC"], res[, "PTF"])
res <- res[ids, , drop = FALSE]
rownames(res) <- NULL
res
}
#' Check whether `PTF` and `SWRC` are compatible and implemented
#'
#' @inheritParams SWRCs
#'
#' @return A logical vector.
#'
#' @examples
#' check_SWRC_vs_PTF("Campbell1974", c("Cosby1984", "Rosetta3"))
#'
#' @md
#' @export
check_SWRC_vs_PTF <- function(swrc_name, ptf_name, fail = FALSE) {
swrc_names <- std_swrc(swrc_name)
ptf_names <- std_ptf(ptf_name)
n_swrcs <- length(swrc_names)
n_ptfs <- length(ptf_names)
if (n_swrcs == 1L && n_ptfs > 1L) {
swrc_names <- rep(swrc_names, n_ptfs)
n_swrcs <- n_ptfs
} else if (n_swrcs > 1L && n_ptfs == 1L) {
ptf_names <- rep(ptf_names, n_swrcs)
n_ptfs <- n_swrcs
}
stopifnot(n_swrcs == n_ptfs)
# Check if SWRC/PTF implemented in rSOILWAT2
has_rSW2 <- rSW2_check_SWRC_vs_PTF(swrc_names, ptf_names)
# Check if SWRC/PTF implemented in SOILWAT2
has_SW2 <- mapply(
function(s, p) .Call(C_sw_check_SWRC_vs_PTF, s, p),
swrc_names,
ptf_names
)
# SWRC/PTF implemented in either rSOILWAT2 or SOILWAT2
res <- unname(has_rSW2 | has_SW2)
if (!all(res) && isTRUE(fail)) {
ids <- which(!res)
stop(
toString(paste(swrc_names[ids], ptf_names[ids], collapse = "-")),
" are not available or incompatible."
)
}
res
}
#' Check whether `PTF` and `SWRC` are compatible and implemented in `rSOILWAT2`
#'
#' @inheritParams SWRCs
#'
#' @md
#' @noRd
rSW2_check_SWRC_vs_PTF <- function(swrc_name, ptf_name) {
swrc_name <- std_swrc(swrc_name)
ptf_name <- std_ptf(ptf_name)
ptf_name %in% ptfs_implemented_by_rSW2() & (
swrc_name == "vanGenuchten1980" & ptf_name == "Rosetta3" |
swrc_name == "FXW" & ptf_name == "neuroFX2021"
)
}
#' Estimate `SWRC` parameters from soil texture with a pedotransfer function
#'
#' @inheritParams SWRCs
#' @param ... Additional parameters passed to selected `PTF` function.
#'
#' @section Notes:
#' [swrc_names()] lists implemented `SWRCs`;
#' [ptf_names()] lists implemented `PTFs`; and
#' [check_ptf_availability()] checks availability of `PTFs`.
#'
#' @section Notes:
#' The soil parameters `sand`, `clay`, `fcoarse`, and `bdensity` must be of
#' the same length, i.e., represent one soil (length 1) or
#' multiple soil (layers) (length > 1); however, `bdensity` may be `NULL`.
#' The arguments selecting `SWRC` (`swrc_name`) and `PTF` (`ptf_name`)
#' are recycled for multiple soil layers.
#'
#' @inherit SWRCs references
#'
#' @return `swrcp`, i.e,.
#' a numeric matrix where rows represent soil (layers) and
#' columns represent a fixed number of `SWRC` parameters.
#' The interpretation is dependent on the selected `SWRC`, see
#' `SOILWAT2` input file `swrc_param.in`
# nolint start: line_length_linter.
#' (
#' `system.file("extdata", "example1", "Input", "swrc_params.in", package = "rSOILWAT2")`
#' ).
# nolint end
#'
#' @examples
#' ptf_estimate(sand = c(0.5, 0.3), clay = c(0.2, 0.1), fcoarse = c(0, 0))
#'
#' soils <- swSoils_Layers(rSOILWAT2::sw_exampleData)
#'
#' # Use PTF "Cosby1984" to estimate parameters of SWRC "Campbell1974"
#' ptf_estimate(
#' sand = soils[, "sand_frac"],
#' clay = soils[, "clay_frac"],
#' fcoarse = soils[, "gravel_content"],
#' swrc_name = "Campbell1974",
#' ptf_name = "Cosby1984"
#' )
#'
#' # Use PTF "Rosetta3" to estimate parameters of SWRC "vanGenuchten1980"
#' if (check_ptf_availability("Rosetta3")) {
#' ptf_estimate(
#' sand = soils[, "sand_frac"],
#' clay = soils[, "clay_frac"],
#' fcoarse = soils[, "gravel_content"],
#' bdensity = soils[, "bulkDensity_g/cm^3"],
#' swrc_name = "vanGenuchten1980",
#' ptf_name = "Rosetta3"
#' )
#' }
#'
#' # Use PTF "neuroFX2021" to estimate parameters of SWRC `FXW`
#' \dontrun{
#' # Set neuroFX2021 file path, see details in `ptf_neuroFX2021_for_FXW()`
#' options(RSW2_FILENEUROFX2021 = "path/to/sscbd.RData")
#' }
#'
#' if (check_ptf_availability("neuroFX2021")) {
#' ptf_estimate(
#' sand = soils[, "sand_frac"],
#' clay = soils[, "clay_frac"],
#' fcoarse = soils[, "gravel_content"],
#' bdensity = soils[, "bulkDensity_g/cm^3"],
#' swrc_name = "FXW",
#' ptf_name = "neuroFX2021"
#' )
#' }
#'
#' @md
#' @export
ptf_estimate <- function(
sand,
clay,
fcoarse,
bdensity = NULL,
swrc_name,
ptf_name,
fail = FALSE,
...
) {
#--- Check for consistency between SWRC and PTF
swrc_name <- std_swrc(swrc_name)[1]
ptf_name <- std_ptf(ptf_name)[1]
check_SWRC_vs_PTF(swrc_name, ptf_name, fail = TRUE)
#--- Determine whether we use a C- or R-implemented PTF
swrcp <- if (ptf_name %in% ptfs_implemented_by_rSW2()) {
rSW2_SWRC_PTF_estimate_parameters(
ptf_name = ptf_name,
sand = sand,
clay = clay,
fcoarse = fcoarse,
bdensity = bdensity,
fail = fail,
...
)
} else {
.Call(
C_rSW2_SWRC_PTF_estimate_parameters,
ptf_type = rep_len(encode_name2ptf(ptf_name), length(sand)),
sand = sand,
clay = clay,
fcoarse = fcoarse,
bdensity = bdensity
)
}
#--- Check validity of estimated SWRCp
if (!all(check_swrcp(swrc_name, swrcp))) {
msg <- "Some estimated parameters failed checks."
if (isTRUE(fail)) stop(msg) else warning(msg)
}
swrcp
}
#' Check availability of `PTFs`
#'
#' `PTFs` implemented in `SOILWAT2` are always available;
#' `PTFs` implemented in `rSOILWAT2` may have additional requirements, e.g.,
#' live internet connection or access to specific data files.
#'
#' @param ptfs A character vector. `PTF` names to be checked;
#' defaults to [ptf_names()].
#' @param verbose A logical value.
#'
#' @return A named logical vector with current availability of `PTFs`;
#' `PTFs` that are not implemented return `NA`.
#'
#' @examples
#' check_ptf_availability()
#' check_ptf_availability("neuroFX2021")
#' check_ptf_availability("nonexistent_PTF")
#'
#' @export
#' @md
check_ptf_availability <- function(
ptfs = names(ptf_names()),
verbose = interactive()
) {
res <- rep(NA, length(ptfs))
names(res) <- ptfs
rptfs <- ptfs_implemented_by_rSW2()
is_rptf <- ptfs %in% rptfs
# PTFs implemented in SOILWAT2 are always available
tmp <- ptfs %in% ptfs_implemented_by_SW2(names_only = TRUE) & !is_rptf
res[tmp] <- TRUE
# Check requested PTFs implemented in R
has_rptfs <- vapply(
ptfs[is_rptf],
function(ptf) {
switch(
EXPR = ptf,
Rosetta3 = ptf_Rosetta3_availability(verbose = verbose),
neuroFX2021 = ptf_neuroFX2021_availability(verbose = verbose),
NA
)
},
FUN.VALUE = NA,
USE.NAMES = TRUE
)
res[names(has_rptfs)] <- has_rptfs
res
}
#' Estimate parameters of selected soil water retention curve (`SWRC`)
#' using selected pedotransfer function (`PTF`) that are implemented in `R`
#'
#' @inheritParams ptf_estimate
#' @param fail A logical value. If requested `PTF` fails,
#' then issue a warning (`FALSE`) or throw an error (`TRUE`, default).
#'
#' @return `swrcp`, i.e,.
#' a numeric matrix where rows represent soil (layers) and
#' columns represent a fixed number of `SWRC` parameters.
#' The interpretation is dependent on the selected `SWRC`.
#' However, return value is `NULL`
#' only if `fail` is `FALSE` and requested `PTF` is not implemented in `R`.
#'
#' @inherit SWRCs references
#'
#' @section Details:
#' [ptf_estimate()] is the function that should be directly called; this here
#' is an internal helper function.
#'
#' @section Notes:
#' See `SWRC_PTF_estimate_parameters()` in `SOILWAT2` for `PTFs`
#' implemented in C.
#'
#' @md
rSW2_SWRC_PTF_estimate_parameters <- function( # nolint: object_length_linter.
ptf_name,
sand,
clay,
fcoarse,
bdensity = NULL,
fail = TRUE,
...
) {
ptf_name <- std_ptf(ptf_name)[[1L]]
has_ptf <- ptf_name %in% ptfs_implemented_by_rSW2()
list_soilargs <- list(
sand = sand,
clay = clay,
bdensity = bdensity
)
if (has_ptf && ptf_name == "Rosetta3") {
dots <- list(...)
dots[["version"]] <- if ("version" %in% names(dots)) {
as.character(dots[["version"]])
} else {
"3"
}
do.call(
ptf_Rosetta_for_vanGenuchten1980,
args = c(list_soilargs, dots)
)
} else if (has_ptf && ptf_name == "neuroFX2021") {
do.call(
ptf_neuroFX2021_for_FXW,
args = c(list_soilargs, list(...))
)[["mean"]]
} else {
msg <- paste("PTF", shQuote(ptf_name), "is not implemented in rSOILWAT2.")
if (isTRUE(fail)) stop(msg) else warning(msg)
NULL
}
}
#' Estimate van Genuchten 1980 `SWRC` parameters using `Rosetta` live `API`
#'
#' @inheritParams SWRCs
#' @param version A character string that selects a `Rosetta` version.
#'
#' @return `swrcp`, i.e,.
#' a numeric matrix where rows represent soil (layers) and
#' columns represent a fixed number of `SWRC` parameters: \itemize{
#' \item `swrcp[0]` (`theta_r`): residual volumetric water content
#' of the matric component (units of `[cm / cm]`)
#' \item `swrcp[1]` (`theta_s`): saturated volumetric water content
#' of the matric component (units of `[cm / cm]`)
#' \item `swrcp[2]` (`alpha`): related to the inverse of
#' air entry suction (units of `[cm-1]`)
#' \item `swrcp[3]` (`n`): measure of the pore-size distribution `[-]`
#' \item `swrcp[4]` (`K_sat`): saturated hydraulic conductivity `[cm / day]`
#' }
#'
#' @references
#' Mualem, Y. 1976. A new model for predicting the hydraulic conductivity of
#' unsaturated porous media.
#' Water Resources Research, 12:513-522, \doi{10.1029/WR012i003p00513}
#' @references
#' van Genuchten, M. T. 1980. A Closed-form Equation for Predicting the
#' Hydraulic Conductivity of Unsaturated Soils.
#' Soil Science Society of America Journal, 44:892-898,
#' \doi{10.2136/sssaj1980.03615995004400050002x}
#' @references
#' Zhang, Y., & Schaap, M. G. 2017. Weighted recalibration of the
#' Rosetta pedotransfer model with improved estimates of
#' hydraulic parameter distributions and summary statistics (Rosetta3).
#' Journal of Hydrology, 547:39-53, \doi{10.1016/j.jhydrol.2017.01.004}
#'
#' @section Details:
#' [ptf_estimate()] is the function that should be directly called; this here
#' is an internal helper function.
#'
#' @section Notes:
#' This function calls [soilDB::ROSETTA()] and
#' a live internet connection is required to access `Rosetta`.
#'
#' @seealso [soilDB::ROSETTA()]
#'
#' @md
ptf_Rosetta_for_vanGenuchten1980 <- function( # nolint: object_length_linter.
sand,
clay,
bdensity = NULL,
version = c("3", "1", "2"),
verbose = interactive(),
...
) {
stopifnot(ptf_Rosetta3_availability(verbose = verbose))
version <- match.arg(version)
if (verbose) {
message("Connecting live to ROSETTA API...")
}
tmp_txt <- 100 * data.frame(
sand = sand,
silt = 1 - (sand + clay),
clay = clay
)
var_txt <- c("sand", "silt", "clay")
tmp <- if (is.null(bdensity)) {
soilDB::ROSETTA(tmp_txt, vars = var_txt, v = version)
} else {
soilDB::ROSETTA(
cbind(tmp_txt, bdensity = bdensity),
vars = c(var_txt, "bdensity"),
v = version
)
}
unname(data.matrix(data.frame(
tmp[, c("theta_r", "theta_s")],
10 ^ tmp[, "alpha"],
10 ^ tmp[, "npar"],
10 ^ tmp[, "ksat"],
0
)))
}
# Checks availability of `Rosetta3` `PTF`
#
# Note: `check_ptf_availability()` requires function name
# to match pattern "ptf_XXX_availability" where XXX = name of PTF
ptf_Rosetta3_availability <- function(verbose = interactive(), ...) {
tmp <- c(
requireNamespace("soilDB"),
requireNamespace("curl") && curl::has_internet()
)
res <- all(tmp)
if (!res && verbose) {
if (!tmp[1]) {
message(
"`ptf_Rosetta3_availability()`: ",
"R package 'soilDB' is not available."
)
}
if (!tmp[2]) {
message(
"`ptf_Rosetta3_availability()`: ",
"R package 'curl' is not available or there is no live internet."
)
}
}
res
}
#' Estimate `FXW` `SWRC` parameters using `neuroFX2021`
#'
#' @inheritParams SWRCs
#' @param file_neuroFX2021 A character string that contains the file name with
#' full path of the `neuroFX2021` R object provided by Rudiyanto et al. 2021;
#' The path to the appropriate file can be set per R session
#' via option `"RSW2_FILENEUROFX2021"`, see additional details.
#'
#' @return `swrcp`, i.e,.
#' a numeric matrix where rows represent soil (layers) and
#' columns represent a fixed number of `SWRC` parameters: \itemize{
#' \item `swrcp[0]` (`theta_s`): saturated volumetric water content
#' of the matric component (units of `[cm / cm]`)
#' \item `swrcp[1]` (`alpha`): shape parameter (units of `[cm-1]`)
#' \item `swrcp[2]` (`n`): shape parameter `[-]`
#' \item `swrcp[3]` (`m`): shape parameter `[-]`
#' \item `swrcp[4]` (`K_sat`): saturated hydraulic conductivity `[cm / day]`
#' \item `swrcp[5]` (`L`): tortuosity/connectivity parameter `[-]`
#' }
#'
#' @references
#' Rudiyanto, Minasny, B., Chaney, N. W., Maggi, F., Goh Eng Giap, S.,
#' Shah, R. M., Fiantis, D., & Setiawan, B. I. 2021.
#' Pedotransfer functions for estimating soil hydraulic properties from
#' saturation to dryness.
#' Geoderma, 403:115194, \doi{10.1016/j.geoderma.2021.115194}
#' @references
#' Fredlund, D. G., & Xing, A. 1994.
#' Equations for the soil-water characteristic curve.
#' Canadian Geotechnical Journal, 31: 512–532, \doi{10.1139/t94-061}
#' @references
#' Wang, Y., Jin, M., & Deng, Z. 2018.
#' Alternative model for predicting soil hydraulic conductivity over
#' the complete moisture range.
#' Water Resources Research, 54:6860–6876, \doi{10.1029/2018WR023037}
#'
#' @section Details:
#' [ptf_estimate()] is the function that should be directly called; this here
#' is an internal helper function.
#'
#' @section Details:
#' This function requires that users download
#' the fitted `neuroFX2021` neural networks published by Rudiyanto et al. 2021
#' in Supplementary Material 1 (resulting in a local file named `xxx_mmc1.zip`).
#' This needs to be unzipped and the resulting `tar` file unpacked;
#' this produces a folder `R code for neuroFX2021`.
#' This folder contains two R data files : `ssc.RData` and `sscbd.RData`.
#' The argument `file_neuroFX2021` is the file name (with path) to `sscbd.RData`
#' if soil density data are available and to `ssc.RData` otherwise
#' (see Rudiyanto et al. 2021).
#' The path to the appropriate file can be set per R session
#' via option `"RSW2_FILENEUROFX2021"`
#' (and avoid passing it directly as argument to the function);
#' this can be useful, for example, if `ptf_estimate()` is used for `FXW`.
#'
#' @md
ptf_neuroFX2021_for_FXW <- function(
sand,
clay,
bdensity = NULL,
file_neuroFX2021 = getOption("RSW2_FILENEUROFX2021", NULL),
...
) {
stopifnot(ptf_neuroFX2021_availability(file_neuroFX2021))
# Load `neuroFX2021`
nfx <- new.env()
load(file_neuroFX2021, envir = nfx)
# Check whether type of neuroFX2021 is
# SSC (sand, silt, clay) or
# SSCBD (sand, silt, clay, bulk density)
is_sscbd <- dim(nfx[["tW1"]])[2] == 4
# Prepare soil data
tmp_txt <- data.frame(
sand = sand,
silt = 1 - (sand + clay),
clay = clay
)
if (!is.null(bdensity)) {
if (is_sscbd) {
tmp_txt[, "bd"] <- bdensity
} else {
warning(
"`ptf_neuroFX2021_for_FXW()`: ",
"`bdensity` ignored because ",
"'neuroFX2021' object is for SSC (sand, silt, clay)."
)
}
} else {
if (is_sscbd) {
stop(
"`ptf_neuroFX2021_for_FXW()`: ",
"'neuroFX2021' object is for SSCBD (sand, silt, clay, bulk density) ",
"but `bdensity` contains no values."
)
}
}
# Evaluate neuroFX2021
res <- iterate_neuroFX(tmp_txt, nfx, niter = dim(nfx[["tW1"]])[1])
# Aggregate across iterations
tmp_res <- lapply(c("mean", "sd"), function(f) apply(res, 2:3, f))
# Backtransformation
for (k in seq_along(tmp_res)) {
# backtransform log(alpha) -> alpha
tmp_res[[k]][, 2] <- exp(tmp_res[[k]][, 2])
# backtransform log(n - 1) -> n
tmp_res[[k]][, 3] <- 1 + exp(tmp_res[[k]][, 3])
# backtransform log10(Ks) -> Ks
tmp_res[[k]][, 5] <- 10 ^ tmp_res[[k]][, 5]
}
list(
mean = tmp_res[[1]],
sd = tmp_res[[2]]
)
}
# Checks availability of `neuroFX2021` `PTF`
#
# Note: `check_ptf_availability()` requires function name
# to match pattern "ptf_XXX_availability" where XXX = name of PTF
ptf_neuroFX2021_availability <- function(
file_neuroFX2021 = getOption("RSW2_FILENEUROFX2021", NULL),
verbose = interactive(),
...
) {
res <- !is.null(file_neuroFX2021) && file.exists(file_neuroFX2021)
if (!res && verbose) {
message(
"`ptf_neuroFX2021_availability()`: ",
"data file 'file_neuroFX2021' does not exist; ",
"see documentation for `ptf_neuroFX2021_for_FXW()` and consider setting ",
"`options(RSW2_FILENEUROFX2021 = \"path/to/sscbd.RData\")`"
)
}
res
}
# Evaluate neural net: code based on Rudiyanto et al. 2021
eval_nnet <- function(X, W1, W2) {
N <- nrow(X)
# from input layer to hidden layer
xt <- rbind(t(X), matrix(1, nrow = 1, ncol = N))
h <- W1 %*% xt
# activation function
y1 <- tanh(h)
# from hidden layer to output layer
t(W2 %*% rbind(y1, matrix(1, nrow = 1, ncol = N)))
}
# Iterate over neuroFX2021: code based on Rudiyanto et al. 2021
iterate_neuroFX <- function(x, nfx, niter = 50) {
res <- array(dim = c(niter, nrow(x), 6))
# loop though n iteration
for (k in seq_len(niter)) {
# Predict theta_sat, log(alpha), log(n-1), m
res[k, , 1:4] <- eval_nnet(
x,
W1 = nfx[["tW1"]][k, , ],
W2 = nfx[["tW2"]][k, , ]
)
# Predict log10(K_sat), L
res[k, , 5:6] <- eval_nnet(
x,
W1 = nfx[["kW1"]][k, , ],
W2 = nfx[["kW2"]][k, , ]
)
}
res
}
#' Check Soil Water Retention Curve parameters
#'
#' @inheritParams SWRCs
#'
#' @section Notes:
#' The argument selecting `SWRC` (`swrc_name`) is recycled
#' for multiple parameter sets, i.e., rows of `swrcp`.
#'
#' @section Details:
#' [swrc_names()] lists implemented `SWRCs`.
#'
#' @seealso [ptf_estimate()]
#'
#' @examples
#' swrc_name <- "Campbell1974"
#' ptf_name <- "Cosby1984AndOthers"
#' swrcp <- ptf_estimate(
#' sand = c(0.5, 0.3),
#' clay = c(0.2, 0.1),
#' fcoarse = c(0, 0),
#' swrc_name = swrc_name,
#' ptf_name = ptf_name
#' )
#'
#' check_swrcp(swrc_name, swrcp)
#' check_swrcp(swrc_name, swrcp[1, ])
#'
#' swrcp2 <- swrcp
#' swrcp2[1, 1] <- -10
#' check_swrcp(swrc_name, swrcp2)
#'
#' @export
#' @md
check_swrcp <- function(swrc_name, swrcp) {
# lengths of arguments are checked by `C_rSW2_SWRC_check_parameters()`
.Call(
C_rSW2_SWRC_check_parameters,
swrc_type = rep_len(
encode_name2swrc(swrc_name)[1],
if (is.matrix(swrcp)) nrow(swrcp) else 1
),
swrcp = swrcp
)
}
#------ Soil Water Retention Curves ------
#' Conversion between bulk soil water content and soil water potential
#'
#' @inheritParams SWRCs
#' @param direction A character string. Indicates the direction of
#' soil water conversion.
#' @param x A numeric value, vector, or matrix.
#' The soil water values to be converted,
#' either soil water potential (units `[MPa]`) of the soil matric component or
#' bulk volumetric water content (units `[cm/cm]`).
#' @param outer_if_equalsize A logical value.
#' Relevant only if `x` of length `l` and soils of length `d` are equal.
#' If `TRUE`, then the returned object has a size of `l x d` = `l x l`
#' where the `d` sets of soil values are repeated for each value of `x`.
#' If `FALSE` (default), then the returned object has a size of `l` = `d`
#' where the the `SWRC` conversion is applied to the
#' first element of `x` and soils, the second elements, and so on.
#'
#' @return The dimensions of the output are a function of `x` and the
#' number of soil values (e.g., rows or length of `swrc[["swrcp"]]`).
#' The returned object has:
#' \itemize{
#' \item length `l` if both `x` and soils are of length `l`.
#' \item length `l` if `x` has length `l` and there is one soil.
#' \item length `d` if `x` is one value and soils are of length `d`.
#' \item size `l x d` if `x` has length `l` and soils are of length `d`
#' (if `l` and `d` are not equal or `outer_if_equalsize` is `TRUE`;
#' cf. the first case);
#' the `d` sets of soil values are repeated for each value of `x`.
#' \item size `l x d` if `x` has size `l x d` and there is one soil.
#' the soil is repeated for each value of `x`.
#' \item size `l x d` if `x` has size `l x d` and soils are of length `d`
#' the `d` sets of soil values are repeated for each row of `x`.
#' }
#'
#'
#' @inherit SWRCs references
#'
#' @section Details:
#' [swrc_names()] lists implemented `SWRCs`;
#' [ptf_names()] lists implemented `PTFs`; and
#' [check_ptf_availability()] checks availability of `PTFs`.
#'
#' @section Details:
#' For backward compatibility, `fcoarse` and `layer_width` may be missing.
#' If they are missing, then the soils are assumed to contain
#' `0%` coarse fragments and be represented by `1 cm` wide soil layers.
#'
#' @section Details:
#' Arguments `sand`, `clay`, and `bdensity` are only required
#' if `SWRC` parameter values need to be estimated on the fly,
#' i.e., if `swrc` does not contain the element `swrcp`
#' (with suitable `SWRC` parameter values).
#' This is handled by [ptf_estimate()] and additionally requires
#' the element `ptf_name` for argument `swrc`.
#'
#' @section Details:
#' If `swrc` contains element `swrcp` with one set of `SWRC` parameters,
#' i.e., one row, then the parameter set is repeated for each value of `x`.
#'
#' @section Details:
#' If `vwc` inputs represent the matric component
#' (instead of expected bulk values), then set `fcoarse` to 0.
#' This works, however, only if `swrcp` are provided or `fcoarse` is not
#' utilized by the requested `PTF`.
#'
#'
#' @seealso
#' [ptf_estimate()],
#' [check_swrcp()],
#' [check_ptf_availability()]
#'
#' @examples
#' fsand <- c(0.5, 0.3)
#' fclay <- c(0.2, 0.1)
#' fcrs1 <- c(0, 0)
#' fcrs2 <- c(0.4, 0.1)
#'
#' swrc1 <- list(
#' name = "Campbell1974",
#' swrcp = ptf_estimate(
#' sand = fsand,
#' clay = fclay,
#' fcoarse = fcrs1,
#' swrc_name = "Campbell1974",
#' ptf_name = "Cosby1984"
#' )
#' )
#' swrc_swp_to_vwc(-1.5, fcoarse = fcrs1, swrc = swrc1)
#' swrc_swp_to_vwc(c(-1.5, NA), fcoarse = fcrs1, swrc = swrc1)
#' swrc_swp_to_vwc(-1.5, fcoarse = fcrs1, sand = fsand, clay = fclay)
#' swrc_vwc_to_swp(c(0.10, 0.15, 0.20), fcoarse = fcrs1, swrc = swrc1)
#' swrc_vwc_to_swp(c(0.10, NA, 0.20), fcoarse = fcrs1, swrc = swrc1)
#'
#' swrc2 <- list(
#' name = "Campbell1974",
#' swrcp = ptf_estimate(
#' sand = fsand,
#' clay = fclay,
#' fcoarse = fcrs2,
#' swrc_name = "Campbell1974",
#' ptf_name = "Cosby1984"
#' )
#' )
#' swrc_swp_to_vwc(-1.5, fcoarse = fcrs2, swrc = swrc2)
#' (1 - fcrs2) * swrc_swp_to_vwc(-1.5, swrc = swrc2)
#' swrc_swp_to_vwc(-1.5, fcoarse = fcrs2, sand = fsand, clay = fclay)
#' swrc_vwc_to_swp(c(0.10, 0.15, 0.20), fcoarse = fcrs2, swrc = swrc2)
#'
#'
#' # Available water holding capacity "AWC"
#' soils <- swSoils_Layers(rSOILWAT2::sw_exampleData)
#' p <- ptf_estimate(
#' sand = soils[, "sand_frac"],
#' clay = soils[, "clay_frac"],
#' fcoarse = soils[, "gravel_content"]
#' )
#' tmp <- swrc_swp_to_vwc(
#' c(-1.5, -0.033),
#' fcoarse = soils[, "gravel_content"],
#' swrc = list(name = "Campbell1974", swrcp = p)
#' )
#' awc <- diff(c(0, soils[, "depth_cm"])) * as.vector(diff(tmp))
#'
#'
#' # Shape of SWRCs
#' theta <- seq(0.05, 0.55, by = 0.001)
#' soils <- data.frame(
#' sand_frac = c(sand = 0.92, silty_loam = 0.17, silty_clay = 0.06),
#' clay_frac = c(0.03, 0.13, 0.58),
#' bd = c(1.614, 1.464, 1.437)
#' )
#' phi <- list(
#' Campbell1974 = swrc_vwc_to_swp(
#' theta,
#' sand = soils[, "sand_frac"],
#' clay = soils[, "clay_frac"],
#' swrc = list(swrc_name = "Campbell1974", ptf_name = "Cosby1984")
#' )
#' )
#'
#' if (check_ptf_availability("Rosetta3")) {
#' phi[["vanGenuchten1980"]] <- swrc_vwc_to_swp(
#' theta,
#' sand = soils[, "sand_frac"],
#' clay = soils[, "clay_frac"],
#' bdensity = soils[, "bd"],
#' swrc = list(swrc_name = "vanGenuchten1980", ptf_name = "Rosetta3")
#' )
#' }
#'
#' # Use PTF "neuroFX2021" to estimate parameters of SWRC `FXW`
#' \dontrun{
#' # Set neuroFX2021 file path, see details in `ptf_neuroFX2021_for_FXW()`
#' options(RSW2_FILENEUROFX2021 = "path/to/sscbd.RData")
#' }
#'
#' if (check_ptf_availability("neuroFX2021")) {
#' phi[["FXW"]] <- swrc_vwc_to_swp(
#' theta,
#' sand = soils[, "sand_frac"],
#' clay = soils[, "clay_frac"],
#' bdensity = soils[, "bd"],
#' swrc = list(swrc_name = "FXW", ptf_name = "neuroFX2021")
#' )
#' }
#'
#' if (interactive() && requireNamespace("graphics")) {
#' par_prev <- graphics::par(mfcol = c(length(phi), 1))
#'
#' for (k in seq_along(phi)) {
#' graphics::matplot(
#' theta, -phi[[k]],
#' type = "l",
#' log = "y",
#' xlim = c(0, max(theta)),
#' xlab = "theta [m/m]",
#' ylim = c(1e-4, 1e6),
#' ylab = "-phi [MPa]",
#' main = paste0("Soil Water Retention Curve (", names(phi)[k], ")")
#' )
#' graphics::abline(h = -c(-1.5, -0.033), col = "gray", lty = 3)
#' graphics::legend("topright", rownames(soils), col = 1:3, lty = 1:3)
#' }
#'
#' graphics::par(par_prev)
#' }
#'
#'
#' @export
#' @md
swrc_conversion <- function(
direction = c("swp_to_vwc", "vwc_to_swp"),
x,
fcoarse,
layer_width,
swrc,
sand = NULL,
clay = NULL,
bdensity = NULL,
outer_if_equalsize = FALSE,
verbose = FALSE
) {
#--- Check inputs
direction <- match.arg(direction)
# `name` can be used as short form of `swrc_name`
if (!("swrc_name" %in% names(swrc)) && "name" %in% names(swrc)) {
swrc[["swrc_name"]] <- swrc[["name"]]
}
stopifnot("swrc_name" %in% names(swrc))
swrc[["swrc_name"]] <- std_swrc(swrc[["swrc_name"]])[1]
swrc[["swrc_type"]] <- encode_name2swrc(swrc[["swrc_name"]])
# Do we need to estimate swrcp?
swrc[["swrcp"]] <- if (
"swrcp" %in% names(swrc) && !is.null(swrc[["swrcp"]])
) {
if (is.null(dim(swrc[["swrcp"]]))) {
matrix(swrc[["swrcp"]], nrow = 1)
} else {
as.matrix(swrc[["swrcp"]])
}
}
# Do we have sufficient information to estimate swrcp?
if (
is.null(swrc[["swrcp"]]) &&
(
!all(c("swrc_name", "ptf_name") %in% names(swrc)) ||
is.null(sand) || is.null(clay)
)
) {
stop("Insufficient information to estimate SWRC parameters.")
}
# Do we have sufficient soil parameters?
if (missing(fcoarse) && missing(layer_width)) {
ntmp <- if (!is.null(swrc[["swrcp"]])) {
nrow(swrc[["swrcp"]])
} else {
if (!is.null(sand)) {
length(sand)
} else if (!is.null(clay)) {
length(clay)
}
}
if (!is.null(ntmp)) {
fcoarse <- rep(0, ntmp)
layer_width <- rep(1, ntmp)
} else {
stop("Insufficient soil parameters to use SWRC.")
}
} else if (missing(fcoarse)) {
fcoarse <- rep(0, length(layer_width))
} else if (missing(layer_width)) {
layer_width <- rep(1, length(fcoarse))
}
# Put together available soil parameters and check for consistency
soils <- list(
fcoarse = fcoarse,
layer_width = layer_width
)
if (is.null(swrc[["swrcp"]])) {
soils <- c(
soils,
list(sand = sand, clay = clay),
if (!is.null(bdensity)) list(bdensity = bdensity)
)
}
nsoils <- unique(lengths(soils))
if (length(nsoils) > 1) {
stop("Soil variables have different lengths.")
}
if (!is.null(swrc[["swrcp"]]) && nrow(swrc[["swrcp"]]) != nsoils) {
stop("Dimensions of `swrcp` and length of soil variables disagree.")
}
#--- Determine dimensions of data and result
nrx <- NROW(x)
ncx <- NCOL(x)
nx <- nrx * ncx
nx1d <- nrx == 1 || ncx == 1
res <- array(dim = c(nrx, ncx))
#--- Prepare inputs and make SWRC conversion
if (
nx1d && (nx == 1 || nsoils == 1 || (nx == nsoils && !outer_if_equalsize))
) {
# 1a. x [len = 1] + soils [len = 1] --> res [len = 1, dim = 1 x 1]
# nothing to prepare
if (nx == 1 && nsoils > 1) {
# 2. x [len = 1] + soils [len = d] --> res [len = d, dim = 1 x d]
x <- rep_len(x, nsoils)
} else if (nx > 1 && nx1d && nsoils == 1) {
# 3. x [len = l] + soils [len = 1] --> res [len = l, dim = l x 1]
soils <- lapply(soils, rep_len, length.out = nx)
if (!is.null(swrc[["swrcp"]])) {
swrc[["swrcp"]] <- swrc[["swrcp"]][rep(1, nx), , drop = FALSE]
}
} else if (nx == nsoils && !outer_if_equalsize) {
# 1b. x [len = l] + soils [len = l] --> res [len = l, dim = l x 1]
x <- as.vector(unlist(x))
}
if (is.null(swrc[["swrcp"]])) {
swrc[["swrcp"]] <- ptf_estimate(
sand = soils[["sand"]],
clay = soils[["clay"]],
fcoarse = soils[["fcoarse"]],
bdensity = soils[["bdensity"]],
swrc_name = swrc[["swrc_name"]],
ptf_name = swrc[["ptf_name"]]
)
}
res <- swrc_conversion_1d(direction, x, soils, swrc, verbose)
} else if (nx1d && nx > 1 && nsoils > 1) {
# 4. x [len = l] + soils [len = d] -> res [dim = l x d]
# (x repeated for each soil)
if (is.null(swrc[["swrcp"]])) {
swrc[["swrcp"]] <- ptf_estimate(
sand = soils[["sand"]],
clay = soils[["clay"]],
fcoarse = soils[["fcoarse"]],
bdensity = soils[["bdensity"]],
swrc_name = swrc[["swrc_name"]],
ptf_name = swrc[["ptf_name"]]
)
}
tmp <- lapply(
x,
function(v) {
swrc_conversion_1d(direction, rep_len(v, nsoils), soils, swrc, verbose)
}
)
res <- matrix(unlist(tmp), nrow = nx, ncol = nsoils, byrow = TRUE)
} else if (nx > 1 && !nx1d && nsoils == 1) {
# 5. x [dim = l x d] + soils [len = 1] --> res [dim = l x d]
soils <- lapply(soils, rep_len, length.out = nrx)
if (is.null(swrc[["swrcp"]])) {
swrc[["swrcp"]] <- ptf_estimate(
sand = soils[["sand"]],
clay = soils[["clay"]],
fcoarse = soils[["fcoarse"]],
bdensity = soils[["bdensity"]],
swrc_name = swrc[["swrc_name"]],
ptf_name = swrc[["ptf_name"]]
)
} else {
swrc[["swrcp"]] <- swrc[["swrcp"]][rep(1, nrx), , drop = FALSE]
}
res <- vapply(
seq_len(ncx),
function(k) swrc_conversion_1d(direction, x[, k], soils, swrc, verbose),
FUN.VALUE = rep(1, nrx),
USE.NAMES = FALSE
)
} else if (nx > 1 && !nx1d && nsoils == ncx) {
# 6. x [dim = l x d] + soils [len = d] --> res [dim = l x d]
# (soils repeated for row of x value)
if (is.null(swrc[["swrcp"]])) {
swrc[["swrcp"]] <- ptf_estimate(
sand = soils[["sand"]],
clay = soils[["clay"]],
fcoarse = soils[["fcoarse"]],
bdensity = soils[["bdensity"]],
swrc_name = swrc[["swrc_name"]],
ptf_name = swrc[["ptf_name"]]
)
}
swrc[["swrc_type"]] <- rep_len(swrc[["swrc_type"]], nsoils)
res <- vapply(
seq_len(ncx),
function(k) {
ids <- rep.int(k, nrx)
swrc_conversion_1d(
direction,
x = x[, k],
soils = lapply(soils, function(sp) sp[ids]),
swrc = list(
swrc_type = swrc[["swrc_type"]][ids],
swrcp = swrc[["swrcp"]][ids, , drop = FALSE]
),
verbose = verbose
)
},
FUN.VALUE = rep(1, nrx),
USE.NAMES = FALSE
)
} else {
stop("Unsuitable inputs.")
}
res
}
#' Helper function of \code{swrc_conversion} to access underlying C code
#' @noRd
swrc_conversion_1d <- function(direction, x, soils, swrc, verbose) {
prev_verbosity <- sw_verbosity(verbose = as.logical(verbose))
on.exit(sw_verbosity(prev_verbosity))
# lengths of arguments are checked by `C_rSW2_SWRC()`
nx <- length(x)
switch(
EXPR = direction,
# C_rSW2_SWRC(direction = 1) returns [cm] convert to [cm/cm]
swp_to_vwc = 1 / soils[["layer_width"]] * .Call(
C_rSW2_SWRC,
# x = SWP [MPa] convert to [-bar]
x = - 10 * x,
direction = 1L,
swrc_type = rep_len(swrc[["swrc_type"]], nx),
swrcp = swrc[["swrcp"]],
fcoarse = soils[["fcoarse"]],
width = soils[["layer_width"]]
),
# C_rSW2_SWRC(direction = 2) returns [-bar] convert to [MPa]
vwc_to_swp = - 0.1 * .Call(
C_rSW2_SWRC,
# x = VWC (bulk) [cm/cm] convert to SWC [cm]
x = x * soils[["layer_width"]],
direction = 2L,
swrc_type = rep_len(swrc[["swrc_type"]], nx),
swrcp = swrc[["swrcp"]],
fcoarse = soils[["fcoarse"]],
width = soils[["layer_width"]]
)
)
}
#' @describeIn swrc_conversion Convenience wrapper
#' to convert from `SWP` to bulk `VWC` with selected `SWRC`
#'
#' @param swp_MPa A numeric object. The soil water potential values
#' (units `[MPa]`) of the soil matric component to be converted to
#' bulk volumetric water content
#' (i.e., relative to the whole soil; units `[cm/cm]`).
#'
#' @export
#' @md
swrc_swp_to_vwc <- function(
swp_MPa,
fcoarse,
layer_width,
swrc = list(swrc_name = NULL, ptf_name = NULL, swrcp = NULL),
sand = NULL,
clay = NULL,
bdensity = NULL,
outer_if_equalsize = FALSE,
verbose = FALSE
) {
swrc_conversion(
direction = "swp_to_vwc",
x = swp_MPa,
sand = sand,
clay = clay,
fcoarse = fcoarse,
bdensity = bdensity,
layer_width = layer_width,
swrc = swrc,
outer_if_equalsize = outer_if_equalsize,
verbose = verbose
)
}
#' @describeIn swrc_conversion Convenience wrapper
#' to convert from bulk `VWC` to matric `SWP` with selected `SWRC`
#'
#' @param vwcBulk A numeric object. The volumetric water content values
#' (relative to the whole soil; units `[cm/cm]`)
#' to be converted to soil water potential (units `[MPa]`)
#' of the soil matric component.
#'
#' @export
#' @md
swrc_vwc_to_swp <- function(
vwcBulk,
fcoarse,
layer_width,
swrc = list(swrc_name = NULL, ptf_name = NULL, swrcp = NULL),
sand = NULL,
clay = NULL,
bdensity = NULL,
outer_if_equalsize = FALSE,
verbose = FALSE
) {
swrc_conversion(
direction = "vwc_to_swp",
x = vwcBulk,
sand = sand,
clay = clay,
fcoarse = fcoarse,
bdensity = bdensity,
layer_width = layer_width,
swrc = swrc,
outer_if_equalsize = outer_if_equalsize,
verbose = verbose
)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.