brownag/mixed-colors-mollic.R

## mixed-colors-mollic.R
library(aqp)
library(soilDB)

## aqp_mixed_colors
#
#  The goal here is to compare dry and moist colors when numerically "mixed" over certain depth intervals
#  `aqp_mixed_colors` is vectorized over profiles -- performs color "mixing" by profile on interval [ztop, zbot]
#
#  And to compare "mixing" methods -- using a simple depth-weighted average / dominant hue approach
#                                   and the powerful, but computationally costly, method in aqp::aggregateColor
## Methods
##   1. depth-weighted average value and chroma + dominant hue in interval ("simple" aggregation)
##   2. conversion of Munsell -> RGB -> CIELAB2000 mix -> RGB -> Munsell (via aqp::aggregateColor)
#
# Returning results in a three part list reflecting the various methods and color data used.
#
aqp_mixed_colors <- function(p, ztop = 0, zbot = 18,
                         m_hue = "m_hue", m_value = "m_value", m_chroma = "m_chroma",
                         d_hue = "d_hue", d_value = "d_value", d_chroma = "d_chroma") {

  # trunc: new wrapper function around glomApply for constant top and bottom depth
  p.trc <- trunc(p, ztop, zbot)

  # calculate depth weights (horizon thickness)
  hzd <- horizonDepths(p)
  p.trc$wt <- p.trc[[hzd[2]]] - p.trc[[hzd[1]]]

  # iterate over profiles, calculating depth-weighted average value/chroma + dominant hues
  mixed <- profileApply(p.trc,
                        frameify = TRUE, # frameify = TRUE means return data.frame result for _each profile_
                        function(p.sub) {

    dwt <- aggregate(p.sub[["wt"]], by=list(p.sub[[d_hue]]), sum)
    mwt <- aggregate(p.sub[["wt"]], by=list(p.sub[[m_hue]]), sum)

    if (!nrow(dwt))
      dwt <- data.frame(Group=NA, x = 1)

    if (!nrow(mwt))
      mwt <- data.frame(Group=NA, x = 1)

    # construct result
    res <- data.frame(
        id = profile_id(p.sub), # profile ID
        n_hz = nrow(p.sub),     # number of horizons
        dominant_d_hue = dwt[which(dwt$x == max(dwt$x, na.rm = TRUE))[1], 1],
        dominant_m_hue = mwt[which(mwt$x == max(mwt$x, na.rm = TRUE))[1], 1],
        wt_d_value  = weighted.mean(p.sub[[d_value]],  p.sub[["wt"]]),
        wt_m_value  = weighted.mean(p.sub[[m_value]],  p.sub[["wt"]]),
        wt_d_chroma = weighted.mean(p.sub[[d_chroma]], p.sub[["wt"]]),
        wt_m_chroma = weighted.mean(p.sub[[m_chroma]], p.sub[["wt"]]))

    # put idname into ID name slot in result
    names(res)[1] <- idname(p.sub)

    return(res)
  })

  ### calculate colors mixed in LAB space

  ## moist
  p.trc$m_lab <- aqp::munsell2rgb(p.trc[[m_hue]], p.trc[[m_value]], p.trc[[m_chroma]])
  m_res <- aqp::aggregateColor(p.trc, groups = idname(p.trc), k = 1, col = "m_lab")$aggregate.data
  m_res_missing <- which(!profile_id(p.trc) %in% m_res[[idname(p.trc)]])

  # deal with missing values
  if(length(m_res_missing)) {
    m_res.tmp <- m_res[NA,][1:length(m_res_missing),]
    m_res.tmp[[idname(p.trc)]] <- profile_id(p.trc)[m_res_missing]
    m_res <- rbind(m_res, m_res.tmp)
    m_res <- m_res[match(profile_id(p.trc), m_res[[idname(p.trc)]]),] # reapply original order
  }

  ## dry
  p.trc$d_lab <- aqp::munsell2rgb(p.trc[[d_hue]], p.trc[[d_value]], p.trc[[d_chroma]])
  d_res <- aqp::aggregateColor(p.trc, groups = idname(p.trc), k = 1, col = "d_lab")$aggregate.data
  d_res_missing <- which(!profile_id(p.trc) %in% d_res[[idname(p.trc)]])

  # deal with missing values
  if(length(d_res_missing)) {
    d_res.tmp <- m_res[NA,][1:length(d_res_missing),]
    d_res.tmp[[idname(p.trc)]] <- profile_id(p.trc)[d_res_missing]
    d_res <- rbind(d_res, d_res.tmp)
    d_res <- d_res[match(profile_id(p.trc), d_res[[idname(p.trc)]]),] # reapply original order
  }

  # construct final result
  res <- list()
  res$depth_weighted <- mixed
  res$m_aggregate <- m_res
  res$d_aggregate <- d_res

  return(res)
}

