1 Introduction

1.1 Perceptual Color Distance

1.2 Soil Color Contrast Class

1.3 Delta-E

2 Applications

2.1 Setup

library(aqp)
library(soilDB)
library(farver)
library(cluster)
library(ape)
library(colorspace)
library(latticeExtra)
library(knitr)
library(MASS)
library(vegan)
library(viridisLite)
library(tactile)

2.2 Core Functions

Soil Survey Technical Note 2

2.2.1 Color Contrast Metrics

Compute several indices of color contrast from 2 vectors of Munsell colors, m1 and m2. Comparisons fully vectorized, e.g. comparisons performed element-wise \(m1_i ~\leftrightarrow~ m2_i\) \(...\) \(m1_n ~\leftrightarrow~ m2_n\).

\(dH\): change in Munsell hue, see Hue Position section below
\(dV\): change in Munsell value (absolute value of difference)
\(dC\): change in Munsell chroma (absolute value of difference)
\(\Delta{E_{00}}\): CIE2000 Delta-E #1, #2
\(cc\): soil color contrast class (Soil Survey Technical Note 2)

Results as a data.frame.

# examples
m1 <- c('10YR 6/3', '7.5YR 3/3', '10YR 2/2', '7.5YR 3/4')
m2 <- c('5YR 3/4', '7.5YR 4/4', '2.5YR 2/2', '7.5YR 6/3')

# result is a data.frame
d <- colorContrast(m1, m2)

# check
kable(d, row.names=FALSE)

m1	m2	dH	dV	dC	dE00	cc
10YR 6/3	5YR 3/4	2	3	1	31.313572	Prominent
7.5YR 3/3	7.5YR 4/4	0	1	1	9.647979	Faint
10YR 2/2	2.5YR 2/2	3	0	0	6.814101	Faint
7.5YR 3/4	7.5YR 6/3	0	3	1	30.265956	Distinct

Results summarized in a figure.

# more examples
m1 <- c('10YR 6/3', '7.5YR 3/3', '10YR 2/2', '7.5YR 3/4', '2.5Y 6/8', '5B 4/6')
m2 <- c('5YR 3/4', '7.5YR 4/4', '2.5YR 2/2', '7.5YR 6/3', '10YR 2/1', '5GY 3/4')

# graphical comparison
colorContrastPlot(m1, m2)

Simulated redoximorphic feature colors, contrast classes and \(\Delta{E_{00}}\).

m1 <- paste0('7.5YR 4/', 2:8)
m2 <- rep('10YR 5/2', times=length(m1))
colorContrastPlot(m1, m2, labels = c('F3M', 'MAT'), d.cex = 0.8, col.cex = 0.8)

m1 <- paste0('5Y 4/', 5:1)
m2 <- rep('10YR 3/5', times=length(m1))
colorContrastPlot(m1, m2, labels = c('RMX/FED', 'MAT'), d.cex = 0.8, col.cex = 0.8)

m1 <- c('2.5Y 3/2', '5Y 3/2', '2.5GY 3/2', '5GY 3/2', '7.5BG 3/2', '2.5B 3/2')
m2 <- rep('10YR 3/5', times=length(m1))
colorContrastPlot(m1, m2, labels = c('RMX', 'MAT'), d.cex = 0.8, col.cex = 0.8)

2.2.2 Color Contrast Class

Test rules and cases outlined in Soil Survey Technical Note 2.

10YR 6/3 vs 5YR 3/4

contrastClass(v1=6, c1=3, v2=3, c2=4, dH=2, dV=3, dC=1, verbose = TRUE)

## $faint
##   v1 c1 v2 c2 dH dV dC f.case1 f.case2 f.case3 low.value.chroma       res
## 1  6  3  3  4  2  3  1   FALSE   FALSE   FALSE            FALSE Prominent
## 
## $distinct
##   v1 c1 v2 c2 dH dV dC d.case1 d.case2 d.case3       res
## 1  6  3  3  4  2  3  1   FALSE   FALSE   FALSE Prominent

7.5YR 3/3 vs 7.5YR 4/4

contrastClass(v1=3, c1=3, v2=4, c2=4, dH=0, dV=1, dC=1, verbose = TRUE)

## $faint
##   v1 c1 v2 c2 dH dV dC f.case1 f.case2 f.case3 low.value.chroma   res
## 1  3  3  4  4  0  1  1    TRUE   FALSE   FALSE            FALSE Faint
## 
## $distinct
##   v1 c1 v2 c2 dH dV dC d.case1 d.case2 d.case3   res
## 1  3  3  4  4  0  1  1   FALSE   FALSE   FALSE Faint

10YR 2/2 vs 2.5YR 2/2

contrastClass(v1=2, c1=2, v2=2, c2=2, dH=0, dV=0, dC=0, verbose = TRUE)

## $faint
##   v1 c1 v2 c2 dH dV dC f.case1 f.case2 f.case3 low.value.chroma   res
## 1  2  2  2  2  0  0  0    TRUE   FALSE   FALSE             TRUE Faint
## 
## $distinct
##   v1 c1 v2 c2 dH dV dC d.case1 d.case2 d.case3   res
## 1  2  2  2  2  0  0  0   FALSE   FALSE   FALSE Faint

7.5YR 3/4 vs 7.5YR 6/3

contrastClass(v1=3, c1=4, v2=5, c2=3, dH=0, dV=3, dC=1, verbose = TRUE)

