tests/testthat/test-slab.R

In ncss-tech/aqp: Algorithms for Quantitative Pedology

context("slab method for SoilProfileCollection objects")

test_that("basic slab functionality", {
  
  data(sp1, package = 'aqp')
  depths(sp1) <- id ~ top + bottom
 
  # aggregate entire collection
  # 1-unit slabs
  a.1 <- slab(sp1, fm = ~ prop, strict=TRUE)
  # 5-unit slabs
  a.2 <- slab(sp1, fm = ~ prop, strict=TRUE, slab.structure=5)
  # custom slabs
  a.3 <- slab(sp1, fm = ~ prop, strict=TRUE, slab.structure=c(0,5,10,25,50))
  
  # did it work?
  expect_true(inherits(a.1, 'data.frame'))
  expect_true(inherits(a.2, 'data.frame'))
  expect_true(inherits(a.3, 'data.frame'))
  
  # number of results, these are in long format
  expect_equal(nrow(a.1), 240)
  expect_equal(nrow(a.2), 48)
  expect_equal(nrow(a.3), 4)
  
  # required column names
  nm <- names(a.1)
  expect_true(any(grepl('variable', nm)))
  expect_true(any(grepl('all.profiles', nm)))
  expect_true(any(grepl('top', nm)))
  expect_true(any(grepl('bottom', nm)))
  expect_true(any(grepl('contributing_fraction', nm)))
  expect_true(any(grepl('p.q', nm)))
})

test_that("slab structure for a single slab", {
  data(sp4, package = 'aqp')
  depths(sp4) <- id ~ top + bottom
  
  sp4$group <- c(rep('A',5), rep('B',5))
  a.1 <- slab(
    sp4,
    fm = ~ sand + silt + clay,
    slab.structure = c(0, 10),
    slab.fun = mean,
    na.rm = TRUE
  )
  expect_equal(nrow(a.1), 3)
  
  # again, this time within groups defined by a site-level attribute:
  a.1 <- slab(
    sp4,
    fm = group ~ sand + silt + clay,
    slab.structure = c(0, 10),
    slab.fun = mean,
    na.rm = TRUE
  )
  expect_equal(nrow(a.1), 6)
})