# a fake profile
dat0 <- data.frame(id       = 1,
                   top      = 0,
                   bottom   = 50,
                   d_hue    = "10YR",
                   m_hue    = "10YR",
                   d_value  = 5,
                   m_value  = 3,
                   d_chroma = 3,
                   m_chroma = 3)

# a real profile
dat1 <- data.frame(id       = 2,
                   top      = c(0, 7, 27, 43),
                   bottom   = c(7, 27, 43, 59),
                   d_hue    = c("10YR","10YR","10YR","10YR"),
                   m_hue    = c("10YR","10YR","10YR","10YR"),
                   d_value  = c(7, 4,     4, 5),
                   m_value  = c(6, 2.5, 2.5, 4),
                   d_chroma = c(2, 3,     3, 3),
                   m_chroma = c(2, 2,     2, 3))

# same color data, shuffled depths a bit so more horizons in surface 0-18cm
dat2 <- data.frame(id       = 3,
                   top      = c(0, 5, 12, 16),
                   bottom   = c(5, 12, 16, 59),
                   d_hue    = c("10YR","10YR","10YR","10YR"),
                   m_hue    = c("10YR","10YR","10YR","10YR"),
                   d_value  = c(7, 4,     4, 5),
                   m_value  = c(6, 2.5, 2.5, 4),
                   d_chroma = c(2, 3,     3, 3),
                   m_chroma = c(2, 2,     2, 3))

dat <- rbind(dat0, dat1, dat2)
depths(dat) <- id ~ top + bottom

dat_test <- aqp_mixed_colors(dat)

data("loafercreek", package = "soilDB")
res <- aqp_mixed_colors(loafercreek)

#f <- fetchNASIS()
#res <- aqp_mixed_colors(f)

#' Construct a Munsell Hue Value/Chroma Code
#'
#' @param the_hue Character vector of hue values (in the Munsell notation e.g. "10YR")
#' @param the_value Numeric vector of value values
#' @param the_chroma Numeric vector of chroma values
#' @param digits Number of digits to round value and chroma off to (default: 0; integer values)
#'
#' @return A character vector of Munsell color codes in the format \code{HUE VALUE/CHROMA}
#' @export hvc_to_munsell
#'
#' @examples
#'
#' hvc_to_munsell("10YR", 2, 2)
#' # [1] "10YR 2/2"
#'
hvc_to_munsell <- function(the_hue, the_value, the_chroma, digits = 0, allow57chroma = FALSE) {
  chroma <- the_chroma
  if(!allow57chroma) {
    # chromas of 5 and 7 are not chips available in soil color book
    # so divide input chromas in two, round,
    # multipy by 2 to turn 5 -> [4 or 6] or 7 -> [6 or 8]
    idx <- which(round(chroma, digits) %in% c(5,7))
    if(length(idx))
      chroma[idx] <- 2 * round(chroma[idx] / 2, digits)
  }
  res <- sprintf("%s %s/%s", the_hue, round(the_value, digits), round(chroma, digits))
  # allow for null chroma, convert other missing values to NA
  res[is.na(the_hue) | is.na(the_value)] <- NA

  return(res)
}

# this will round to nearest integer as needed
d_depthweight <- hvc_to_munsell(res$depth_weighted$dominant_d_hue,
                                res$depth_weighted$wt_d_value,
                                res$depth_weighted$wt_d_chroma)

m_depthweight <- hvc_to_munsell(res$depth_weighted$dominant_m_hue,
                                res$depth_weighted$wt_m_value,
                                res$depth_weighted$wt_m_chroma)

good.m.idx <- which(!grepl("NA", m_depthweight) &
                      !is.na(m_depthweight))

good.d.idx <- which(!grepl("NA", d_depthweight) &
                      !is.na(d_depthweight))


# compare depth-weighted dry versus moist -- just as a general sanity check
idx <- good.m.idx[good.m.idx %in% good.d.idx][5:10]

# idx <- c(7, 23, 70, 78, 96) # nice test IDs for loafercreek that demonstrate some things

test <- data.frame(id = res$depth_weighted$peiid[idx],
                   dry = d_depthweight[idx],
                   moist = m_depthweight[idx])

aqp::colorContrastPlot(test$dry, test$moist, # the moist colors (m2) are darker, as expected
                       labels = c("Dry", "Moist"))