## $faint
##   v1 c1 v2 c2 dH dV dC f.case1 f.case2 f.case3 low.value.chroma      res
## 1  3  4  5  3  0  3  1   FALSE   FALSE   FALSE            FALSE Distinct
## 
## $distinct
##   v1 c1 v2 c2 dH dV dC d.case1 d.case2 d.case3      res
## 1  3  4  5  3  0  3  1    TRUE   FALSE   FALSE Distinct

2.2.3 Hue Position

Differences in Munsell hue (\(dH\)) are computed according to radial position (clock-wise) from 5R to 5PB.

par(mar = c(0, 0, 0, 0))
huePositionCircle()

Use huePosition() to extract hues in order, or convert a vector of hues into an index of positions.

# hues used in the description of soil color (except for neutral hues)
hues <- huePosition(x = NULL, returnHues = TRUE)
# convert hue names into position
hue.pos <- huePosition(hues)
# check
kable(head(cbind(hue.pos, hues)))

hue.pos	hues
1	5R
2	7.5R
3	10R
4	2.5YR
5	5YR
6	7.5YR

Arranged in CIELAB colorspace, the ordering of hues looks like this. The huePositionPlot() function in {sharpshootR} will make this kind of figure for any combination of value and chroma.

# make these into real colors by fixing value and chroma at '6'
colors <- paste0(hues, ' 6/6')
# convenient labels for figure
hue.labels <- sprintf("%s\n%s", hues, hue.pos)

# combine colors names, hex representation of colors, and CIELAB coordinates
x <- data.frame(
  colors, 
  hex=parseMunsell(colors), 
  parseMunsell(colors, returnLAB=TRUE), 
  stringsAsFactors = FALSE
)

# make the figure
plot(B ~ A, data=x, type='n', asp=0.75, xlab='A-coordinate', ylab='B-coordinate', ylim=c(-25, 45), main='Hue Order per TN #2\nCIELAB Colorspace')
abline(h=0, v=0, lty=2)
points(B ~ A, data=x, pch=15, col=x$hex, cex=4.5)
text(x$A, x$B, labels = hue.labels, cex=0.75)

Demo sharpshootR::huePostitionPlot()

sharpshootR::huePositionPlot(contour.dE00 = TRUE, origin = '10YR 4/4')

2.3 Color Contrast Charts

Inspect the data used to generate the figures with the additional argument returnData=TRUE.

Basic chart, single reference hue.

contrastChart(m = '7.5YR 4/3', hues = c('7.5YR'))

Threshold color chips by \(\Delta{E_{00}}\).

contrastChart(m = '7.5YR 4/3', hues = c('7.5YR'), thresh = 15)

Multiple reference hues.

contrastChart(m = '7.5YR 4/3', hues = c('5YR', '7.5YR', '10YR'))

Split chips by contrast class, single reference hue.

contrastChart(m = '7.5YR 4/3', hues = c('7.5YR'), style = 'CC')

Split chips by contrast class, multiple reference hues. This works better when run through latticeExtra::useOuterStrips().

fig <- contrastChart(m = '7.5YR 4/3', hues = c('5YR', '7.5YR'), style = 'CC')
useOuterStrips(fig, strip = strip.custom(bg=grey(0.85)), strip.left = strip.custom(bg=grey(0.85)) )

2.4 Example 1

# load full LUT from {aqp}
data(munsell)

# subset some common colors
hues <- c('10R', '2.5YR', '5YR', '7.5YR', '10YR', '2.5Y', '5Y')
x <- subset(munsell, subset=value %in% 3:5 & chroma %in% 1:3 & hue %in% hues)

# convert into an ordered factor, according to hue position
huePos <- huePosition(x=NULL, returnHues = TRUE)
x$hue <- factor(x$hue, levels=huePos, ordered = TRUE)
x$hue <- droplevels(x$hue)


# ok
table(x$hue)

## 
##   10R 2.5YR   5YR 7.5YR  10YR  2.5Y    5Y 
##     9     9     9     9     9     9     9

nrow(x)

## [1] 63

# convert into hex notation for plotting
x$color <- munsell2rgb(x$hue, x$value, x$chroma)
x$munsell <- sprintf("%s %s/%s", x$hue, x$value, x$chroma)

Arrange colors in a familiar manner.

xyplot(value ~ chroma | hue, data=x,
       xlim=c(0.5, 3.5), ylim=c(2.5, 5.5), 
       scales=list(alternating=3, x=list(at=1:3), y=list(at=3:5)),
       as.table=TRUE, strip=strip.custom(bg='grey'),
       subscripts=TRUE,
       panel=function(xx, yy, subscripts, ...) {
         panel.grid(-1, -1)
         panel.points(xx, yy, col=x$color[subscripts], pch=15, cex=6)
       }
)

Evaluate \(\Delta~E_{00}\) via {farver} pacakge.

# dE00
# pair-wise distances are only present in the upper triangle
# lower triangle is 0
d <- compare_colour(x[, c('L', 'A', 'B')], from_space='lab', white_from = 'D65', method='cie2000')
# copy color codes for convenience
dimnames(d) <- list(x$munsell, x$munsell)

# convert to dist object, which expects distances in the lower triangle
d.dist <- as.dist(t(d))
# hierarchical clustering
h <- as.phylo(as.hclust(diana(d.dist)))
# nMDS
mds <- MASS::sammon(d.dist, trace = FALSE)

