Merge pull request #306 from ncss-tech/thicknessOf

Add `thicknessOf()`

NAMESPACE

@@ -151,6 +151,7 @@ export(texmod_to_fragvoltot)
 export(textureTriangleSummary)
 export(texture_to_taxpartsize)
 export(texture_to_texmod)
+export(thicknessOf)
 export(thompson.bell.darkness)
 export(unroll)
 export(warpHorizons)

R/thicknessOf.R

@@ -0,0 +1,116 @@
+#' Calculate Thickness of Horizons Matching Logical Criteria
+#'
+#' This function calculates the cumulative (`method="cumulative"`, default) or maximum difference between (`method="minmax"`) horizons within a profile that match a defined pattern (`pattern`) or, more generally, any set of horizon-level logical expressions encoded in a function (`FUN`).
+#' 
+#' @param x A _SoilProfileCollection_
+#' @param pattern _character_. A pattern to match in `hzdesgn`; used with the default `FUN` definition for regular expression pattern matching on horizons.
+#' @param hzdesgn _character_. A column containing horizon designations or other horizon-level character label used to identify matches; used with the default `FUN` definition.
+#' @param method _character_. Either `"cumulative"` (default) or `"minmax"`. See details.
+#' @param prefix _character_. Column prefix for calculated `thickvar` (and `depthvar` for `method="minmax"`) column results. Default: `""`.
+#' @param thickvar _character_ Length `1`. Column name to use for calculated thickness column. Default: `"thickness"`
+#' @param depthvars _character_. Length `2`. Column names to use for calculated minimum top depth and maximum bottom depth in `method="minmax"`. Default: `horizonDepths(x)`
+#' @param FUN _function_. A function that returns a _logical_ vector equal in length to the number of horizons in `x`. See details.
+#' @param na.rm _logical_. Omit `NA` values in summaries of thickness and in matching? Default: `FALSE`
+#' @param ... Additional arguments passed to the matching function `FUN`.
+#'
+#' @return A _data.frame_-like object (corresponding to the `aqp_df_class()` of `x`) with one row per profile in `x`. First column is always the profile ID which is followed by `"thickness"`. In `method="minmax"` the upper and lower boundaries used to calculate `"thickness"` are also returned as `"tmin"` and `"tmax"` columns, respectively.
+#' 
+#' @details
+#' The two thickness methods currently available are:
+#'  - `method="cumulative"` (default): cumulative thickness of horizons where `FUN` returns true
+#'  - `method="minmax"`: maximum bottom depth minus minimum top depth of horizons where `FUN` returns true
+#' 
+#' If a custom function (`FUN`) is used, it should accept arbitrary additional arguments via an ellipsis (`...`).
+#' It is not necessary to do anything with arguments, but the result should match the number of horizons found in the input SoilProfileCollection `x`.
+#' 
+#' @export
+#'
+#' @examples
+#' data("jacobs2000")
+#' 
+#' # cumulative thickness of horizon designations matching "Bt"
+#' thicknessOf(jacobs2000, "Bt")
+#' 
+#' # maximum bottom depth minus minimum top depth of horizon designations matching "Bt"
+#' thicknessOf(jacobs2000, "Bt", prefix = "Bt_", method = "minmax")
+#' 
+#' # cumulative thickness of horizon designations matching "A|B"
+#' thicknessOf(jacobs2000, "A|B", prefix = "AorB_")
+#' 
+#' # maximum bottom depth minus minimum top depth of horizon designations matching "A|B"
+#' thicknessOf(jacobs2000, "A|B", method = "minmax", prefix = "AorB_")
+#' # note that the latter includes the thickness of E horizons between the A and the B
+#' 
+#' # when using a custom function (be sure to accept ... and consider the effect of NA values)
+#' 
+#' # calculate cumulative thickness of horizons containing >18% clay
+#' thicknessOf(jacobs2000, prefix = "claygt18_", 
+#'             FUN = function(x, ...) !is.na(x[["clay"]]) & x[["clay"]] > 18)
+#' 
+thicknessOf <- function(x, 
+                        pattern = NULL, 
+                        hzdesgn = hzdesgnname(x, required = TRUE), 
+                        method = "cumulative",
+                        prefix = "",
+                        thickvar = "thickness",
+                        depthvars = horizonDepths(x),
+                        FUN = function(x, pattern, hzdesgn, ...) grepl(pattern, x[[hzdesgn]]),
+                        na.rm = FALSE,
+                        ...) {
+  .internalTHK <- NULL
+  .interalHZM <- NULL
+
+  if (is.null(hzdesgn) || !hzdesgn %in% horizonNames(x)) {
+    stop("Horizon designation column (", hzdesgn, ") does not exist.")
+  }
+
+  # check inputs
+  method <- match.arg(tolower(gsub("/", "", method)), c("cumulative", "minmax"))
+  stopifnot(length(thickvar) == 1)
+  stopifnot(length(depthvars) == 2)
+
+  # extract SPC column names and horizon data
+  hzd <- horizonDepths(x)
+  idn <- idname(x)
+  h <- data.table::data.table(horizons(x))
+
+  # determine horizons matching criteria of interest
+  h$.internalHZM <- FUN(x, pattern = pattern, hzdesgn = hzdesgn, na.rm = na.rm, ...)
+
+  # create a named list for data.table aggregation 
+  lid <- list(h[[idn]])
+  names(lid) <- idn
+
+  if (method == "cumulative") {
+
+    # sum thicknesses of all matching horizons
+    h$.internalTHK <- x[[hzd[2]]] - x[[hzd[1]]]
+    res <- h[, list(thickness = sum(.internalTHK[.internalHZM], na.rm = na.rm)), by = lid]
+    colnames(res)[2] <- thickvar
+
+  } else if (method == "minmax") {
+
+    # determine minimum top depth and maximum bottom depth of matching horizons
+    res <- h[, list(`min` = suppressWarnings(min(.SD[[1]][.internalHZM], na.rm = na.rm)),
+                    `max` = suppressWarnings(max(.SD[[2]][.internalHZM], na.rm = na.rm))),
+              by = lid, .SDcols = c(hzd, ".internalHZM")]
+
+    # calculate thickness as MAX(bottom) - MIN(top)
+    res$thickness <- res$`max` - res$`min`
+
+    # if a profile has no matching horizons min/max results will be +/- infinity 
+    # this means the profile has 0 thickness of matching horizons
+    res$thickness[!is.finite(res$thickness)] <- 0L
+
+    # use user-defined labels
+    colnames(res)[2:3] <- depthvars
+    colnames(res)[4] <- thickvar
+  } else stop("unknown thickness method: ", shQuote(method), call. = FALSE)
+
+  # add prefix
+  colnames(res)[2:ncol(res)] <- paste0(prefix, colnames(res)[2:ncol(res)])
+
+  # return as SPC data.frame class type
+  .as.data.frame.aqp(res, aqp_df_class(x))
+}
+

