Skip to content

Commit

Permalink
Merge pull request #211 from ncss-tech/estimateAWC
Browse files Browse the repository at this point in the history 
estimateAWC and correctAWC using RO2 component lookuptable
  • Loading branch information
@dylanbeaudette
dylanbeaudette authored Mar 3, 2021
2 parents e2f48bc + 943d456 commit 3e17a22
Show file tree
Hide file tree
Showing 4 changed files with 282 additions and 0 deletions.
195 changes: 195 additions & 0 deletions R/estimateAWC.R
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,195 @@
#' Estimate available water capacity for fine-earth fraction
#'
#' @param texcl character, USDA textural class fine earth fraction
#' @param omcl integer, Organic matter class. 1: less than 1.5 percent, 2: less than 5, 3: greate than 5
#' @param precision integer, Number of decimal places in result default: 2
#' @param FUN Function for interpolating between table values default: \code{mean}
#' @param ... Additional arguments to \code{FUN}
#'
#' @return A numeric vector double containing estimated available water capacities for fine-earth fraction.
#' @export estimateAWC
#'
#' @examples
#'
#' # organic matter, loam texture, low medium and high OM
#' base.awc <- estimateAWC(c("l","l","l"), c(1, 2, 3), na.rm = TRUE)
#' base.awc
#'
estimateAWC <- function(texcl, omcl, precision = 2,
FUN = mean, ...) {

#awc_lut, loaded from region2 textureclass-awc lookup table
awc_lut <- structure(list(texcl = structure(c(4L, 8L, 16L, 6L, 24L, 11L, 13L, 12L, 14L, 5L, 22L, 7L,
25L, 10L, 21L, 20L, 18L, 3L, 19L, 17L, 1L, 2L, 9L, 15L, 23L),
.Label = c("C", "C_heavy", "CL", "COS", "COSL", "FS", "FSL", "GR", "HPM", "L", "LCOS", "LFS", "LS",
"LVFS", "MPM", "S", "SC", "SCL", "SIC", "SICL", "SIL", "SL", "SPM", "VFS", "VFSL"), class = "factor"),
loOM_l = c(0.02, 0.02, 0.05, 0.05, 0.05, 0.05, 0.06, 0.08, 0.09, 0.09, 0.1, 0.13, 0.14, 0.17, 0.15, 0.17, 0.14, 0.17, 0.14, 0.14, 0.14, 0.12, NA, NA, NA),
loOM_h = c(0.04, 0.04, 0.08, 0.08, 0.08, 0.07, 0.08, 0.11, 0.11, 0.12, 0.13, 0.15, 0.17, 0.18, 0.2, 0.21, 0.18, 0.21, 0.17, 0.16, 0.16, 0.15, NA, NA, NA),
mdOM_l = c(0.03, 0.03, 0.06, 0.06, 0.06, 0.06, 0.07, 0.09, 0.1, 0.1, 0.12, 0.15, 0.16, 0.17, 0.2, 0.18, 0.16, 0.15, 0.11, 0.15, 0.15, 0.13, NA, NA, NA),
mdOM_h = c(0.05, 0.05, 0.08, 0.08, 0.08, 0.07, 0.09, 0.12, 0.12, 0.13, 0.14, 0.17, 0.19, 0.19, 0.22, 0.2, 0.18, 0.19, 0.13, 0.17, 0.17, 0.16, NA, NA, NA),
hiOM_l = c(0.04, 0.04, 0.07, 0.07, 0.07, 0.07, 0.08, 0.1, 0.11, 0.11, 0.13, 0.16, 0.18, 0.2, 0.22, 0.21, 0.18, 0.17, 0.12, 0.13, 0.13, 0.15, 0.35, 0.45, 0.55),
hiOM_h = c(0.06, 0.06, 0.09, 0.09, 0.09, 0.08, 0.1, 0.13, 0.13, 0.14, 0.15, 0.18, 0.2, 0.22, 0.24, 0.23, 0.2, 0.19, 0.14, 0.15, 0.15, 0.18, 0.45, 0.55, 0.65)),
class = "data.frame", row.names = c(NA, -25L))

if (length(texcl) != length(omcl))
stop("Error: Mismatch in length of input vectors `texcl` and `omcl`")

split.lut <- split(awc_lut, f = awc_lut$texcl)
out <- rep(NA, length(texcl))

for (i in 1:length(texcl)) {
res <- lapply(split.lut, function(x) {
if (grepl(pattern = sprintf("^%s$", texcl[i]),
x$texcl, ignore.case = TRUE)) {
idx.lo <- 2*omcl[i]
idx.hi <- idx.lo + 1
if (!is.na(idx.hi) & !is.na(idx.lo)) {
return(x[,idx.lo:idx.hi])
} else return(NA)
}
})
out[i] <- FUN(unlist(res), ...)
}
return(round(out, precision))
}