# remove 0's in lower triangle and diagonal
d[lower.tri(d, diag = TRUE)] <- NA
hist(d)

# double-check that ordering is correct
# yes
table(h$tip.label == x$munsell)

## 
## TRUE 
##   63

par(mar=c(0.5,0.5,1,0.5), bg='black', fg='white')

plot(h, label.offset=0.55, edge.color='white', tip.color='white', cex=0.66, direction='downwards', font=1)
tiplabels(pch=15, cex=1.5, col=x$color, offset=0.23)
title('Divisive Hierarchical Clustering', col.main='white', line=0)
mtext(text = expression("Distance metric:"~Delta*E['00']), side = 3, adj = 0, line=-0.5)

plot(h, edge.color='white', tip.color='white', cex=0.5, direction='downwards', show.tip.label=FALSE)
tiplabels(pch=15, cex=1.5, col=x$color, offset=0.25)
tiplabels(text=x$value, cex=0.55, frame='none', col='white', offset=0.23)
title('Divisive Hierarchical Clustering', col.main='white', line=0)
mtext(text = expression("Distance metric:"~Delta*E['00']), side = 3, adj = 0, line=-0.5)

par(mar=c(0.5,0.5,1.5,0.5), bg='black', fg='white')

plot(mds$points, type='n', axes=FALSE)
abline(h=0, v=0, col=grey(0.85), lty=3)
points(mds$points, pch=15, col=x$color, cex=3.5)
text(mds$points[, 1],  mds$points[, 2], x$value, cex=0.75)
mtext(text = expression("Distance metric:"~Delta*E['00']), side = 3, adj = 0)
mtext(text = 'Munsell value', side = 1, adj = 1, line=-1)
title('nMDS', col.main='white', line=0.75)

plot(mds$points, type='n', axes=FALSE)
abline(h=0, v=0, col=grey(0.85), lty=3)
points(mds$points, pch=15, col=x$color, cex=3.5)
text(mds$points[, 1],  mds$points[, 2], paste0(x$value, '/', x$chroma), cex=0.66)
mtext(text = expression(Delta*E['00']), side = 3, adj = 0)
mtext(text = 'Munsell value/chroma', side = 1, adj = 1, line=-1)
title('nMDS', col.main='white', line=0.75)

plot(mds$points, type='n', axes=FALSE)
abline(h=0, v=0, col=grey(0.85), lty=3)
points(mds$points, pch=15, col=x$color, cex=3.5)
text(mds$points[, 1],  mds$points[, 2], x$hue, cex=0.66)
mtext(text = expression(Delta*E['00']), side = 3, adj = 0)
mtext(text = 'Munsell hue', side = 1, adj = 1, line=-1)
title('nMDS', col.main='white', line=0.75)

plot(mds$points, type='n', axes=FALSE)
abline(h=0, v=0, col=grey(0.85), lty=3)
points(mds$points, pch=15, col=x$color, cex=3.5)
text(mds$points[, 1],  mds$points[, 2], x$munsell, cex=0.66)
mtext(text = expression(Delta*E['00']), side = 3, adj = 0)
mtext(text = 'Munsell color', side = 1, adj = 1, line=-1)
title('nMDS', col.main='white', line=0.75)

2.4.1 Nearly Imperceptible Differences

TODO

2.5 Color Chips Arranged by Perceptual Differences

# subset some common colors
hues <- c('7.5YR')
x <- subset(munsell, subset=value %in% 3:8 & chroma %in% c(1,2,3,4,6,8) & hue %in% hues)

# convert into hex notation for plotting
x$color <- munsell2rgb(x$hue, x$value, x$chroma)
x$munsell <- sprintf("%s/%s", x$value, x$chroma)

# dE00
# pair-wise distances are only present in the upper triangle
# lower triangle is 0
d <- compare_colour(x[, c('L', 'A', 'B')], from_space='lab', white_from = 'D65', method = 'cie2000')

# copy color codes for convenience
dimnames(d) <- list(x$munsell, x$munsell)

# convert to dist object, which expects distances in the lower triangle
d.dist <- as.dist(t(d))

# nMDS
mds <- MASS::sammon(d.dist, trace = FALSE)

# rotate 90 degrees CCW
# to roughly follow Munsell color book page layout
# column-order
m <- matrix(
  c(0, -1, 
    1, 0), 
  byrow = FALSE, ncol = 2
)

# apply transformation
mds.trans <- mds$points %*% m

par(mar=c(2.5, 0.5, 3, 0.5), bg='black', fg='white')

plot(mds.trans, type='n', axes=FALSE, asp=1)
abline(h=seq(-25, 25, 5), v=seq(-25, 25, 5), col='white', lty=3)
points(mds.trans, pch=22, bg=x$color, cex=6)
text(mds.trans[, 1],  mds.trans[, 2], paste0(x$value, '/', x$chroma), cex=0.66)
# mtext(text = expression(Delta*E['00']), side = 3, adj = 0)
# mtext(text = '7.5YR', side = 1, adj = 1, line=-1)
title('Perceptual Distances\n7.5YR Page', col.main='white', line=0.75)
mtext(text = expression(Delta*E['00']%->%nMDS%->%90*degree~CCW~rotation), side = 1, adj = 0)

Map perceptual arrangement to original representation in color book.

p <- procrustes(mds.trans, x[, c('value', 'chroma')], scale = TRUE, symmetric = TRUE)