# this is a known "easter egg" in the loafercreek dataset
# representing a fairly easy-to-commit data entry error (v.s. 7.5YR)
# note that the test indices includes one bogus "7.5R" hue in dry and has prominent contrast/high

# get the aggregated munsell chips from the 0-18cm interval
d_aqp <- hvc_to_munsell(res$d_aggregate$munsell.hue,
                        res$d_aggregate$munsell.value,
                        res$d_aggregate$munsell.chroma)

m_aqp <- hvc_to_munsell(res$m_aggregate$munsell.hue,
                        res$m_aggregate$munsell.value,
                        res$m_aggregate$munsell.chroma)

# compare depth-weighted dry versus aqp-mixed dry
moist_data <- data.frame(id = res$depth_weighted$peiid,
                         dwt = m_depthweight,
                         aqp = m_aqp,
                         ord = res$m_aggregate$munsell.sigma)

dry_data <- data.frame(id = res$depth_weighted$peiid,
                       dwt = d_depthweight,
                       aqp = d_aqp,
                       ord = res$d_aggregate$munsell.sigma)

# remove NA values (mostly from weighted.average missing values)
moist_data <- moist_data[good.m.idx,]
dry_data <- dry_data[good.d.idx,]

# order the data based on munsell.sigma
moist_data <- moist_data[order(moist_data$ord),]
dry_data <- dry_data[order(dry_data$ord),]

idx2 <- seq(1,nrow(moist_data),8)
aqp::colorContrastPlot(moist_data$dwt[idx2], moist_data$aqp[idx2],
                       printMetrics = F, col.cex = 0.0,
                       labels = c("nearest chip\ndepth-weighted mean V + C\nw/ dominant H",
                                  "aggregateColor"))
text(0.3,0.42, "Comparison of 0-18cm Numerical Mixing Methods for the Mollic/Umbric Epipedon")

idx3 <- seq(1,nrow(dry_data),8)
aqp::colorContrastPlot(dry_data$dwt[idx3], dry_data$aqp[idx3], printMetrics = F, col.cex = 0.0,
                       labels = c("nearest chip\ndepth-weighted mean V + C\nw/ dominant H",
                                  "aggregateColor"))
