testing related to overlapping horizons #296

NEWS.md

@@ -1,3 +1,6 @@
+# aqp 2.01 (2023-08-29)
+
+
 # aqp 2.0 (2023-08-28)
 This is a major update to aqp that may create some issues for code depending on specific inputs/outputs in aqp < 1.42, particularly those relying on `slice()`, `slab()`, and `profile_compare()`. `slice()` and `profile_compare()` are now deprecated, but will continue to work for the rest of calendar year 2023. There are no plans to maintain these functions beyond aqp 2.0. The new version of `slab()` is a drop-in replacement for the previous version of the function.
 

R/dice.R

@@ -1,9 +1,11 @@
 
-## https://github.com/ncss-tech/aqp/issues/115
+## details on slice() -> dice()
+# https://github.com/ncss-tech/aqp/issues/115
 
 ## TODO: 
 ##   * DT full outer join ideas
 ##     https://stackoverflow.com/questions/15170741/how-does-one-do-a-full-join-using-data-table
+##   * fully test / document compatibility with overlapping horizons: https://github.com/ncss-tech/aqp/issues/296
 
 setGeneric("dice", function(x,
                             fm = NULL,
@@ -93,6 +95,8 @@ setGeneric("dice", function(x,
     if (length(z) == 0) {
       z <- NULL
     } else {
+      ## TODO: adding new bogus horizons when there are overlapping source horizons
+
       # z-index is specified
       # note z-index defines slice tops, lower limit is (z + 1)
       # must fill from min(z) --- [gaps] --- max(z) + 1

R/fillHzGaps.R

@@ -5,6 +5,8 @@
 #' Gaps are defined as:
 #'  * within each profile, for horizons `i` to `n_hz`:
 #'  * `bottom_i != top_i+1 (but only to i = 1:(n_hz - 1)`
+#'  
+#'  @note This function cannot yet handle overlapping horizons.
 #'
 #' @param x `SoilProfileCollection` object
 #' 
@@ -77,6 +79,8 @@ fillHzGaps <- function(x, flag = TRUE, to_top = 0, to_bottom = max(x)) {
 
   h <- data.table::as.data.table(horizons(x))
 
+  ## TODO: adapt for use with overlapping horizons
+  # https://github.com/ncss-tech/aqp/issues/296
   lead.idx <- 2:nrow(h)
   lag.idx <- 1:(nrow(h) - 1)
 

misc/overlapping-horizons/notes.R

@@ -0,0 +1,96 @@
+library(aqp)
+
+## keeping track of progress here:
+# https://github.com/ncss-tech/aqp/issues/296
+
+
+
+## example data from Sharon
+
+# x <- read.csv('e:/temp/2020GRR014SPCdata.csv')
+# 
+# depths(x) <- upedonid ~ hzdept + hzdepb
+# hzdesgnname(x) <- 'hzname'
+# site(x) <- ~ NasisSiteName + taxonname + dspplotid + earthcovkind1 + earthcovkind2
+# 
+# x <- x[c(1, 2, 3, 8), ]
+# 
+# saveRDS(x, file = 'misc/overlapping-horizons/example-DSP-data.rds')
+
+
+# example DSP data
+x <- readRDS('misc/overlapping-horizons/example-DSP-data.rds')
+
+# all 4 profiles have overlap
+checkHzDepthLogic(x)
+
+## TODO: there should be an overlap error reported here:
+checkHzDepthLogic(x, byhz = TRUE)
+
+## consider a new function flagOverlap()
+
+# crude prototype, single profile at a time
+flagOverlap <- function(x) {
+
+  .fo <- function(i) {
+    .x <- i[, , .TOP]
+    .o <- overlapMetrics(.x, thresh = 0.1)
+
+    .res <- rep(FALSE, times = length(.x))
+    .res[.o$idx] <- TRUE
+
+    return(.res)
+  }
+
+   profileApply(x, .fo, simplify = TRUE)
+}
+
+x$.overlapFlag <- flagOverlap(x)
+
+par(mar = c(2, 0, 3, 2))
+options(.aqp.plotSPC.args = list(name.style = 'center-center', color = 'hzID', hz.depths = TRUE, depth.axis = FALSE, cex.names = 0.85))
+
+plotSPC(x, color = '.overlapFlag')
+
+plotSPC(x, color = 'hzID', show.legend = FALSE)
+plotSPC(x, color = 'claytotest')
+
+
+## TODO: fillHzGaps() adding bogus horizons
+xx <- fillHzGaps(x[1, ])
+plotSPC(xx)
+
+
+
+## test dice()
+
+# ok
+(d <- dice(x, fm = 30 ~ claytotest + .overlapFlag, SPC = FALSE))
+
+# ok: no error thrown
+(d <- dice(x, fm = 30 ~ claytotest + .overlapFlag, SPC = FALSE, strict = TRUE))
+
+# not right: extra NA slices...?
+d <- dice(x, fm = 25:45 ~ claytotest + .overlapFlag)
+plotSPC(d, color = 'claytotest', cex.names = 0.5, hz.depths = FALSE, depth.axis = TRUE)
+horizons(d[1, ])
+
+d <- dice(x, fm = 28:40 ~ claytotest + .overlapFlag)
+plotSPC(d, color = 'claytotest', cex.names = 0.5, hz.depths = FALSE, depth.axis = TRUE)
+horizons(d)
+
+
+
+
+## test slab()
+# wt.mean is correct
+(a <- slab(x, upedonid ~ claytotest, slab.structure = c(0, 50), slab.fun = mean, na.rm = TRUE))
+
+profileApply(x, function(i) {
+  i <- trunc(i, 0, 50)
+  .v <- i$claytotest
+  .w <- i$hzdepb - i$hzdept
+  idx <- which(!is.na(.v) & !is.na(.w))
+  weighted.mean(.v[idx], w = .w[idx])
+})
+