par(mar=c(2.5, 0, 3, 0), bg='black', fg='white')

plot(p, pch=22, kind=0, axes=FALSE, xlab='', ylab='')
points(p, display='rotated', pch=22, bg=x$color, cex=5)
points(p, display='target', pch=21, bg=x$color, cex=4)
lines(p, type='arrows', length=0.1)
text(p, display='rotated', cex=0.66, labels=paste0(x$value, '/', x$chroma))
title('Perceptual vs Color Book Arrangement', col.main='white', line=1)
title(sub='7.5YR Page', col.sub='white', line=0)
legend('bottomright', legend=c('color book', 'perceptual'), pch=c(22, 21), pt.bg=parseMunsell('7.5YR 4/6'), pt.cex=2, bty='n', inset=c(0.25, 0))

2.6 Moist vs. Dry Colors

x <- fetchOSD('musick', colorState = 'dry')
y <- fetchOSD('musick', colorState = 'moist')

m1 <- sprintf("%s %s/%s", x$hue, x$value, x$chroma)
m2 <- sprintf("%s %s/%s", y$hue, y$value, y$chroma)

cc <- colorContrast(m1, m2)

colorContrastPlot(m1, m2, labels=c('Dry', 'Moist'), d.cex = 0.9)

Tinker with plot settings and optimize for dark background.

x <- fetchOSD('Leefield', colorState = 'dry')
y <- fetchOSD('Leefield', colorState = 'moist')

m1 <- sprintf("%s %s/%s", x$hue, x$value, x$chroma)
m2 <- sprintf("%s %s/%s", y$hue, y$value, y$chroma)

par(bg='black', fg='white')
colorContrastPlot(m1, m2, labels=c('Dry', 'Moist'), printMetrics = TRUE, d.cex = 1.25, col.cex = 1.5, label.cex = 2, label.font = 2)

2.7 Differences Between Adjacent Chips

# fresh start
data(munsell)

# for now, most commonly used hue pages, but expanded to full range in value / chroma
x <- subset(munsell, subset = hue %in% c('10R', '2.5YR', '5YR', '7.5YR', '10YR', '2.5Y', '5Y') & value %in% 2:8 & chroma %in% 2:8)

# create all pair-wise combinations
cols <- paste0(x$hue, ' ', x$value, '/', x$chroma)
z <- combn(cols, 2)

# convert to 2-column data.frame
# ~ 674,541 combinations!
z <- data.frame(t(z), stringsAsFactors = FALSE)

## TODO: make this a new function

# borrowed from aqp::colorContrast
m1.pieces <- parseMunsell(z[[1]], convertColors = FALSE)
m2.pieces <- parseMunsell(z[[2]], convertColors = FALSE)

# convert to value and chroma to numeric
m1.pieces[[2]] <- as.numeric(m1.pieces[[2]])
m1.pieces[[3]] <- as.numeric(m1.pieces[[3]])
m2.pieces[[2]] <- as.numeric(m2.pieces[[2]])
m2.pieces[[3]] <- as.numeric(m2.pieces[[3]])

# difference in number of hue chips, clock-wise, as specified in:
# https://www.nrcs.usda.gov/wps/portal/nrcs/detail/soils/ref/?cid=nrcs142p2_053569
dH <- abs(huePosition(m1.pieces[[1]]) - huePosition(m2.pieces[[1]]))
# difference in number of value chips
dV <- abs(m1.pieces[[2]] - m2.pieces[[2]])
# difference in number of chroma chips
dC <- abs(m1.pieces[[3]] - m2.pieces[[3]])

## ^^^^ need a new function

# just adjacent chips
# ~ 11,498 combinations
idx <- which(dH %in% c(0, 1) & dV %in% c(0, 1) & dC %in% c(0, 1))
z <- z[idx, ]

# compute color contrast metrics on adjacent chips
d <- colorContrast(z[[1]], z[[2]])

# summaries
quantile(d$dE00)

##         0%        25%        50%        75%       100% 
##  0.5991672  4.5772217  8.5009734  9.7734551 12.6498032

# encode combinations
d$code <- sprintf("%s-%s-%s", d$dH, d$dV, d$dC)

# compute select quantiles of dE00 by combination
adj <- lapply(split(d, d$code), function(i) {
  qq <- round(quantile(i$dE00, probs = c(0.05, 0.5, 0.95)))
  data.frame(
    dH = i$dH[1],
    dV = i$dV[1],
    dC = i$dC[1],
    q05 = qq[1],
    q50 = qq[2],
    q95 = qq[3]
  )
    
}
)

adj <- do.call('rbind', adj)

# all combinations
knitr::kable(adj)

	dH	dV	dC	q05	q50	q95
0-0-1	0	0	1	1	3	4
0-1-0	0	1	0	7	9	10
0-1-1	0	1	1	7	9	11
1-0-0	1	0	0	2	4	5
1-0-1	1	0	1	4	5	6
1-1-0	1	1	0	8	9	11
1-1-1	1	1	1	8	10	11

# just changes along a single dimension 
idx <- which(row.names(adj) %in% c('1-0-0', '0-1-0', '0-0-1'))
knitr::kable(adj[idx[3:1], ])

	dH	dV	dC	q05	q50	q95
1-0-0	1	0	0	2	4	5
0-1-0	0	1	0	7	9	10
0-0-1	0	0	1	1	3	4