text(0.3,0.42, "Comparison of 0-18cm Numerical Mixing Methods for the Mollic/Umbric Epipedon")
	library(aqp)
	library(soilDB)

	## aqp_mixed_colors
	#
	# The goal here is to compare dry and moist colors when numerically "mixed" over certain depth intervals
	# `aqp_mixed_colors` is vectorized over profiles -- performs color "mixing" by profile on interval [ztop, zbot]
	#
	# And to compare "mixing" methods -- using a simple depth-weighted average / dominant hue approach
	# and the powerful, but computationally costly, method in aqp::aggregateColor
	## Methods
	## 1. depth-weighted average value and chroma + dominant hue in interval ("simple" aggregation)
	## 2. conversion of Munsell -> RGB -> CIELAB2000 mix -> RGB -> Munsell (via aqp::aggregateColor)
	#
	# Returning results in a three part list reflecting the various methods and color data used.
	#
	aqp_mixed_colors <- function(p, ztop = 0, zbot = 18,
	m_hue = "m_hue", m_value = "m_value", m_chroma = "m_chroma",
	d_hue = "d_hue", d_value = "d_value", d_chroma = "d_chroma") {

	# trunc: new wrapper function around glomApply for constant top and bottom depth
	p.trc <- trunc(p, ztop, zbot)

	# calculate depth weights (horizon thickness)
	hzd <- horizonDepths(p)
	p.trc$wt <- p.trc[[hzd[2]]] - p.trc[[hzd[1]]]

	# iterate over profiles, calculating depth-weighted average value/chroma + dominant hues
	mixed <- profileApply(p.trc,
	frameify = TRUE, # frameify = TRUE means return data.frame result for _each profile_
	function(p.sub) {

	dwt <- aggregate(p.sub[["wt"]], by=list(p.sub[[d_hue]]), sum)
	mwt <- aggregate(p.sub[["wt"]], by=list(p.sub[[m_hue]]), sum)

	if (!nrow(dwt))
	dwt <- data.frame(Group=NA, x = 1)

	if (!nrow(mwt))
	mwt <- data.frame(Group=NA, x = 1)

	# construct result
	res <- data.frame(
	id = profile_id(p.sub), # profile ID
	n_hz = nrow(p.sub), # number of horizons
	dominant_d_hue = dwt[which(dwt$x == max(dwt$x, na.rm = TRUE))[1], 1],
	dominant_m_hue = mwt[which(mwt$x == max(mwt$x, na.rm = TRUE))[1], 1],
	wt_d_value = weighted.mean(p.sub[[d_value]], p.sub[["wt"]]),
	wt_m_value = weighted.mean(p.sub[[m_value]], p.sub[["wt"]]),
	wt_d_chroma = weighted.mean(p.sub[[d_chroma]], p.sub[["wt"]]),
	wt_m_chroma = weighted.mean(p.sub[[m_chroma]], p.sub[["wt"]]))

	# put idname into ID name slot in result
	names(res)[1] <- idname(p.sub)

	return(res)
	})

	### calculate colors mixed in LAB space

	## moist
	p.trc$m_lab <- aqp::munsell2rgb(p.trc[[m_hue]], p.trc[[m_value]], p.trc[[m_chroma]])
	m_res <- aqp::aggregateColor(p.trc, groups = idname(p.trc), k = 1, col = "m_lab")$aggregate.data
	m_res_missing <- which(!profile_id(p.trc) %in% m_res[[idname(p.trc)]])

	# deal with missing values
	if(length(m_res_missing)) {
	m_res.tmp <- m_res[NA,][1:length(m_res_missing),]
	m_res.tmp[[idname(p.trc)]] <- profile_id(p.trc)[m_res_missing]
	m_res <- rbind(m_res, m_res.tmp)
	m_res <- m_res[match(profile_id(p.trc), m_res[[idname(p.trc)]]),] # reapply original order
	}

	## dry
	p.trc$d_lab <- aqp::munsell2rgb(p.trc[[d_hue]], p.trc[[d_value]], p.trc[[d_chroma]])
	d_res <- aqp::aggregateColor(p.trc, groups = idname(p.trc), k = 1, col = "d_lab")$aggregate.data
	d_res_missing <- which(!profile_id(p.trc) %in% d_res[[idname(p.trc)]])

	# deal with missing values
	if(length(d_res_missing)) {
	d_res.tmp <- m_res[NA,][1:length(d_res_missing),]
	d_res.tmp[[idname(p.trc)]] <- profile_id(p.trc)[d_res_missing]
	d_res <- rbind(d_res, d_res.tmp)
	d_res <- d_res[match(profile_id(p.trc), d_res[[idname(p.trc)]]),] # reapply original order
	}

	# construct final result
	res <- list()
	res$depth_weighted <- mixed
	res$m_aggregate <- m_res
	res$d_aggregate <- d_res

	return(res)
	}

	# a fake profile
	dat0 <- data.frame(id = 1,
	top = 0,
	bottom = 50,
	d_hue = "10YR",
	m_hue = "10YR",
	d_value = 5,
	m_value = 3,
	d_chroma = 3,
	m_chroma = 3)

	# a real profile
	dat1 <- data.frame(id = 2,
	top = c(0, 7, 27, 43),
	bottom = c(7, 27, 43, 59),
	d_hue = c("10YR","10YR","10YR","10YR"),
	m_hue = c("10YR","10YR","10YR","10YR"),
	d_value = c(7, 4, 4, 5),
	m_value = c(6, 2.5, 2.5, 4),
	d_chroma = c(2, 3, 3, 3),
	m_chroma = c(2, 2, 2, 3))

	# same color data, shuffled depths a bit so more horizons in surface 0-18cm
	dat2 <- data.frame(id = 3,
	top = c(0, 5, 12, 16),
	bottom = c(5, 12, 16, 59),
	d_hue = c("10YR","10YR","10YR","10YR"),
	m_hue = c("10YR","10YR","10YR","10YR"),
	d_value = c(7, 4, 4, 5),
	m_value = c(6, 2.5, 2.5, 4),
	d_chroma = c(2, 3, 3, 3),
	m_chroma = c(2, 2, 2, 3))

	dat <- rbind(dat0, dat1, dat2)
	depths(dat) <- id ~ top + bottom

	dat_test <- aqp_mixed_colors(dat)

	data("loafercreek", package = "soilDB")
	res <- aqp_mixed_colors(loafercreek)

	#f <- fetchNASIS()
	#res <- aqp_mixed_colors(f)

	#' Construct a Munsell Hue Value/Chroma Code
	#'
	#' @param the_hue Character vector of hue values (in the Munsell notation e.g. "10YR")
	#' @param the_value Numeric vector of value values
	#' @param the_chroma Numeric vector of chroma values
	#' @param digits Number of digits to round value and chroma off to (default: 0; integer values)
	#'
	#' @return A character vector of Munsell color codes in the format \code{HUE VALUE/CHROMA}
	#' @export hvc_to_munsell
	#'
	#' @examples
	#'
	#' hvc_to_munsell("10YR", 2, 2)
	#' # [1] "10YR 2/2"
	#'
	hvc_to_munsell <- function(the_hue, the_value, the_chroma, digits = 0, allow57chroma = FALSE) {
	chroma <- the_chroma
	if(!allow57chroma) {
	# chromas of 5 and 7 are not chips available in soil color book
	# so divide input chromas in two, round,
	# multipy by 2 to turn 5 -> [4 or 6] or 7 -> [6 or 8]
	idx <- which(round(chroma, digits) %in% c(5,7))
	if(length(idx))
	chroma[idx] <- 2 * round(chroma[idx] / 2, digits)
	}
	res <- sprintf("%s %s/%s", the_hue, round(the_value, digits), round(chroma, digits))
	# allow for null chroma, convert other missing values to NA
	res[is.na(the_hue) \| is.na(the_value)] <- NA

	return(res)
	}

	# this will round to nearest integer as needed
	d_depthweight <- hvc_to_munsell(res$depth_weighted$dominant_d_hue,
	res$depth_weighted$wt_d_value,
	res$depth_weighted$wt_d_chroma)

	m_depthweight <- hvc_to_munsell(res$depth_weighted$dominant_m_hue,
	res$depth_weighted$wt_m_value,
	res$depth_weighted$wt_m_chroma)

	good.m.idx <- which(!grepl("NA", m_depthweight) &
	!is.na(m_depthweight))

	good.d.idx <- which(!grepl("NA", d_depthweight) &
	!is.na(d_depthweight))


	# compare depth-weighted dry versus moist -- just as a general sanity check
	idx <- good.m.idx[good.m.idx %in% good.d.idx][5:10]

	# idx <- c(7, 23, 70, 78, 96) # nice test IDs for loafercreek that demonstrate some things

	test <- data.frame(id = res$depth_weighted$peiid[idx],
	dry = d_depthweight[idx],
	moist = m_depthweight[idx])

	aqp::colorContrastPlot(test$dry, test$moist, # the moist colors (m2) are darker, as expected
	labels = c("Dry", "Moist"))

	# this is a known "easter egg" in the loafercreek dataset
	# representing a fairly easy-to-commit data entry error (v.s. 7.5YR)
	# note that the test indices includes one bogus "7.5R" hue in dry and has prominent contrast/high

	# get the aggregated munsell chips from the 0-18cm interval
	d_aqp <- hvc_to_munsell(res$d_aggregate$munsell.hue,
	res$d_aggregate$munsell.value,
	res$d_aggregate$munsell.chroma)

	m_aqp <- hvc_to_munsell(res$m_aggregate$munsell.hue,
	res$m_aggregate$munsell.value,
	res$m_aggregate$munsell.chroma)

	# compare depth-weighted dry versus aqp-mixed dry
	moist_data <- data.frame(id = res$depth_weighted$peiid,
	dwt = m_depthweight,
	aqp = m_aqp,
	ord = res$m_aggregate$munsell.sigma)

	dry_data <- data.frame(id = res$depth_weighted$peiid,
	dwt = d_depthweight,
	aqp = d_aqp,
	ord = res$d_aggregate$munsell.sigma)

	# remove NA values (mostly from weighted.average missing values)
	moist_data <- moist_data[good.m.idx,]
	dry_data <- dry_data[good.d.idx,]

	# order the data based on munsell.sigma
	moist_data <- moist_data[order(moist_data$ord),]
	dry_data <- dry_data[order(dry_data$ord),]

	idx2 <- seq(1,nrow(moist_data),8)
	aqp::colorContrastPlot(moist_data$dwt[idx2], moist_data$aqp[idx2],
	printMetrics = F, col.cex = 0.0,
	labels = c("nearest chip\ndepth-weighted mean V + C\nw/ dominant H",
	"aggregateColor"))
	text(0.3,0.42, "Comparison of 0-18cm Numerical Mixing Methods for the Mollic/Umbric Epipedon")

	idx3 <- seq(1,nrow(dry_data),8)
	aqp::colorContrastPlot(dry_data$dwt[idx3], dry_data$aqp[idx3], printMetrics = F, col.cex = 0.0,
	labels = c("nearest chip\ndepth-weighted mean V + C\nw/ dominant H",
	"aggregateColor"))
	text(0.3,0.42, "Comparison of 0-18cm Numerical Mixing Methods for the Mollic/Umbric Epipedon")