proper documentation, housecleaning

NAMESPACE

@@ -153,6 +153,7 @@ export(texture_to_taxpartsize)
 export(texture_to_texmod)
 export(thompson.bell.darkness)
 export(unroll)
+export(warpHorizons)
 export(xtableTauW)
 exportClasses(SoilProfileCollection)
 exportMethods("GHL<-")

R/warpHorizons.R

@@ -0,0 +1,167 @@
+
+
+# https://www.fao.org/3/cb0509en/CB0509EN.pdf
+
+
+#' @title Inflate / Deflate Horizon Thickness
+#' 
+#' @description This function applies a warping factor to the horizons of a single-profile `SoilProfileCollection` object. Warping values >1 will inflate horizon thickness, values <1 will deflate horizon thickness.
+#'
+#' @param x a `SoilProfileCollection` object
+#' @param fact numeric or character; warping factor specified as a single numeric value, vector of numeric values (length = nrow(x)), or column name of a horizon level attribute containing numeric values
+#' @param updateProfileID logical; modify profile IDs
+#' @param suffix character; suffix added to profile IDs when `updateProfileID = TRUE`
+#' 
+#' @author D.E. Beaudette and S.W. Salley
+#'
+#' @return a modified version of `x`, `SoilProfileCollection` object
+#' @export
+#'
+#' @examples 
+#' 
+#' # create an example profile
+#' s <- quickSPC('p1:AA|Bt1Bt1Bt1|Bt2Bt2B|Bt3|Cr|RRRRR')
+#' 
+#' # warp each horizon
+#' # values > 1: inflation
+#' # values < 1: deflation (erosion / compaction)
+#' s.w <- warpHorizons(s, fact = c(1.3, 0.7, 0.8, 1, 1, 1))
+#' 
+#' # combine original + warped
+#' x <- combine(s, s.w)
+#' 
+#' # compute profile bottom depths
+#' .bottoms <- x[, , .LAST, .BOTTOM]
+#' 
+#' # change in total depth after warping
+#' # used to vertically offset the warped profile
+#' .yoff <- c(0, .bottoms[1] - .bottoms[2])
+#' 
+#' # depths for line segments connecting horizon tops
+#' .y1 <- x[1, , .TOP]
+#' .y2 <- x[2, , .TOP] + .yoff[2]
+#' 
+#' # sketches
+#' # can't automatically add a depth axis
+#' par(mar = c(0.5, 0, 0, 2))
+#' plotSPC(
+#'   x, 
+#'   name.style = 'center-center', 
+#'   cex.names = 0.8, 
+#'   width = 0.2, 
+#'   max.depth = 150, 
+#'   depth.axis = list(line = -3), 
+#'   y.offset = .yoff
+#' )
+#' 
+#' # illustrate warping with arrows
+#' arrows(x0 = 1 + 0.25, y0 = .y1, x1 = 2 - 0.25, y1 = .y2, len = 0.1, col = 2)
+#' 
+#' # manually add depth axis
+#' axis(side = 4, line = -3.5, las = 1, at = seq(from = 0, to = 150, by = 25))
+#' 
+#' 
+#' # apply to multiple profiles
+#' # text-based template
+#' .template <- c(
+#' 'P1:AAA|BwBwBwBw|CCCCCCC|CdCdCdCd',
+#' 'P2:ApAp|AA|E|BhsBhs|Bw1Bw1|CCCCC',
+#' 'P3: A|Bt1Bt1Bt1|Bt2Bt2Bt2|Bt3|Cr|RRRRR'
+#' )
+#' 
+#' # each horizon label is '10' depth-units (default)
+#' s <- quickSPC(.template)
+#' 
+#' # random warping factor, by horizon
+#' s$f <- runif(n = nrow(s), min = 0.8, max = 1.2)
+#' 
+#' # warp horizons by profile, result is a list of SPCs
+#' s.w <- profileApply(s, FUN = warpHorizons, fact = 'f')
+#'
+#' # flatten list -> SoilProfileCollection
+#' s.w <- combine(s.w)
+#' 
+#' # combine with original SPC
+#' x <- combine(s, s.w)
+#' 
+#' # sketches 
+#' par(mar = c(0.5, 0, 0, 2.5))
+#' plotSPC(
+#'   x,
+#'   name.style = 'center-center',
+#'   cex.names = 0.8,
+#'   width = 0.3,
+#'   max.depth = 165,
+#'   depth.axis = list(line = -2)
+#' )
+#' 
+#'
+warpHorizons <- function(x, fact, updateProfileID = TRUE, suffix = '-w') {
+
+  ## TODO: 
+  ## * vectorize over profiles, and make more efficient
+  ## * support the more general, Ax + c affine transformation for 1D geometry
+  ##   -> https://www.cs.utexas.edu/users/fussell/courses/cs384g-fall2011/lectures/lecture07-Affine.pdf
+
+
+
+  # extract parts and pieces of the SPC
+  .h <- horizons(x)
+  .htb <- horizonDepths(x)
+  .n <- nrow(.h)
+
+  # sanity check: fact should be:
+  # * numeric, length == 1, used by all horizons
+  # * numeric, length == nrow(x), each horizon has its own factor
+  # * hz column name
+
+  if(inherits(fact, 'character')) {
+    fact <- x[[fact]]
+    if(is.null(fact)) {
+      stop('fact must name a column in horizons')
+    }
+  } else if(inherits(fact, 'numeric')) {
+
+    if(length(fact) > 1 && length(fact) != .n) {
+      stop('fact must be length 1 or nrow(x)')
+    }
+
+  } else {
+    stop('fact must be either character vector or numeric')
+  }
+
+  # warping is applied to horizon thickness
+  .thick <- .h[[.htb[2]]] - .h[[.htb[1]]]
+
+  # apply inflation/deflation factor to horizon thickness
+  # round to integers
+  .thick <- round(.thick * fact)
+
+  # future: apply offset here
+
+  # generate new horizon depth sequence
+  # starting from original topmost depth
+  .start <- .h[1,  .htb[1]]
+
+  # tops / bottoms
+  .tops <- c(.start, cumsum(.thick[-.n]))
+  .bottoms <- c(cumsum(.thick))
+
+  # replace original values
+  .h[[.htb[1]]] <- .tops
+  .h[[.htb[2]]] <- .bottoms
+
+  # re-pack horizons
+  replaceHorizons(x) <- .h
+
+  # optionally update profile ID
+  # this is handy when trying to combine(old, new)
+  if(updateProfileID) {
+    .pID <- profile_id(x)
+    profile_id(x) <- sprintf("%s%s", .pID, suffix)
+  }
+
+  return(x)
+}
+
+