2.8 dE00 Between Hues

Including N chip.

library(corrplot)

# fresh start
data(munsell)

# for now, most commonly used hue pages, but expanded to full range in value / chroma
x1 <- subset(munsell, subset = value == 6 & chroma == 6)

# neutral
x2 <- subset(munsell, subset = hue == 'N' & value == 6 )

# combine
x <- rbind(x1, x2)

# re-arrange into "hue order"
idx <- match(huePosition(returnHues = TRUE, includeNeutral = TRUE), x$hue)
x <- x[idx, ]

x$m <- sprintf("%s %s/%s", x$hue, x$value, x$chroma)
x$col <- rgb(x$r, x$g, x$b, maxColorValue = 1)


d <- farver::compare_colour(
  from = x[, c('L', 'A', 'B')],
  to = x[, c('L', 'A', 'B')], 
  from_space = 'lab', 
  method = 'CIE2000'
)

dimnames(d) <- list(x$m, x$m)

par(fg = 'white', bg = 'black')
corrplot(
    d, 
    col = mako(25), 
    is.corr = FALSE, 
    col.lim = range(d), 
    method = "color", 
    order = "original",
    type = "upper",
    diag = FALSE,
    tl.cex = 0.66,
    tl.col = x$col,
    # tl.col = 'white',
    mar = c(0.1, 0, 0, 0.8),
    # addgrid = TRUE,
  )

Without N chip.

x <- x[x$hue != 'N', ]

d <- farver::compare_colour(
  from = x[, c('L', 'A', 'B')],
  to = x[, c('L', 'A', 'B')], 
  from_space = 'lab', 
  method = 'CIE2000'
)

dimnames(d) <- list(x$m, x$m)

par(fg = 'white', bg = 'black')
corrplot(
    d, 
    col = mako(25), 
    is.corr = FALSE, 
    col.lim = range(d), 
    method = "color", 
    order = "original",
    type = "upper",
    diag = FALSE,
    tl.cex = 0.66,
    tl.col = x$col,
    # tl.col = 'white',
    mar = c(0.1, 0, 0, 0.8),
    # addgrid = TRUE,
  )

2.9 Example 2

data(munsell)
x1 <- subset(munsell, subset = value %in% 2:6 & chroma %in% c(0, 2:6) & hue != 'N')
x2 <- subset(munsell, subset = hue == 'N')


cols1 <- paste0(x1$hue, ' ', x1$value, '/', x1$chroma)
cols2 <- paste0(x2$hue, ' ', x2$value, '/', x2$chroma)

z <- expand.grid(colors = cols1, neutral = cols2)

d <- colorContrast(z$neutral, z$colors)

# extract hue
.m1 <- parseMunsell(d$m1, convertColors = FALSE)
.m2 <- parseMunsell(d$m2, convertColors = FALSE)

d$.hue1 <- factor(.m1$hue, levels = huePosition(returnHues = TRUE, includeNeutral = TRUE))
d$.hue2 <- factor(.m2$hue, levels = huePosition(returnHues = TRUE, includeNeutral = TRUE))

d$.value1 <- .m1$value
d$.value2 <- .m2$value

bwplot(.hue2 ~ dE00, data = d, 
       par.settings = tactile.theme(), 
       xlab = expression(Delta*E['00']),
       main = 'CIE Delta-E00 vs. Soil Color Contrast Class\nhues: 10R-5Y, N; value 2-8; chroma 2-8, 0',
       sub = 'Soil Survey Technical Note #2',
       scales = list(x = list(tick.number = 10)), 
       panel = function(...) {
         panel.grid(h = 3, v = -1)
         panel.bwplot(...)
       })

bwplot(.value2  ~ dE00, data = d, 
       par.settings = tactile.theme(), 
       xlab = expression(Delta*E['00']),
       main = 'CIE Delta-E00 vs. Soil Color Contrast Class\nhues: 10R-5Y, N; value 2-8; chroma 2-8, 0',
       sub = 'Soil Survey Technical Note #2',
       scales = list(x = list(tick.number = 10)), 
       panel = function(...) {
         panel.grid(h = 3, v = -1)
         panel.bwplot(...)
       })

bwplot(.value2  ~ dE00 | .hue2, data = d, 
       par.settings = tactile.theme(), 
       xlab = expression(Delta*E['00']),
       main = 'CIE Delta-E00 vs. Soil Color Contrast Class\nhues: 10R-5Y, N; value 2-8; chroma 2-8, 0',
       sub = 'Soil Survey Technical Note #2',
       scales = list(x = list(tick.number = 10)), 
       panel = function(...) {
         panel.grid(h = 3, v = -1)
         panel.bwplot(...)
       })

levelplot(
  dE00 ~ .value2 * .hue2, data = d,
  subset = .value1 == 2,
  par.settings = tactile.theme()
  
)

bwplot(cc ~ dE00, data = d, 
       par.settings = tactile.theme(), 
       xlab = expression(Delta*E['00']),
       main = 'CIE Delta-E00 vs. Soil Color Contrast Class\nhues: 10R-5Y, N; value 2-8; chroma 2-8, 0',
       sub = 'Soil Survey Technical Note #2',
       scales = list(x = list(tick.number = 10)), 
       panel = function(...) {
         panel.grid(h = 3, v = -1)
         panel.bwplot(...)
       })

library(rms)
library(tactile)

