estimatePSCS: bug fix for shallow soils with restriction >36cm from s…

…oil surface, but <36cm from MINERAL soil surface

 - thanks to Tim Riebe for reporting

R/estimatePSCS.R

@@ -96,10 +96,13 @@ estimatePSCS = function(p, hzdesgn = "hzname", clay.attr = "clay",
   default_t <- rep(25, length(soildepth))
   default_b <- rep(100, length(soildepth))
 
+  odepth <- getMineralSoilSurfaceDepth(p, hzdesgn, simplify = FALSE)[[hz.depths[2]]]
+  ohzidx <- which(odepth > 0)
+
   # Key part A (soils with restriction in shallow depth)
-  lt36idx <- which(soildepth <= 36)
+  lt36idx <- which(soildepth - odepth <= 36)
   if (length(lt36idx) > 0) {
-    default_t[lt36idx] <- 0
+    default_t[lt36idx] <- 0 # O horizon correction applied below
     default_b[lt36idx] <- soildepth[lt36idx]
     shallow_flag[lt36idx] <- TRUE
   }
@@ -115,8 +118,6 @@ estimatePSCS = function(p, hzdesgn = "hzname", clay.attr = "clay",
   }
 
   # Adjust PSCS range downward if organic soil material is present at surface (i.e. mineral soil surface depth > 0)
-  odepth <- getMineralSoilSurfaceDepth(p, hzdesgn, simplify = FALSE)[[hz.depths[2]]]
-  ohzidx <- which(odepth > 0)
   if (length(ohzidx) > 0) {
 
     default_t[ohzidx] <- default_t[ohzidx] + odepth[ohzidx]

tests/testthat/test-estimatePSCS.R

@@ -6,10 +6,22 @@ data(sp1, package = 'aqp')
 depths(sp1) <- id ~ top + bottom
 site(sp1) <- ~ group
 
-p <- sp1[1]
+p <- sp1[1,]
 attr <- 'prop' # clay contents %
 
-q <- sp1[2]
+q <- sp1[2,]
+
+x <- data.frame(
+  peiid = 706300,
+  taxsubgrp = "Lithic Humicryods",
+  top = c(0, 13, 16, 18, 24, 40),
+  bottom = c(13, 16, 18, 24, 40, 65),
+  name = c("Oi", "A", "E", "Bhs", "2C", "2R"),
+  texture = c("SPM", "SIL", "SIL", "SIL", "SIL", "BR"),
+  prop = c(0, 6, 6, 6, 6, 6)
+)
+depths(x) <- peiid ~ top + bottom
+site(x) <- ~ taxsubgrp
 
 test_that("estimatePSCS()", {
 
@@ -64,6 +76,14 @@ test_that("estimatePSCS()", {
   expect_error(estimatePSCS(q2, clay.attr = 'prop', texcl.attr = "texture", hzdesgn = 'foo'))
 })
 
+test_that("estimatePSCS() thin soil profile with O horizon", {
+  expect_equal(estimatePSCS(x, clay.attr = 'prop', texcl.attr = "foo", hzdesgn = 'name'), c(13, 40))
+  expect_equal(estimatePSCS(c(q,x), clay.attr = 'prop', texcl.attr = "foo", hzdesgn = 'name'),
+               data.frame(id = c("706300", "P002"), 
+                          pscs_top = c(13, 30),
+                          pscs_bottom = c(40, 59)))
+})
+
 test_that("estimatePSCS() multiple profiles",{
   e <- estimatePSCS(sp1, clay.attr = 'prop', texcl.attr = "texture", hzdesgn = 'name')
   expect_equal(e$pscs_top, c(49, 30, 2, 32, 5, 31, 25, 27, 28))