tests/testthat/test-thicknessOf.R

@@ -0,0 +1,44 @@
+test_that("thicknessOf works", {
+  data("jacobs2000")
+
+  ## method="cumulative"
+
+  # cumulative thickness of horizon designations matching "A|B"
+  x1 <- thicknessOf(jacobs2000, "A|B", prefix = "AorB_")
+  expect_equal(nrow(x1), length(jacobs2000))
+  expect_equal(x1$AorB_thickness, c(131, 117, 136, 20, 54, 110, 43))
+
+  ## method="minmax"
+
+  # maximum bottom depth minus minimum top depth of horizon designations matching "A|B"
+  x2 <- thicknessOf(jacobs2000, "A|B", method = "minmax", prefix = "AorB_")
+  expect_equal(ncol(x2), 4)
+  expect_equal(x2$AorB_top, rep(0, nrow(x2)))
+  expect_equal(x2$AorB_thickness, c(156, 145, 175, 20, 135, 168, 140))
+  expect_true(all(x2$AorB_thickness >= x1$AorB_thickness))
+
+  ## custom logical function
+
+  # calculate cumulative thickness of horizons containing >18% clay
+  x3 <- thicknessOf(jacobs2000, FUN = function(x, ...) !is.na(x[["clay"]]) & x[["clay"]] > 18)
+  expect_equal(x3$thickness, c(170, 167, 81, 0, 0, 49, 0))
+
+  ## missing property and or depth data
+
+  # function without NA handling, and na.rm=FALSE
+  x4 <- thicknessOf(jacobs2000, FUN = function(x, ...) x[["clay"]] > 18)
+  expect_equal(x4$thickness, c(170, 167, 81, 0, NA_integer_, 49, 0))
+
+  # function without NA handling, and na.rm=TRUE
+  x5 <- thicknessOf(jacobs2000, FUN = function(x, ...) x[["clay"]] > 18, na.rm = TRUE)
+  expect_equal(x5$thickness, c(170, 167, 81, 0, 0, 49, 0))
+
+  # missing horizon depths, and na.rm=FALSE
+  jacobs2000@horizons$top[1] <- NA_integer_
+  x6 <- thicknessOf(jacobs2000, "A|B")
+  expect_equal(x6$thickness, c(NA_integer_, 117, 136, 20, 54, 110, 43))
+
+  # missing horizond depths, and na.rm = TRUE
+  x7 <- thicknessOf(jacobs2000, "A|B", na.rm = TRUE)
+  expect_equal(x7$thickness, c(113, 117, 136, 20, 54, 110, 43))
+})