#' Apply rock fragment or salt correction to available water content
#'
#' @param awc Numeric vector of available water capacities (e.g. from \code{estimateAWC})
#' @param total_rf Numeric vector of rock fragment volume percentage, 0 - 100
#' @param gravel Numeric vector of gravel volume percentage, 0 - 100
#' @param ec Numeric vector of electrical conductivity, mmhos/cm
#'
#' @return A numeric vector (double) containing estimated available water capacities corrected for rock fragments and salts
#' @export correctAWC
#'
#' @examples
#'
#' # medium organic matter, loam texture
#' base.awc <- 0.18 # estimateAWC(texcl = "l", omcl = 2, na.rm = TRUE)
#'
#' # medium organic matter, loam texture w/ 23% rock fragments by volume
#' corrected.awc <- correctAWC(base.awc, total_rf = 23)
#' corrected.awc
#'
#' # medium organic matter, loam texture w/ 0% frags by volume and 8 mmhos/cm salts
#' salty.awc <- correctAWC(base.awc, total_rf = 0, ec = 8)
#' salty.awc
correctAWC <- function(awc, total_rf, gravel = NULL, ec = NULL) {

# rf_lut from region2 rock fragment/salt correction lookup table
rf_lut <- .get_rf_ec_lut_RO2()

if (length(awc) != length(total_rf)) {
if (length(awc) == 1)
awc <- rep(awc, length(total_rf))
else if (length(total_rf) == 1)
total_rf <- rep(total_rf, length(awc))
else
stop("Error: Mismatch in length of input vectors.")
}
use_salt_correction <- FALSE

lut.row.idx <- .get_rf_lut_row(rf_lut, total = total_rf, gravel = gravel)
lut.col.idx <- .get_rf_lut_column(rf_lut, awc)

awc_c <- length(lut.row.idx)
awc_c[is.na(lut.col.idx)] <- NA

if (length(lut.row.idx) != length(lut.col.idx)) {
stop("unable to lookup rock fragment correction", call. = FALSE)
}

for (i in 1:length(lut.row.idx)) {
if (!is.na(lut.col.idx[i])) {
awc_c[i] <- rf_lut[lut.row.idx[i], lut.col.idx[i]]
} else {
awc_c[i] <- NA
}
}

if (!is.null(ec) & (length(ec) == 1 | length(ec) == length(awc)) & any(ec >= 2))
use_salt_correction = TRUE

if (length(ec) == 1)
ec <- rep(ec, length(awc))

if (use_salt_correction) {
# use EC to identify degree of reduction
lut.row.idx2 <- .get_rf_lut_row(rf_lut, ec = ec)

# use RF-reduced AWC for salt correction
lut.col.idx2 <- .get_rf_lut_column(rf_lut, awc_c)

awc_c2 <- length(lut.row.idx2)
awc_c2[is.na(lut.col.idx2)] <- NA

# use lookup table and return result
for (i in 1:length(lut.row.idx2)) {
if (!is.na(lut.col.idx2[i])) {
awc_c2[i] <- rf_lut[lut.row.idx2[i], na.omit(lut.col.idx2[i])]
} else {
awc_c2[i] <- NA
}
}
return(unlist(awc_c2))
}
return(unlist(awc_c))
}

# returns most limiting RF-reduction based on total RF, gravel or electrical conductivity
.get_rf_lut_row <- function(rf_lut, total = NULL, gravel = NULL, ec = NULL) {
n <- max(length(total), length(gravel), length(ec))
buf <- numeric(n)
for (i in 1:n) {
idx1 <- idx2 <- idx3 <- -2
if (!is.null(total[i]))
idx1 <- which(rf_lut$rf_byvol_l <= total[i] & rf_lut$rf_byvol_h > total[i])
if (!is.null(gravel[i]))
idx2 <- which(rf_lut$gravel_l <= gravel[i] & rf_lut$gravel_h > gravel[i])
if (!is.null(ec[i]))
idx3 <- which(rf_lut$salts_mmhos_l <= ec[i] & rf_lut$salts_mmhos_h > ec[i])
buf[i] <- max(idx1, idx2, idx3) + 1
}
return(buf)
}

