add checks for illogical slab() structure

[brownag]

Encountered this error when trying to calcualate slab means on the control section. Need to add checks in slab to verify:

• PSCS top depth is above bottom depth
• PSCS top and bottom depths are not the same
• PSCS bottom depth is not greater than the maximum horizon depth

Reproducible example (until the data are fixed in NASIS):

slabbing function for control section

f.pcs.prop <- function(x, prop) {
   sv <- c(x$pscs_top, x$pscs_bottom)
   if(any(is.na(sv)))
     return(NA)
   fm <- as.formula(paste('~', prop))
   s <- slab(x, fm, slab.structure=sv, slab.fun=mean, na.rm=TRUE)
   return(s$value)
 }

Case 1

library(soilDB)
testdata=fetchKSSL(pedlabsampnum = 'UCD08JER010')
profileApply(testdata,FUN=f.pcs.prop,prop='clay')

Case 2

testdata=fetchKSSL(pedlabsampnum = 'UCD08JER010')
testdata@horizons[5,]$pscs_top=testdata@horizons[5,]$pscs_bot
profileApply(testdata,FUN=f.pcs.prop,prop='clay')

both result in:

Error in `$<-.data.frame`(`*tmp*`, "seg.label", value = integer(0)) : 
  replacement has 0 rows, data has 2

Case 3

testdata=fetchKSSL(pedlabsampnum = 'UCD08JER010')
testdata@horizons[5,]$hzn_bot=82
profileApply(testdata,FUN=f.pcs.prop,prop='clay')

results in: Error in rep(NA, times = max.d - slab.structure[2]) : invalid 'times' argument

[brownag]

Updated reprex. 2 of 3 items from OP are resolved now after prior slab upgrades.

library(aqp)
#> This is aqp 2.0.3
data("jacobs2000", package="aqp")

# same depth: properly errors
slab(jacobs2000, ~clay, slab.structure = c(50, 50))
#> Error: Invalid slab.structure

# large bottom depth: works as it should
slab(jacobs2000, ~clay, slab.structure = c(50, 10000))
#>   variable all.profiles p.q5 p.q25 p.q50 p.q75 p.q95 contributing_fraction top
#> 1     clay            1  0.4     2    20    39    47            0.01305097  50
#>   bottom
#> 1  10000

The depth logic error still needs a message/handling:

# depth logic error: needs informative message
slab(jacobs2000, ~clay, slab.structure = c(75, 50))
#> Error in rep(1, j[x]): invalid 'times' argument

Also, related to the very deep bottom depth above, with very large numbers, there is a limit in the formula parsing code used in dice(). The scipen option affects the formula that is created.

This can be fixed in the regex dice uses to get terms of the formula by ignoring the "e+" of scientific notation when finding the property variables being summarized.

# integer values: scientific notation breaks parsing depending on scipen option (>=1e5)
slab(jacobs2000, ~clay, slab.structure = c(50, 100000))
#> Error in parse(text = fm[[1]]): <text>:1:4: unexpected input
#> 1: 50:1e
#>        ^

slab(jacobs2000, ~clay, slab.structure = c(50, 1e5))
#> Error in parse(text = fm[[1]]): <text>:1:4: unexpected input
#> 1: 50:1e
#>        ^

slab(jacobs2000, ~clay, slab.structure = c(50, 100000L))
#> Error in parse(text = fm[[1]]): <text>:1:4: unexpected input
#> 1: 50:1e
#>       

Setting scipen higher works fine:

options(scipen=99)
slab(jacobs2000, ~clay, slab.structure = c(50, 100000))
#>   variable all.profiles p.q5 p.q25 p.q50 p.q75 p.q95 contributing_fraction top
#> 1     clay            1  0.4     2    20    39    47           0.001299221  50
#>   bottom
#> 1 100000

Related to https://github.com/ncss-tech/aqp/issues/106, although considering slab function arguments rather than input data depths. slab() can partially take decimal slab.structure.

I don't know whether slab should force slab.structure to be an integer. The following appears to work, but output depths are converted to integer.

# decimal values: constant slab interval, works but output depths converted to integer
slab(jacobs2000[1,], ~clay, slab.structure = c(5.5)) |> head()
#>   variable all.profiles p.q5 p.q25 p.q50 p.q75 p.q95 contributing_fraction top
#> 1     clay            1    7     7     7     7     7                     1   0
#> 2     clay            1    7     7     7     7     7                     1   5
#> 3     clay            1    7     7     7     7     7                     1  11
#> 4     clay            1    6     6     7     7     7                     1  16
#> 5     clay            1    6     6     6     6     6                     1  22
#> 6     clay            1    6     6     6     6     6                     1  27
#>   bottom
#> 1      5
#> 2     11
#> 3     16
#> 4     22
#> 5     27
#> 6     33

For the more general ideas in #106, I think changing depth_units() to something with appropriate precision to represent the decimal is the way to go.

The following still don't work:

# decimal values: single interval slab structure, index breaks
slab(jacobs2000, ~clay, slab.structure = c(5.5, 10))
#> Error in `[<-.data.frame`(`*tmp*`, , g, value = 1): replacement has 1 row, data has 0
slab(jacobs2000, ~clay, slab.structure = c(5.5, 10.5))
#> Error in `[<-.data.frame`(`*tmp*`, , g, value = 1): replacement has 1 row, data has 0

# decimal values: multiple interval slab structure with decimals returns nothing
slab(jacobs2000, ~clay, slab.structure = c(5.5, 10.5, 15.5))
#>  [1] variable              all.profiles          p.q5                 
#>  [4] p.q25                 p.q50                 p.q75                
#>  [7] p.q95                 contributing_fraction top                  
#> [10] bottom               
#> <0 rows> (or 0-length row.names)