From 1511c6228a448f9a91f3b6b049b904a591ae22e5 Mon Sep 17 00:00:00 2001
From: Beaudette <dylan.beaudette@usda.gov>
Date: Wed, 9 Oct 2024 13:19:08 -0700
Subject: [PATCH] updated experiments, and a new one related to soil color

---
 misc/sandbox/mps-smooth-soil-color.R |  81 +++++++++++++
 misc/sandbox/mps2-lambda-eval.R      | 174 ++++++++++++++++++++++-----
 2 files changed, 228 insertions(+), 27 deletions(-)
 create mode 100644 misc/sandbox/mps-smooth-soil-color.R

diff --git a/misc/sandbox/mps-smooth-soil-color.R b/misc/sandbox/mps-smooth-soil-color.R
new file mode 100644
index 000000000..ed93a4a1d
--- /dev/null
+++ b/misc/sandbox/mps-smooth-soil-color.R
@@ -0,0 +1,81 @@
+library(aqp)
+library(mpspline2)
+
+# example soil profile with some wild colors
+x <- list(
+  id = 'P1',
+  depths = c(5, 25, 33, 100, 150),
+  name = c('A', 'Bt1', 'Bt2', 'BC', 'Cr'),
+  p1 = c(5, 25, 35, 10, 8),
+  color = c('10YR 2/1', '7.5YR 3/3', '2.5Y 8/2', '2.5YR 4/6', '5G 6/3'),
+  hzd = c('clear', 'clear', 'abrupt', 'gradual', NA)
+)
+
+# init SPC
+x <- quickSPC(x)
+x$hzd <- hzDistinctnessCodeToOffset(x$hzd)
+
+# convert Munsell -> sRGB in hex notation
+x$col_source <- parseMunsell(x$color)
+
+# convert Munsell -> CIELAB
+.lab <- parseMunsell(x$color, convertColors = TRUE, returnLAB  = TRUE)
+
+# shortcut to splice-in CIELAB color coordinates
+replaceHorizons(x) <- cbind(horizons(x), .lab)
+
+# check
+plotSPC(x, color = 'L', hz.distinctness.offset = 'hzd')
+
+# hack to smooth multiple variables
+# future enhancement to spc2mpspline()
+.lambda <- 0.1
+.spcL <- spc2mpspline(x, var_name = 'L', lam = .lambda, method = 'est_1cm')
+.spcA <- spc2mpspline(x, var_name = 'A', lam = .lambda, method = 'est_1cm')
+.spcB <- spc2mpspline(x, var_name = 'B', lam = .lambda, method = 'est_1cm')
+
+m <- .spcL
+m$A_spline <- .spcA$A_spline
+m$B_spline <- .spcB$B_spline
+
+
+# check
+# ... negative numbers truncated at 0
+par(mar = c(0, 0, 3, 3))
+plotSPC(m, color = 'L_spline', name = NA, lwd = 0, divide.hz = FALSE)
+plotSPC(m, color = 'A_spline', name = NA, lwd = 0, divide.hz = FALSE)
+plotSPC(m, color = 'B_spline', name = NA, lwd = 0, divide.hz = FALSE)
+
+# back-transform to Munsell at this point
+.lab <- horizons(m)[, c('L_spline', 'A_spline', 'B_spline')]
+names(.lab) <- c('L', 'A', 'B')
+
+# interesting...
+.mun <- col2Munsell(.lab, space = 'CIELAB')
+table(.mun$hue)
+table(.mun$value)
+table(.mun$chroma)
+
+# convert smoothed CIELAB -> sRGB
+.srgb <- convertColor(horizons(m)[, c('L_spline', 'A_spline', 'B_spline')], from = 'Lab', to = 'sRGB', from.ref.white = 'D65', to.ref.white = 'D65')
+
+# sRGB -> hex notation
+m$col_spline <- rgb(.srgb, maxColorValue = 1)
+
+# ok
+plotSPC(m, color = 'col_spline', name = NA, lwd = 0, divide.hz = FALSE)
+
+# normalize names and combine SPCs
+m$soil_color <- m$col_spline
+x$soil_color <- x$col_source
+
+profile_id(m) <- 'P1-EA Spline'
+
+z <- combine(x, m)
+
+# compare side by side
+par(mar = c(0, 0, 0, 3))
+plotSPC(z, color = 'soil_color', name = NA, lwd = 0, divide.hz = FALSE, cex.names = 1)
+
+# green hues lost due to truncation of smoothed values at x>=0
+
diff --git a/misc/sandbox/mps2-lambda-eval.R b/misc/sandbox/mps2-lambda-eval.R
index 22019e072..724dcc734 100644
--- a/misc/sandbox/mps2-lambda-eval.R
+++ b/misc/sandbox/mps2-lambda-eval.R
@@ -1,26 +1,104 @@
 # install mpsline2 from CRAN
 
 library(aqp)