data(munsell)
x <- subset(munsell, subset = value %in% 2:8 & chroma %in% c(0, 2:8) & hue %in% c('10R', '2.5YR', '5YR', '7.5YR', '10YR', '2.5Y', '5Y', '10GY', 'N'))

## TODO: more expansive set
# x <- subset(munsell, subset = value %in% 2:6 & chroma %in% c(0, 2:6))

cols <- paste0(x$hue, ' ', x$value, '/', x$chroma)
z <- combn(cols, 2)

z <- data.frame(t(z), stringsAsFactors = FALSE)
d <- colorContrast(z[[1]], z[[2]])

# extract hue
.m1 <- parseMunsell(d$m1, convertColors = FALSE)
.m2 <- parseMunsell(d$m2, convertColors = FALSE)

d$.hue1 <- factor(.m1$hue, levels = huePosition(returnHues = TRUE, includeNeutral = TRUE))
d$.hue2 <- factor(.m2$hue, levels = huePosition(returnHues = TRUE, includeNeutral = TRUE))

kable(head(d), row.names = FALSE)

m1	m2	dC	dE00	cc	.hue1	.hue2
10GY 2/2	10GY 2/3	1	4.690099	Faint	10GY	10GY
10GY 2/2	10GY 2/4	2	8.228123	Distinct	10GY	10GY
10GY 2/2	10GY 2/5	3	11.022977	Distinct	10GY	10GY
10GY 2/2	10GY 2/6	4	13.297008	Prominent	10GY	10GY
10GY 2/2	10GY 2/7	5	14.568608	Prominent	10GY	10GY
10GY 2/2	10GY 2/8	6	15.457711	Prominent	10GY	10GY

bwplot(cc ~ dE00, data = d, 
       par.settings = tactile.theme(), 
       xlab = expression(Delta*E['00']),
       main = 'CIE Delta-E00 vs. Soil Color Contrast Class\nhues: 10R-5Y, N; value 2-8; chroma 2-8, 0',
       sub = 'Soil Survey Technical Note #2',
       scales = list(x = list(tick.number = 10)), 
       panel = function(...) {
         panel.grid(h = 3, v = -1)
         panel.bwplot(...)
       })

d.sub <- subset(d, subset = .hue1 == '10YR' | .hue2 == '10YR')

bwplot(dE00 ~ .hue1, data = d.sub, 
       as.table = TRUE,
       par.settings = tactile.theme(), 
       xlab = expression(Delta*E['00']),
       main = 'CIE Delta-E00 vs. Soil Color Contrast Class\nhues: 10R-5Y, N; value 2-8; chroma 2-8, 0',
       sub = 'Soil Survey Technical Note #2',
       scales = list(alternating = 3, x = list(tick.number = 10)), 
       panel = function(...) {
         panel.grid(h = 3, v = -1)
         panel.bwplot(...)
       })

# only viz comparisons between colors with same hue
d.sub <- subset(d, subset = .hue1 == .hue2)

# panel by hue
bwplot(cc ~ dE00 | .hue1, data = d.sub, 
       as.table = TRUE,
       par.settings = tactile.theme(), 
       xlab = expression(Delta*E['00']),
       main = 'CIE Delta-E00 vs. Soil Color Contrast Class\nhues: 10R-5Y, value 2-8, chroma 2-8',
       sub = 'Soil Survey Technical Note #2',
       scales = list(alternating = 1, x = list(tick.number = 10)), 
       panel = function(...) {
         panel.grid(h = 3, v = -1)
         panel.bwplot(...)
       })

bwplot(dV ~ dE00, data=d, 
       subset = dH < 1 & dC < 1,
       par.settings=tactile.theme(), 
       xlab=expression(Delta*E['00']),
       main='CIE Delta-E00 vs. Change in Munsell Value\nhue 10R-5Y, value 2-8, chroma 2-8',
       sub = 'Constant Hue and Chroma',
       scales=list(x=list(tick.number=10)), 
       panel=function(...) {
         panel.grid(h=3, v=-1)
         panel.bwplot(...)
       })

bwplot(dH ~ dE00, data=d,
       subset = dV < 1 & dC < 1,
       par.settings=tactile.theme(), 
       xlab=expression(Delta*E['00']),
       main='CIE Delta-E00 vs. Change in Munsell Hue\nhue 10R-10R-5Y, value 2-8, chroma 2-8',
       sub = 'Constant Value and Chroma',
       scales=list(x=list(tick.number=10)), 
       panel=function(...) {
         panel.grid(h=3, v=-1)
         panel.bwplot(...)
       })

bwplot(dC ~ dE00, data=d, 
       subset = dV < 1 & dH < 1,
       par.settings=tactile.theme(), 
       xlab=expression(Delta*E['00']),
       main='CIE Delta-E00 vs. Change in Munsell Chroma\nhue 10R-10R-5Y, value 2-8, chroma 2-8',
       sub = 'Constant Hue and Value',
       scales=list(x=list(tick.number=10)), 
       panel=function(...) {
         panel.grid(h=3, v=-1)
         panel.bwplot(...)
       })

idx.F <- which(d$cc == 'Faint' & d$dE00 > 12)
idx.D <- which(d$cc == 'Distinct' & d$dE00 > 20)
idx.P <- which(d$cc == 'Prominent' & d$dE00 < 20)


kable(head(d[idx.F, ]), row.names = FALSE)