misc/sandbox/warp-examples-testing.R

@@ -0,0 +1,81 @@
+library(aqp)
+library(soilDB)
+
+# function example
+
+# create an example profile
+s <- quickSPC('p1:AA|Bt1Bt1Bt1|Bt2Bt2B|Bt3|Cr|RRRRR')
+
+# warp each horizon
+# values > 1: inflation
+# values < 1: deflation (erosion / compaction)
+s.w <- warpHorizons(s, fact = c(1.3, 0.7, 0.8, 1, 1, 1))
+
+# combine original + warped
+x <- combine(s, s.w)
+
+# compute profile bottom depths
+.bottoms <- x[, , .LAST, .BOTTOM]
+
+# change in total depth after warping
+# used to vertically offset the warped profile
+.yoff <- c(0, .bottoms[1] - .bottoms[2])
+
+# depths for line segments connecting horizon tops
+.y1 <- x[1, , .TOP]
+.y2 <- x[2, , .TOP] + .yoff[2]
+
+# sketches
+# can't automatically add a depth axis
+par(mar = c(0.5, 0, 0, 2))
+plotSPC(
+  x, 
+  name.style = 'center-center', 
+  cex.names = 0.8, 
+  width = 0.2, 
+  max.depth = 150, 
+  depth.axis = list(line = -3), 
+  y.offset = .yoff
+)
+
+# illustrate warping with arrows
+arrows(x0 = 1 + 0.25, y0 = .y1, x1 = 2 - 0.25, y1 = .y2, len = 0.1, col = 2)
+
+# manually add depth axis
+axis(side = 4, line = -3.5, las = 1, at = seq(from = 0, to = 150, by = 25))
+
+
+## TODO: back-calculate 1D affine transform from A -> B, or from A -> collection
+# (s.w[, , .BOTTOM] - s.w[, , .TOP]) / (s[, , .BOTTOM] - s[, , .TOP])
+# 
+# lm(s.w[, , .BOTTOM] - s.w[, , .TOP] ~ s[, , .BOTTOM] - s[, , .TOP])
+
+
+
+
+## real data
+
+
+o <- fetchOSD('zook')
+oo <- warpHorizons(o, fact = c(1.8, 1.3, 0.6, 0.75, 0.8, 1, 1, 1))
+x <- combine(o, oo)
+
+.y1 <- x[1, , .TOP]
+.y2 <- x[2, , .TOP]
+
+par(mar = c(1, 0, 0 , 2))
+plotSPC(x, name.style = 'center-center', cex.names = 0.8, width = 0.2, max.depth = 200, depth.axis = list(line = -3))
+arrows(x0 = 1 + 0.25, y0 = .y1, x1 = 2 - 0.25, y1 = .y2, len = 0.1, col = 2)
+
+
+
+o$fact <- c(1, 1, 1, 1, 1, 1, 1, 1)
+oo$fact <- c(1.8, 1.3, 0.6, 0.75, 0.8, 1, 1, 1)
+x <- combine(o, oo)
+
+par(mar = c(1, 0, 3 , 1))
+plotSPC(x, name.style = 'center-center', cex.names = 0.8, width = 0.2, max.depth = 200, depth.axis = FALSE, hz.depths = TRUE, color = 'fact')
+arrows(x0 = 1 + 0.33, y0 = .y1, x1 = 2 - 0.22, y1 = .y2, len = 0.1, col = 'black')
+
+
+