Soil Colors from Spectra

VNIR Spectra

library(aqp)
library(soilDB)
library(reshape2)
library(lattice)
library(tactile)

## gather some example spectra
x <- fetchRaCA('redding', get.vnir = TRUE)
x <- as.data.frame(x$spectra[1:pmin(25, nrow(x$spectra)), ])

y <- fetchRaCA('zook', get.vnir = TRUE)
y <- as.data.frame(y$spectra[1:pmin(25, nrow(y$spectra)), ])

z <- fetchRaCA('drummer', get.vnir = TRUE)
z <- as.data.frame(z$spectra[1:pmin(25, nrow(z$spectra)), ])

# combine
x <- rbind(x, y, z)

# add a temporary ID
x$.id <- 1:nrow(x)

# long format simpler to work with
m <- melt(x, id.vars = '.id')

# remove leading 'X' from wavelength
m$variable <- as.character(m$variable)
m$v <- as.numeric(gsub(pattern = 'X', replacement = '', x = m$variable))

# remove left-over columns
m$variable <- NULL

# check: ok
head(m)

##   .id    value   v
## 1   1 0.195590 350
## 2   2 0.126400 350
## 3   3 0.137835 350
## 4   4 0.157308 350
## 5   5 0.150098 350
## 6   6 0.185702 350

# subset to visible part of the spectrum that spec2Munsell can use
m.sub <- m[which(m$v >= 380 & m$v <= 730), ]

# "round" 1nm -> 10nm resolution
m.sub$v10 <- round(m.sub$v, -1)

# check: ok
head(m.sub)

##      .id    value   v v10
## 2251   1 0.162216 380 380
## 2252   2 0.082300 380 380
## 2253   3 0.096800 380 380
## 2254   4 0.120641 380 380
## 2255   5 0.109728 380 380
## 2256   6 0.149557 380 380

# aggregate to 10nm res via mean
a <- aggregate(value ~ .id + v10, data = m.sub, FUN = mean)

# factor
a$.id <- factor(a$.id)

# check: ok
head(a)

##   .id v10      value
## 1   1 380 0.16242617
## 2   2 380 0.08151667
## 3   3 380 0.09708333
## 4   4 380 0.11996617
## 5   5 380 0.10874017
## 6   6 380 0.14964150

z <- split(a, a$.id)

# TODO: may need some error handling of NA colors
z <- lapply(z, function(i) {
  # assumes 380--730nm at 10nm resolution
  sRGB <- spec2Munsell(i$value, convert = FALSE)
  # pack into df
  res <- data.frame(
    .id = i$.id[1],
    col = rgb(sRGB),
    stringsAsFactors = FALSE
  )
  #
  return(res)
})

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

###

## big assumption: colors are in the same order as levels of .id

# plot spectra
tps <- tactile.theme(
  superpose.line = list(lwd = 2, col = z$col)
)


# annotate with visible portion of the spectrum + colors
xyplot(
  value ~ v, 
  groups = .id, 
  data = m, 
  type = 'l', 
  par.settings = tps, 
  scales = list(x = list(tick.number = 20), y = list(tick.number = 10)), 
  xlab = 'Wavelength (nm)', 
  ylab = 'Reflectance', 
  panel = function(x, y, ...) {
    panel.grid(-1, -1)
    panel.abline(v = c(380, 730), lty = 2)
    panel.xyplot(x = x, y = y, ...)
    
    # position via panel geometry
    yb <- current.panel.limits()$ylim[1]
    yt <- yb + 0.01
    panel.rect(xleft = 380, xright = 435, ybottom = yb, ytop = yt, border = 'violet', col = 'violet')
    panel.rect(xleft = 435, xright = 500, ybottom = yb, ytop = yt, border = 'blue', col = 'blue')
    panel.rect(xleft = 500, xright = 520, ybottom = yb, ytop = yt, border = 'cyan', col = 'cyan')
    panel.rect(xleft = 520, xright = 565, ybottom = yb, ytop = yt, border = 'green', col = 'green')
    panel.rect(xleft = 565, xright = 590, ybottom = yb, ytop = yt, border = 'yellow', col = 'yellow')
    panel.rect(xleft = 590, xright = 625, ybottom = yb, ytop = yt, border = 'orange', col = 'orange')
    panel.rect(xleft = 625, xright = 730, ybottom = yb, ytop = yt, border = 'red', col = 'red')
    
  })