+library(soilDB)
 library(sharpshootR)
 library(mpspline2)
 library(reshape2)
 library(lattice)
 
-# make some example data
-ids <- LETTERS[1]
 
-set.seed(10101)
-x <- lapply(ids, random_profile, n = c(6, 7, 8), n_prop = 1, method = 'LPP', SPC = TRUE, lpp.a = 5, lpp.b = 25, lpp.d = 10, lpp.e = 5, lpp.u = 25)
-x <- combine(x)
+# x <- fetchKSSL(series = 'clarksville')
+# 
+# plotSPC(x[2, ], color = 'clay')
+# 
+# x <- x[2, ]
+# x$p1 <- x$clay
+# x$top <- x$hzn_top
+# x$bottom <- x$hzn_bot
+# profile_id(x) <- 'clarksville'
+# horizonDepths(x) <- c('top', 'bottom')
+
+
+# # # make some example data
+# ids <- LETTERS[1:3]
+# 
+# set.seed(10101)
+# x <- lapply(
+#   ids,
+#   random_profile,
+#   n = c(6, 7, 8),
+#   n_prop = 1,
+#   method = 'LPP',
+#   SPC = TRUE,
+#   lpp.a = 5,
+#   lpp.b = 30,
+#   lpp.d = 10,
+#   lpp.e = 5,
+#   lpp.u = 25
+# )
+# 
+# x <- combine(x)
+# horizons(x)$hzd <- 0
+# site(x)$group <- profile_id(x)
+
+
+
+## example data
+
+x <- list(
+  id = 'P1',
+  depths = c(5, 25, 33, 100, 150),
+  name = c('A', 'Bt1', 'Bt2', 'BC', 'Cr'),
+  p1 = c(5, 25, 35, 10, 8),
+  soil_color = c('10YR 2/2', '10YR 3/3', '10YR 4/4', '10YR 4/6', '5G 6/2'),
+  hzd = c('clear', 'clear', 'abrupt', 'gradual', NA)
+)
+
+x <- quickSPC(x)
+x$hzd <- hzDistinctnessCodeToOffset(x$hzd)
+site(x)$group <- 'A'
+
+
+## wait for this to be fixed
+# https://github.com/obrl-soil/mpspline2/issues/9
+# x$p1 <- scale(x$p1)
+
+
+
+## PAWS example
+# 
+# x <- list(
+#   id = 'P1',
+#   depths = c(5, 25, 33, 100, 150),
+#   name = c('A', 'Bt1', 'Bt2', 'BC', 'Cr'),
+#   awc_r = c(0.11, 0.15, 0.18, 0.08, 0.05),
+#   soil_color = c('10YR 2/2', '10YR 3/3', '10YR 4/4', '10YR 4/6', '5G 6/2'),
+#   hzd = c('clear', 'clear', 'abrupt', 'gradual', NA)
+# )
+# 
+# 
+# x <- quickSPC(x)
+# x$hzd <- hzDistinctnessCodeToOffset(x$hzd)
+# 
+# x$p1 <- (x$bottom - x$top) * x$awc_r
+
+
+
 
 # fake site data
 site(x)$fake_site_attr <- 1:length(x)
 
+
+.cols <- hcl.colors(50, 'spectral')
+
 # check source data
-par(mar=c(0,0,3,1))
-plotSPC(x, color='p1', col.palette = hcl.colors(10, 'viridis'))
+par(mar = c(0, 0, 3, 1))
+plotSPC(x, color = 'p1', col.palette = .cols, hz.distinctness.offset = 'hzd')
 