test_that("extended slab functionality: weighted aggregation", {
  
  data(sp1, package = 'aqp')
  depths(sp1) <- id ~ top + bottom
  
  sp1$weights <- rep(1:2, length(sp1))[1:length(sp1)]
  sp1$wtgrp <- rep(1, length(sp1))
  
  # component weighted averages
  # mukey=461268; Doemill-Jokerst, 3 to 8 percent slopes (615)
  x <- data.frame(cokey = c(21469960L, 21469960L, 21469960L, 21469960L, 
                            21469960L, 21469961L, 21469961L, 21469961L), 
                  comppct_r = c(50L,  50L, 50L, 50L, 50L, 40L, 40L, 40L),
                  hzdept_r = c(0L, 3L, 13L, 23L, 36L, 0L, 3L, 10L), 
                  hzdepb_r = c(3L, 13L, 23L, 36L, 61L, 3L, 10L, 35L), 
                  ph1to1h2o_r = c(6.1, 6.8, 6.7, 6.7, NA, 6.4, 6.5, NA))
  depths(x) <- cokey ~ hzdept_r + hzdepb_r
  site(x) <- ~ comppct_r
  
  # custom function, with actual comppct_r
  a.0 <- slab(x, ~ ph1to1h2o_r,
      slab.structure = c(0, 10),
      weights = "comppct_r",
      slab.fun = weighted.mean,
      na.rm = TRUE
    )
  expect_equal(a.0$value, 6.54, tolerance = 0.005) 
  
  skip_if_not_installed("Hmisc")
  
  # we expect quantile estimates to vary (given weighted v.s. unweighted)
  a.0 <- slab(sp1, fm = ~ prop, weights = "weights")
  a.1 <- slab(sp1, fm = ~ prop, strict = TRUE, weights = "weights")
  a.2 <- slab(sp1, fm = wtgrp ~ prop, strict = TRUE, weights = "weights")
  a.3 <- slab(sp1, fm = wtgrp ~ prop, strict = TRUE)
  
  # expect consistent structure for weighted/unweighted: same column names except for group
  ungroupcols <- colnames(a.3)
  ungroupcols[2] <- "all.profiles"
  expect_equal(colnames(a.0), ungroupcols)
  expect_equal(colnames(a.1), ungroupcols)
  expect_equal(colnames(a.2), colnames(a.3))
  
  # contributing fractions should be identical
  expect_true(all(a.1$contributing_fraction == a.2$contributing_fraction))
  expect_true(all(a.2$contributing_fraction == a.3$contributing_fraction))
  
  # should match this within the tolerance:
  #   soilDB::get_SDA_property(property = 'ph1to1h2o_r', method = "weighted average",
  #                            mukeys = 461268, miscellaneous_areas = FALSE, include_minors = TRUE,
  #                            top_depth = 0, bottom_depth = 10)
                                  
  # comppct_r adjusted to get higher result more like doemill
  x$comppct_r <- c(90, 10)
  a.1 <- slab(x, ~ ph1to1h2o_r, slab.structure = c(0, 10), weights = "comppct_r")
  expect_equal(a.1$p.q50, 6.8)
  
  # comppct_r adjusted to get lower result more like jokerst
  x$comppct_r <- c(10, 90)
  a.2 <- slab(x, ~ ph1to1h2o_r, slab.structure = c(0, 10), weights = "comppct_r")
  expect_equal(a.2$p.q50, 6.5)
  
  # comppct_r adjusted to get lower result more like jokerst
  x$comppct_r <- c(NA, 90)
  na.1 <- slab(x, ~ ph1to1h2o_r, slab.structure = c(0, 10), weights = "comppct_r", na.rm = TRUE)
  expect_equal(na.1$p.q50, 6.5)
  
  na.2 <- slab(x, ~ ph1to1h2o_r, slab.structure = c(0, 10), weights = "comppct_r", slab.fun = weighted.mean, na.rm = TRUE)
  expect_equal(na.2$value, NA_real_)
  
})

test_that("slab calculations: mean, single profile", {
  
  data(sp1, package = 'aqp')
  depths(sp1) <- id ~ top + bottom
  
  # aggregate single profile
  # custom slabs
  a <- slab(
      sp1[1,],
      fm = ~ prop,
      strict = TRUE,
      slab.structure = c(0, 5, 10, 25, 100),
      slab.fun = mean,
      na.rm = TRUE
    )
  
  # using mean and similar functions will cause slab to store single result / slab into 'value' column
  expect_true(any(grepl('value', names(a))))
  
  # calculations done by hand (DEB)
  # weighted mean, single profile
  # 0-5: 9.400
  expect_equal(a$value[1], 9.4, tolerance = 0.0001)
  # 5-10: 7.000
  expect_equal(a$value[2], 7, tolerance = 0.0001)
  # 10-25: 8.466
  expect_equal(a$value[3], 8.4666, tolerance = 0.0001)
  # 25-100: 15.625
  expect_equal(a$value[4], 15.625, tolerance = 0.0001)
  
})


test_that("slab calculations: mean, several profiles", {
  
  data(sp1, package = 'aqp')
  depths(sp1) <- id ~ top + bottom
  
  # aggregate single profile
  # custom slabs
  a <- slab(sp1, fm = ~ prop, strict = TRUE, slab.structure = c(0, 5, 10, 25, 100), slab.fun = mean, na.rm = TRUE)
  
  # weighted mean calculations done by 1-unit slices
  # note slice formula notation works from horizon "tops"
  s.1 <- dice(sp1, 0:4 ~ prop, SPC = FALSE)$prop
  s.2 <- dice(sp1, 5:9 ~ prop, SPC = FALSE)$prop
  s.3 <- dice(sp1, 10:24 ~ prop, SPC = FALSE)$prop
  s.4 <- dice(sp1, 25:99 ~ prop, SPC = FALSE)$prop
  
  # 0-5
  expect_equal(a$value[1], mean(s.1, na.rm=TRUE), tolerance=0.0001)
  # 5-10
  expect_equal(a$value[2], mean(s.2, na.rm=TRUE), tolerance=0.0001)
  # 10-25
  expect_equal(a$value[3], mean(s.3, na.rm=TRUE), tolerance=0.0001)
  # 25-100
  expect_equal(a$value[4], mean(s.4, na.rm=TRUE), tolerance=0.0001)
  
})