m1	m2	dH	dV	dC	dE00	cc	.hue1	.hue2
10GY 2/2	10GY 4/2	0	2	0	16.73171	Faint	10GY	10GY
10GY 2/2	10GY 4/3	0	2	1	17.91859	Faint	10GY	10GY
10GY 2/2	10R 2/2	12	0	0	21.48346	Faint	10GY	10R
10GY 2/2	10R 3/2	12	1	0	23.08248	Faint	10GY	10R
10GY 2/2	10YR 2/2	8	0	0	14.84723	Faint	10GY	10YR
10GY 2/2	10YR 3/2	8	1	0	16.86496	Faint	10GY	10YR

# plot(dE00 ~ dH, data=d)
# plot(dE00 ~ dV, data=d)
# plot(dE00 ~ dC, data=d)

dd <- datadist(d)
options(datadist = "dd")

# there is a little bit of curvature
(m <- ols(dE00 ~ rcs(dH, 3) + rcs(dV, 3) + rcs(dC, 3), data = d))

## Linear Regression Model
##  
##  ols(formula = dE00 ~ rcs(dH, 3) + rcs(dV, 3) + rcs(dC, 3), data = d)
##  
##                     Model Likelihood    Discrimination    
##                           Ratio Test           Indexes    
##  Obs   78606    LR chi2    187346.38    R2       0.908    
##  sigma4.5380    d.f.               6    R2 adj   0.908    
##  d.f.  78599    Pr(> chi2)    0.0000    g       15.700    
##  
##  Residuals
##  
##       Min       1Q   Median       3Q      Max 
##  -27.7705  -2.4446  -0.4395   1.9460  34.7195 
##  
##  
##            Coef    S.E.   t      Pr(>|t|)
##  Intercept  8.8965 0.0590 150.83 <0.0001 
##  dH         0.8184 0.0166  49.24 <0.0001 
##  dH'        1.3236 0.0251  52.70 <0.0001 
##  dV         3.6561 0.0286 127.81 <0.0001 
##  dV'        6.0015 0.0395 151.81 <0.0001 
##  dC        -0.1174 0.0279  -4.21 <0.0001 
##  dC'        1.3793 0.0367  37.61 <0.0001 
##

# (m <- ols(dE00 ~ dH + dV + dC, data = d))

plot(Predict(m, dH=NA, dV=0, dC=0))

plot(Predict(m, dH=0, dV=NA, dC=0))

plot(Predict(m, dH=0, dV=0, dC=NA))

plot(Predict(m, dH=NA, dV=0:2, dC=0))

plot(Predict(m, dV=NA, dH=0:2, dC=0))

plot(Predict(m, dC=NA, dV=0:2, dH=0))

plot(summary(m, dH=c(0,1), dV=c(0,1), dC=c(0,1)))

plot(summary(m))

plot(anova(m), what='partial R2')

plot(nomogram(m))

Add intercept to LP

print(nomogram(m, dH=c(0, 1, 2, 3)))

## Points per unit of linear predictor: 2.068549 
## Linear predictor units per point   : 0.4834307 
## 
## 
##  dH Points
##  0  0     
##  1  2     
##  2  4     
##  3  6     
## 
## 
##  dV Points
##  0    0   
##  1    8   
##  2   19   
##  3   35   
##  4   55   
##  5   78   
##  6  100   
## 
## 
##  dC Points
##  0   0    
##  1   0    
##  2   1    
##  3   2    
##  4   5    
##  5   8    
##  6  11    
##  7  14    
##  8  18

2.10 TODO

Add N hues
include value and chroma of starting chip

library(terra)
library(splines)

# re-fit simpler model, without intercept
summary(m <- lm(dE00 ~ dH * bs(dV, 3) * dC - 1, data = d))

## 
## Call:
## lm(formula = dE00 ~ dH * bs(dV, 3) * dC - 1, data = d)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -28.246  -2.072  -0.458   1.398  33.657 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## dH                3.168009   0.012710 249.252  < 2e-16 ***
## bs(dV, 3)1        8.405325   0.190675  44.082  < 2e-16 ***
## bs(dV, 3)2       36.782454   0.269613 136.427  < 2e-16 ***
## bs(dV, 3)3       59.832350   0.164462 363.806  < 2e-16 ***
## dC                2.631444   0.018068 145.638  < 2e-16 ***
## dH:bs(dV, 3)1    -0.913509   0.049513 -18.450  < 2e-16 ***
## dH:bs(dV, 3)2    -2.273132   0.055273 -41.125  < 2e-16 ***
## dH:bs(dV, 3)3    -2.276508   0.036709 -62.015  < 2e-16 ***
## dH:dC            -0.322343   0.005315 -60.646  < 2e-16 ***
## bs(dV, 3)1:dC    -0.919512   0.076476 -12.024  < 2e-16 ***
## bs(dV, 3)2:dC    -2.349342   0.089353 -26.293  < 2e-16 ***
## bs(dV, 3)3:dC    -2.337495   0.058495 -39.961  < 2e-16 ***
## dH:bs(dV, 3)1:dC  0.100799   0.018526   5.441 5.32e-08 ***
## dH:bs(dV, 3)2:dC  0.281085   0.019078  14.734  < 2e-16 ***
## dH:bs(dV, 3)3:dC  0.275420   0.013151  20.942  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.071 on 78591 degrees of freedom
## Multiple R-squared:  0.9841, Adjusted R-squared:  0.9841 
## F-statistic: 3.253e+05 on 15 and 78591 DF,  p-value: < 2.2e-16