-s <- seq(0, 1, by = 0.25)
+# iterate over possible lambda values
+# 0.1 is default
+s <- c(0.1, 0.25, 0.3, 0.5, 1)
 m <- lapply(s, function(i) {
   .spc <- spc2mpspline(x, var_name = 'p1', lam = i, method = 'est_1cm')
   profile_id(.spc) <- sprintf("%s-0%s", profile_id(.spc), i)
@@ -32,35 +110,77 @@ z$p1 <- z$p1_spline
 
 xx <- combine(z, x)
 
+# site(xx)$id <- profile_id(xx)
+
+
 par(mar = c(0, 0, 3, 1))
-plotSPC(xx, color = 'p1', col.palette = hcl.colors(10, 'viridis'), divide.hz = FALSE, width = 0.35)
+
+plotSPC(xx, color = 'p1', col.palette = .cols, divide.hz = FALSE, width = 0.35, name = NA, lwd = 0, hz.distinctness.offset = 'hzd')
 
 
 o <- order(xx$p1_rmse)
-plotSPC(xx, color = 'p1', col.palette = hcl.colors(10, 'viridis'), divide.hz = FALSE, width = 0.35, plot.order = o)
+plotSPC(xx, color = 'p1', col.palette = .cols, divide.hz = FALSE, width = 0.35, plot.order = o, name = NA, lwd = 0, hz.distinctness.offset = 'hzd')
+
+
+.cols <- c(hcl.colors(n = 5), 'firebrick')
+tps <- tactile::tactile.theme(superpose.line = list(lwd = 2, col = .cols))
 
+site(xx)$lambda <- factor(c('original', s))
 
-tps <- tactile::tactile.theme(superpose.line = list(lwd = 2))
 
 # compare depth-functions by method, no aggregation
-xyplot(cbind(top, bottom) ~ p1, id=xx$id, groups = factor(id),
-       data=as(xx, 'data.frame'),
-       par.settings=tps,
-       auto.key=list(columns=3, lines=TRUE, points=FALSE),
-       strip=strip.custom(bg=grey(0.85)),
-       ylim=c(125, -5), as.table=TRUE, panel=panel.depth_function,
-       scales = list(alternating = 1),
-       asp = 1
+xyplot(
+  cbind(top, bottom) ~ p1 | group, 
+  id = factor(xx$id), 
+  groups = lambda,
+  data = as(xx, 'data.frame'),
+  par.settings = tps,
+  auto.key = list(columns = 1, lines = TRUE, points = FALSE, title = 'lambda', space = 'right', cex = 0.8),
+  strip = strip.custom(bg = grey(0.85)),
+  ylim = c(160, -5), 
+  ylab = 'Depth (cm)',
+  # main = '',
+  as.table = TRUE, 
+  panel = panel.depth_function,
+  scales = list(alternating = 1, y = list(tick.number = 15)),
+  asp = 2
 )
 
-xyplot(cbind(top, bottom) ~ p1 | id, id=xx$id,
-       data=as(xx, 'data.frame'),
-       par.settings=tps,
-       auto.key=list(columns=4, lines=TRUE, points=FALSE),
-       strip=strip.custom(bg=grey(0.85)),
-       ylim=c(125, -5), as.table=TRUE, panel=panel.depth_function,
-       layout = c(6, 1),
-       scales = list(alternating = 1)
+
+# note: can only be used on single-profile groups
+slabInterval <- function(i, .f = sum) {
+  .ss <- as.numeric(strsplit(i, '-', fixed = TRUE)[[1]])
+  
+  .a <- slab(xx, lambda ~ p1,  slab.structure = .ss, slab.fun = .f, na.rm = TRUE)
+  .a <- dcast(.a, lambda ~ variable, value.var = 'value')
+  names(.a)[2] <- i
+  
+  return(.a)
+}
+
+l <- lapply(
+  c('0-150', '25-50', '75-100', '0-25', '0-75'), 
+  slabInterval
 )
 
+a <- Reduce(merge, l)
+
+
+knitr::kable(a, row.names = FALSE, digits = 0, caption = 'Sum of values over depth interval.')
+
+
 
+# xyplot(cbind(top, bottom) ~ p1 | id, id=xx$id,
+#        data=as(xx, 'data.frame'),
+#        par.settings=tps,
+#        auto.key=list(columns=4, lines=TRUE, points=FALSE),
+#        strip=strip.custom(bg=grey(0.85)),
+#        ylim=c(150, -5), as.table=TRUE, panel=panel.depth_function,
+#        layout = c(6, 1),
+#        scales = list(alternating = 1)
+# )
+# 
+# 
+# 
+# 
+#