adding iteration over groups

R/simulateColor.R

@@ -10,8 +10,8 @@
     stop('package `mvtnorm` is required for multivariate simulation', call. = FALSE)
   }
 
-  ## TODO: iterate over list elements
-  .hvc <- parameters[[1]][['hvc']]
+  # extract parameters
+  .hvc <- parameters[['hvc']]
 
   # convert Munsell -> CIELAB
   .lab <- munsell2rgb(the_hue = .hvc$hue, the_value = .hvc$value, the_chroma = .hvc$chroma, returnLAB = TRUE)
@@ -20,6 +20,9 @@
   .lab <- na.omit(.lab)
 
   ## TODO: stop if nrow(.lab) < 3
+  if(nrow(.lab) < 3) {
+    return(NULL)
+  }
 
   # multivariate simulation
   # assuming approx. joint normal distribution over L, A, B coordinates
@@ -29,9 +32,12 @@
     sigma = cov(.lab),
   )
 
+
   # CIELAB -> Munsell hue, value, chroma
   m <- col2Munsell(s, space = 'CIELAB')
 
+  ## TODO: consider including only hues in reference set
+
   # flatten to standard notation
   m <- sprintf('%s %s/%s', m$hue, m$value, m$chroma)
 
@@ -111,12 +117,12 @@
   }
 
 
-  ## TODO: think about alternatives
+  ## TODO: think about alternatives:
+  #
+  #  * ? --> perform conversion without RV, then re-add just before sampling 
+  #  * z <- z[z$munsell != m$queryColor, ]
 
-  # ? --> perform conversion without RV, then re-add just before sampling 
-  # z <- z[z$munsell != m$queryColor, ]
-
-  # convert distances -> similarities
+  ## convert distances -> similarities, interpret as sampling weights
 
   # standard conversion
   # too fast of a drop off between RV and simulated values
@@ -126,16 +132,15 @@
   # simulated values too close to RV
   # s <- 1 - (z$dE00 / max(z$dE00))
 
-  ## according to ?sample there is no need to convert weights -> probabilities
 
-
-  # ## diagnostics for dE00 -> probability
+  ## diagnostics for dE00 -> probability
   # plot(s, z$dE00, type = 'n', las = 1)
   # points(s, z$dE00, col = z$color, pch = 15)
   # text(s, 0, z$munsell, cex = 0.5, srt = 90)
 
 
   # sample with replacement
+  # according to ?sample, there is no need to convert weights -> probabilities
   # using translated dE00 as prior probabilities
   res <- sample(z$munsell, replace = TRUE, size = n, prob = s)
 
@@ -147,14 +152,21 @@
 
 #' @title Simulate Soil Colors
 #' 
-#' @description Simulate plausible soil colors based on proportions by Munsell "chip", or using a seed Munsell chip and threshold specified via CIE2000 color contrast metric.
+#' @description Simulate plausible soil colors based on several possible parameterization of a "range in characteristics" (RIC). Soil color RIC can be specified by a list of parameters:
+#'  * soil color proportions, as output from [aggregateColor()] -- `method = 'proportions'`
+#'  * most likely Munsell color, CIE2000 threshold, and vector of acceptable hues -- `method = 'dE00'`
+#'  * `data.frame` of Munsell hue, value, and chroma representing observed soil colors -- `method = 'mvnorm'`
+#' 
+#' 
 #' 
 #' @author D.E. Beaudette
 #'
 #' @param method simulation method, see details
-#' @param n number of simulated colors per horizon
+#' 
+#' @param n number of simulated colors per group
+#' 
 #' @param parameters a `list`, format depends on `method`:
-#'   * `proportions`: output from [`aggregateColor`] 
+#'   * `proportions`: output from [aggregateColor()] 
 #'   * `dE00`: formatted as `list(m = '7.5YR 3/3', thresh = 5, hues = c('7.5YR'))`
 #'   * `mvnorm`: formatted as `list(hvc = x)`
 #'   
@@ -220,18 +232,26 @@ simulateColor <- function(method = c('dE00', 'proportions', 'mvnorm'), n, parame
     parameters <- list(parameters)
   }
 
+  ## TODO: basic error checking, depends on method
+
+  # select method
   res <- switch(
     method,
     'dE00' = {
+      # manual iteration over parameters
       lapply(parameters, function(i) {
         .simulateColorFromDE00(n = n, parameters = i)
       })
     },
+    # automatic iteration over output from aggregateColor()
     'proportions' = {
       .simulateColorFromProportions(n = n, parameters = parameters)
     },
+    # manual iteration over parameters
     'mvnorm' = {
-      .simulateColorFromMV(n = n, parameters = parameters)
+      lapply(parameters, function(i) {
+        .simulateColorFromMV(n = n, parameters = i)
+      })
     }
   )
 