.get_rf_lut_column <- function(rf_lut, awc_) {
# TODO: relies on 1:1 match of AWC to column in LUT - in future, use ranges and interpolate
return(match(round(awc_,2), rf_lut[1,7:ncol(rf_lut)], nomatch = NA) + 7)
}

.get_rf_ec_lut_RO2 <- function() {
# In this form the relationship is apparent
structure(list( rf_byvol_l = c( 0, 5, 10, 20, 30, 40, 45, 50, 55, 60, 65, 70, 75, 80, 80, 82.5, 85, 87.5, 90, 95, 100),
rf_byvol_h = c( 5, 10, 20, 30, 40, 45, 50, 55, 60, 65, 70, 75, 80, 80, 82.5, 85, 87.5, 90, 95, 99, 100),
salts_mmhos_l = c( 0, 2, 4, 6, 8, 10, 12, 14, 16, 17, 18, 19, 20, 85/4, 22.5, 95/4, 25, 27.5, 30, 35, 37.5),
salts_mmhos_h = c( 2, 4, 6, 8, 10, 12, 14, 16, 17, 18, 19, 20, 85/4, 22.5, 95/4, 25, 27.5, 30, 35, 37.5, 40),
gravel_l = c( 0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100),
gravel_h = c( 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 99, 100),
c1 = c(0.03, 0.03, 0.03, 0.03, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0, 0, 0),
c2 = c(0.04, 0.04, 0.04, 0.03, 0.03, 0.03, 0.03, 0.03, 0.02, 0.02, 0.02, 0.02, 0.02, 0.01, 0.01, 0.01, 0.01, 0.01, 0, 0, 0),
c3 = c(0.05, 0.05, 0.05, 0.04, 0.04, 0.04, 0.04, 0.03, 0.03, 0.03, 0.03, 0.02, 0.02, 0.02, 0.02, 0.01, 0.01, 0.01, 0.01, 0, 0),
c4 = c(0.06, 0.06, 0.05, 0.05, 0.05, 0.05, 0.04, 0.04, 0.04, 0.03, 0.03, 0.03, 0.03, 0.02, 0.02, 0.02, 0.01, 0.01, 0.01, 0, 0),
c5 = c(0.07, 0.07, 0.06, 0.06, 0.06, 0.05, 0.05, 0.05, 0.04, 0.04, 0.04, 0.03, 0.03, 0.02, 0.02, 0.02, 0.01, 0.01, 0.01, 0, 0),
c6 = c(0.08, 0.08, 0.07, 0.07, 0.06, 0.06, 0.06, 0.05, 0.05, 0.04, 0.04, 0.04, 0.03, 0.03, 0.02, 0.02, 0.02, 0.01, 0.01, 0, 0), c7 = c(0.09, 0.09, 0.08, 0.08, 0.07, 0.07, 0.06, 0.06, 0.05, 0.05, 0.05, 0.04, 0.04, 0.03, 0.03, 0.02, 0.02, 0.01, 0.01, 0.01, 0),
c8 = c(0.1, 0.1, 0.09, 0.09, 0.08, 0.08, 0.07, 0.07, 0.06, 0.06, 0.05, 0.05, 0.04, 0.04, 0.03, 0.03, 0.02, 0.02, 0.01, 0.01, 0),
c9 = c(0.11, 0.11, 0.1, 0.09, 0.09, 0.08, 0.08, 0.07, 0.07, 0.06, 0.06, 0.05, 0.04, 0.04, 0.03, 0.03, 0.02, 0.02, 0.01, 0.01, 0),
c10 = c(0.12, 0.11, 0.11, 0.1, 0.1, 0.09, 0.08, 0.08, 0.07, 0.07, 0.06, 0.05, 0.05, 0.04, 0.04, 0.03, 0.02, 0.02, 0.01, 0.01, 0),
c11 = c(0.13, 0.12, 0.12, 0.11, 0.1, 0.1, 0.09, 0.09, 0.08, 0.07, 0.07, 0.06, 0.05, 0.05, 0.04, 0.03, 0.03, 0.02, 0.01, 0.01, 0),
c12 = c(0.14, 0.13, 0.13, 0.12, 0.11, 0.11, 0.1, 0.09, 0.08, 0.08, 0.07, 0.06, 0.05, 0.05, 0.04, 0.04, 0.03, 0.02, 0.01, 0.01, 0),
c13 = c(0.15, 0.14, 0.14, 0.13, 0.12, 0.11, 0.1, 0.1, 0.09, 0.08, 0.08, 0.07, 0.05, 0.05, 0.05, 0.04, 0.03, 0.02, 0.02, 0.01, 0),
c14 = c(0.16, 0.15, 0.14, 0.14, 0.13, 0.12, 0.11, 0.1, 0.1, 0.09, 0.08, 0.07, 0.06, 0.06, 0.05, 0.04, 0.03, 0.02, 0.02, 0.01, 0),
c15 = c(0.17, 0.16, 0.15, 0.15, 0.14, 0.13, 0.12, 0.11, 0.1, 0.1, 0.09, 0.08, 0.07, 0.06, 0.05, 0.04, 0.03, 0.03, 0.02, 0.01, 0),
c16 = c(0.18, 0.17, 0.16, 0.15, 0.14, 0.14, 0.13, 0.12, 0.11, 0.1, 0.09, 0.08, 0.07, 0.06, 0.05, 0.05, 0.04, 0.03, 0.02, 0.01, 0),
c17 = c(0.19, 0.18, 0.17, 0.16, 0.15, 0.14, 0.13, 0.12, 0.11, 0.11, 0.1, 0.09, 0.08, 0.07, 0.06, 0.05, 0.04, 0.03, 0.02, 0.01, 0),
c18 = c(0.2, 0.19, 0.18, 0.17, 0.16, 0.15, 0.14, 0.13, 0.12, 0.11, 0.1, 0.09, 0.08, 0.07, 0.06, 0.05, 0.04, 0.03, 0.02, 0.01, 0),
c19 = c(0.21, 0.2, 0.19, 0.18, 0.17, 0.16, 0.15, 0.14, 0.13, 0.12, 0.11, 0.1, 0.09, 0.08, 0.06, 0.05, 0.04, 0.03, 0.02, 0.01, 0),
c20 = c(0.22, 0.21, 0.2, 0.19, 0.18, 0.17, 0.16, 0.15, 0.13, 0.13, 0.11, 0.11, 0.09, 0.08, 0.07, 0.05, 0.04, 0.03, 0.02, 0.01, 0),
c21 = c(0.23, 0.22, 0.21, 0.2, 0.19, 0.18, 0.17, 0.15, 0.14, 0.13, 0.12, 0.11, 0.09, 0.08, 0.07, 0.05, 0.04, 0.03, 0.02, 0.01, 0),
c22 = c(0.24, 0.23, 0.22, 0.2, 0.19, 0.19, 0.17, 0.16, 0.14, 0.14, 0.12, 0.11, 0.1, 0.09, 0.08, 0.06, 0.04, 0.03, 0.02, 0.01, 0),
c23 = c(0.25, 0.24, 0.23, 0.21, 0.2, 0.19, 0.18, 0.16, 0.15, 0.14, 0.13, 0.11, 0.1, 0.09, 0.08, 0.06, 0.05, 0.04, 0.03, 0.01, 0),
c24 = c(0.3, 0.29, 0.27, 0.26, 0.24, 0.23, 0.21, 0.2, 0.18, 0.17, 0.15, 0.14, 0.12, 0.11, 0.09, 0.08, 0.06, 0.05, 0.03, 0.02, 0)),
class = "data.frame", row.names = c(NA, -21L))
}
36 changes: 36 additions & 0 deletions man/correctAWC.Rd
View file

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

32 changes: 32 additions & 0 deletions man/estimateAWC.Rd
View file

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

19 changes: 19 additions & 0 deletions tests/testthat/test-awc.R
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,19 @@
context("estimateAWC and correctAWC - legacy AWC lookup tables")

test_that('basic functionality', {

# test range of organic matter, loam FEF textures
base.awc <- estimateAWC(c("l","l","l"), c(3, 2, 1), na.rm = TRUE)
expect_equal(base.awc, c(0.21, 0.18, 0.17))

# auto-recycling of singletons, 3 AWCs all corrected with 25% rf
corrected.awc <- correctAWC(base.awc, total_rf = 25)
expect_equal(corrected.awc, c(0.18, 0.15, 0.14))

# using gravel-size fragments often gives somewhat more restrictive AWC
corrected.awc2 <- correctAWC(base.awc, total_rf = 25, gravel = 25)
expect_equal(corrected.awc2, c(0.16, 0.15, 0.14))

salty.awc <- correctAWC(base.awc, total_rf = 0, ec = 8)
expect_equal(salty.awc, c(0.17, 0.14, 0.14))
})

0 comments on commit 3e17a22

Please sign in to comment.