# RMSE
# ~ 2.6
sqrt(mean((d$dE00 - predict(m))^2))

## [1] 4.07097

nd <- expand.grid(dH = 0:3, dV = 0:4, dC = 0:4)
p <- data.frame(predict(m, newdata = nd, interval = 'prediction', level = 0.8))
nd$dE00 <- p$fit
nd$interval <- p$upr - p$lwr 

# factors
nd$dH <- factor(nd$dH, labels = sprintf("dH = %s", sort(unique(nd$dH))))


head(nd)

##       dH dV dC     dE00 interval
## 1 dH = 0  0  0 0.000000 10.43540
## 2 dH = 1  0  0 3.168009 10.43545
## 3 dH = 2  0  0 6.336018 10.43560
## 4 dH = 3  0  0 9.504027 10.43586
## 5 dH = 0  1  0 5.749854 10.43628
## 6 dH = 1  1  0 8.432277 10.43605

tps <- tactile.theme(regions = list(col = viridis(200, alpha = 0.85, direction = -1)))

levelplot(
  dE00 ~ dC * dV | dH, 
  data = nd,
  main = expression(Estimated~Delta*E['00']),
  sub = 'Common Soil Hues, Value 2-8, Chroma 2-8',
  ylab = expression(Delta*Value),
  xlab = expression(Delta*Chroma),
  as.table = TRUE,
  par.settings = tps,
  scales = list(alternating = 1),
  subscripts = TRUE,
  panel = function(x, y, z, subscripts, ...) {
    # background colors / shading ~ dE00
    panel.levelplot(x = x, y = y, z = z, subscripts = subscripts, ...)
    
    
    # labels for dE00 values
    # convert NA -> ''
    .NAidx <- which(is.na(z[subscripts]))
    .txt <- format(round(z[subscripts]), trim = TRUE)
    .txt[.NAidx] <- ''
    
    # legible labels based reasonable threshold 
    .col <- ifelse(z[subscripts] > 25, 'white', 'black')
    
    # annotate with dE00
    panel.text(
      x = x[subscripts],
      y = y[subscripts],
      labels = .txt,
      col = .col,
      cex = 0.8,
      font = 2
    )
    
    # mark regions of current contrast classes
    # polygon outlines perhaps
    .n <- length(x[subscripts])
    .bogus <- rep(4, times = .n)
    
    # simple proxy for dH by panel
    # requires correct ordering of panels
    .dH <- panel.number() - 1
    
    # compute contrast class within panels
    # assume starting value and chroma are not "low"
    .cc <- contrastClass(
      v1 = .bogus, 
      c1 = .bogus, 
      v2 = .bogus, 
      c2 = .bogus, 
      dH = rep(.dH, times = .n), 
      dV = y[subscripts], 
      dC = x[subscripts]
    )
    
    panel.text(
      x = x[subscripts] + 0.33,
      y = y[subscripts] - 0.33,
      labels = abbreviate(.cc, minlength = 1),
      col = .col,
      cex = 0.66,
      font = 3
    )
    
    # assemble contrast class grid / overlay
    .xyz <- cbind(x = x[subscripts], y = y[subscripts], cc = .cc)
    .r <- rast(.xyz, type = 'xyz')
    .p <- as.polygons(.r, dissolve = TRUE)
    
    # polygon coordinates
    .coords <- crds(.p)
    
    # debug
    # browser()    
    
    # annotate contrast class boundaries
    panel.polygon(.coords, lwd = 1, lty = 3, border = 'white')
    
    # could have done this manually with 4 line segments per panel...
    
  }
)

## TODO: think about how to use / display prediction interval
levelplot(
  interval ~ dC * dV | dH, 
  data = nd,
  main = expression(Estimated~Delta*E['00']),
  ylab = expression(Delta*Value),
  xlab = expression(Delta*Chroma),
  as.table = TRUE,
  par.settings = tps,
  scales = list(alternating = 1)
)

Hmmm.

d.m <- parseMunsell(d$m1, convertColors = FALSE)

d$v1 <- d.m$value

xyplot(
  dE00 ~ dV | factor(v1),
  groups = factor(dC),
  data = d, 
  subset = dH == 1,
  par.settings = tactile.theme(), 
  scales = list(alternating = 1),
  type = c('smooth', 'p')
)

du <- unique(d[, c('dH', 'dV', 'dC', 'cc')])
du <- du[which(du$dH < 4 & du$dV < 4 & du$dC < 5), ]

du$dH <- factor(du$dH, labels = sprintf("dH = %s", sort(unique(du$dH))))

tps <- tactile.theme(regions = list(col = grey.colors))


levelplot(
  cc ~ dV * dC | dH, 
  data = du,
  as.table = TRUE,
  par.settings = tps,
  scales = list(alternating = 1),
  subscripts = TRUE,
  panel = function(x, y, z, subscripts, ...) {
    panel.levelplot(x = x, y = y, z = z, subscripts = subscripts, ...)
    
    .NAidx <- which(is.na(z[subscripts]))
    .txt <- abbreviate(z[subscripts], minlength = 1)

    panel.text(x[subscripts], y[subscripts], labels = .txt, col = 'black', cex = 0.85, font = 2)
  }
)

This document is based on aqp version 2.0.

Soil Color Contrast

D.E. Beaudette

2023-02-22