test_that("edge case: slab.structure[2] > max(x)", {
  
  data(sp1, package = 'aqp')
  depths(sp1) <- id ~ top + bottom
  
  # single profile, to max(x)
  a.max <- slab(sp1[1, ], fm = ~ prop, strict=TRUE, slab.structure=c(0,max(sp1[1, ])), slab.fun = mean, na.rm=TRUE)
  
  # single profile
  # custom interval, exceeding max(sp1)
  a <- slab(sp1[1, ], fm = ~ prop, strict=TRUE, slab.structure=c(0,300), slab.fun = mean, na.rm=TRUE)
  
  # weighted mean should be the same
  expect_true(a.max$value == a$value)
  
  # contributing fractions will be different
  # smaller value when slab.structure[2] > max.d
  expect_true(a$contributing_fraction < a.max$contributing_fraction)
  
  # compare with hand-calculated value
  expect_equivalent(a$value, 13.58427, tolerance=0.001)
})

test_that("edge case: slab.structure[1] > 0 (w/ custom slab function)", {
  data(sp3, package = "aqp")
  depths(sp3) <- id ~ top + bottom
  
  # custom 'slab' function, returning mean +/- 1SD
  mean.and.sd <- function(values) {
    m <- mean(values, na.rm = TRUE)
    s <- sd(values, na.rm = TRUE)
    upper <- m + s
    lower <- m - s
    res <- c(mean = m,
             lower = lower,
             upper = upper)
    return(res)
  }
  
  ## this time, compute the weighted mean of selected properties, by profile ID
  a <- slab(sp3,
            fm = id ~ L + A + B,
            slab.structure = c(40, 60), 
            slab.fun = mean.and.sd
  )
  
  # convert long -> wide
  res <- data.table::dcast(
    data.table::as.data.table(a),
    formula = id + top + bottom ~ variable,
    value.var = 'mean'
  )
  
  expect_equal(nrow(res), length(sp3))
  expect_equal(ncol(res), 6L)
})

test_that("overlapping horizons", {
  data(sp4, package = 'aqp')
  depths(sp4) <- id ~ top + bottom
  
  # overlapping horizons results in increased number of slices (according to overlapping thickness)
  x1 <- slab(sp4, ~ K + Mg + Ca + CEC_7 + ex_Ca_to_Mg)
  
  # create overlap
  sp4@horizons[2,]$bottom <- sp4@horizons[2,]$bottom + 12
  
  # strict=TRUE byhz=FALSE removes overlap via whole profile
  expect_message({x2 <- slab(sp4, ~ K + Mg + Ca + CEC_7 + ex_Ca_to_Mg, strict = TRUE, byhz = FALSE)})
  expect_equal(nrow(x2), 245)
  
  # strict=TRUE / strict=FALSE (with default byhz) -- we get same number of rows with overlap as without
  expect_silent(x3 <- slab(sp4, ~ K + Mg + Ca + CEC_7 + ex_Ca_to_Mg, strict = FALSE))
  expect_equal(nrow(x3), 245)
  
  # try with a larger constant slab structure
  expect_silent(x3 <- slab(sp4, ~ K + Mg + Ca + CEC_7 + ex_Ca_to_Mg, slab.structure = 10, strict = FALSE))
  expect_equal(nrow(x3), 25)
})