@@ -260,6 +280,5 @@ simulateColor <- function(method = c('dE00', 'proportions', 'mvnorm'), n, parame
     return(combine(l))
   }
 
-
 }
 

misc/sandbox/simulateColor.R

@@ -3,20 +3,62 @@ library(mvtnorm)
 
 data(loafercreek, package = 'soilDB')
 
+# generalize horizon names using REGEX rules
+n <- c('Oi', 'A', 'BA','Bt1','Bt2','Bt3','Cr','R')
+p <- c('O', '^A$|Ad|Ap|AB','BA$|Bw',
+       'Bt1$|^B$','^Bt$|^Bt2$','^Bt3|^Bt4|CBt$|BCt$|2Bt|2CB$|^C$','Cr','R')
+loafercreek$genhz <- generalize.hz(loafercreek$hzname, n, p)
+
+# remove non-matching generalized horizon names
+loafercreek$genhz[loafercreek$genhz == 'not-used'] <- NA
+loafercreek$genhz <- factor(loafercreek$genhz)
+
+# all colors
 .hvc <- data.frame(
   hue = loafercreek$m_hue,
   value = loafercreek$m_value,
-  chroma = loafercreek$m_chroma,
-  stringsAsFactors = FALSE
+  chroma = loafercreek$m_chroma
 )
 
 p <- list(
   list(hvc = .hvc)
 )
 
-m <- simulateColor(method = 'mvnorm', n = 10, parameters = p)
+# result is a list
+m <- simulateColor(method = 'mvnorm', n = 100, parameters = p)
+
+colorChart(m[[1]])
+
+
+
+# by genhz
+h <- horizons(loafercreek)
+h <- split(h, h$genhz)
+
+p <- lapply(h, function(i) {
+  .res <- data.frame(
+    hue = i$m_hue,
+    value = i$m_value,
+    chroma = i$m_chroma
+  )
+
+  return(list(hvc = na.omit(.res)))
+})
+
+# some genhz have less than required (3) number of rows
+sapply(p, sapply, nrow)
 
-colorChart(m)
+# safely handle parameters without enough data
+# 25 simulations of each
+m <- simulateColor(method = 'mvnorm', n = 25, parameters = p)
+
+# invert list -> labeled rows in data.frame
+# NULL elements dropped
+m <- stack(m)
+
+# inspect results
+colorChart(m$values, m$ind)
+colorChart(m$values, m$ind, annotate = TRUE, size = FALSE, chip.cex = 2)
 
 
 
@@ -60,17 +102,6 @@ update(pp, asp = 1)
 
 
 
-data(loafercreek, package = 'soilDB')
-
-# generalize horizon names using REGEX rules
-n <- c('Oi', 'A', 'BA','Bt1','Bt2','Bt3','Cr','R')
-p <- c('O', '^A$|Ad|Ap|AB','BA$|Bw',
-       'Bt1$|^B$','^Bt$|^Bt2$','^Bt3|^Bt4|CBt$|BCt$|2Bt|2CB$|^C$','Cr','R')
-loafercreek$genhz <- generalize.hz(loafercreek$hzname, n, p)
-
-# remove non-matching generalized horizon names
-loafercreek$genhz[loafercreek$genhz == 'not-used'] <- NA
-loafercreek$genhz <- factor(loafercreek$genhz)
 
 
 cols <- data.frame(
@@ -128,7 +159,7 @@ zz <- combine(z, s)
 
 # cool
 par(mar = c(0, 0, 1, 0))
-plotSPC(zz, name.style = 'center-center', hz.depths = TRUE, plot.depth.axis = FALSE, width = 0.3)
+plotSPC(zz, name.style = 'center-center', hz.depths = TRUE, depth.axis = FALSE, width = 0.3)
 title('aggregateColor based simulation')
 
 
@@ -185,8 +216,12 @@ z <- simulateColor(method = 'dE00', n = n.sim, parameters = p, SPC = z)
 zz <- combine(z, s)
 
 # cool
-par(mar = c(0, 0, 1, 0))
-plotSPC(zz, name.style = 'center-center', hz.depths = TRUE, plot.depth.axis = FALSE, width = 0.3)
-title('dE00 based simulation')
+par(mar = c(0, 0, 0.5, 0))
+plotSPC(zz, name.style = 'center-center', hz.depths = TRUE, depth.axis = FALSE, width = 0.3, lwd = 0.5)
+title('dE00 based simulation', line = -2)
+
+par(mar = c(0, 0, 0.5, 2))
+plotSPC(zz, name.style = 'center-center', width = 0.35, lwd = 0.5, cex.names = 0.7, cex.id = 0.5, max.depth = 100)
+title('dE00 based simulation', line = -2)