Introduction

# load required libaries
library(aqp)
library(soilDB)
library(sharpshootR)
library(cluster)
library(lattice)
library(reshape2)
library(colorspace)

# init vector of taxonnames to summarize
soils <- c('Whiterock', 'Copperopolis', 'Amador', 'Auburn', 'Loafercreek', 'Argonaut', 'Crimeahouse')

Optionally re-make cached data. Requires NASIS and relevant soils in the selected set.

# get all pedons from the selected set
x <- fetchNASIS(from='pedons', rmHzErrors = FALSE, nullFragsAreZero = FALSE)

# convert vector of taxonnames into REGEX pattern for matching
pat <- paste0(soils, collapse = '|')

# subset pedons that match our REGEX pattern
idx <- grep(pat, x$taxonname, ignore.case = TRUE)
x <- x[idx, ]

# save for later
save(x, file='profile-summary-SPC.rda')

Get the cached copy of the SoilProfileCollection object used in this tutorial.

# load cached copy
# load(url('http://ncss-tech.github.io/AQP/aqp/profile-summary-SPC.rda'))
load(file='profile-summary-SPC.rda')

# normalize taxonname via REGEX
for(i in soils) {
  x$taxonname[grep(i, x$taxonname, ignore.case = TRUE)] <- i
}
  

# tabulate number of profiles
table(x$taxonname)
## 
##       Amador     Argonaut       Auburn Copperopolis  Crimeahouse  Loafercreek    Whiterock 
##           15           31           63           27           29          109           34
# convert soil colors to CIE LAB colorspace
# moist colors
cols.lab <- munsell2rgb(x$m_hue, x$m_value, x$m_chroma, returnLAB = TRUE)
x$m_L <- cols.lab$L
x$m_A <- cols.lab$A
x$m_B <- cols.lab$B

# dry colors
cols.lab <- munsell2rgb(x$d_hue, x$d_value, x$d_chroma, returnLAB = TRUE)
x$d_L <- cols.lab$L
x$d_A <- cols.lab$A
x$d_B <- cols.lab$B


# aggregate data by normalized taxonname, via slice-wise median
a.colors <- slab(x, taxonname ~ d_r + d_g + d_b + d_L + d_A + d_B + clay + phfield + total_frags_pct)

# throw out aggregate data that are deeper than 150cm
a.colors <- subset(a.colors, subset=bottom < 150)

# convert long -> wide format
x.colors <- dcast(a.colors, taxonname + top + bottom ~ variable, value.var = 'p.q50')

# check
head(a.colors)
##   variable taxonname      p.q5     p.q25     p.q50     p.q75     p.q95 top bottom
## 1      d_r    Amador 0.5087045 0.5695402 0.6617205 0.6617205 0.7136251   0      1
## 2      d_r    Amador 0.5087045 0.5695402 0.6617205 0.6617205 0.7136251   1      2
## 3      d_r    Amador 0.5087045 0.5695402 0.6617205 0.6617205 0.7136251   2      3
## 4      d_r    Amador 0.5087045 0.5695402 0.6617205 0.6617205 0.7136251   3      4
## 5      d_r    Amador 0.5087045 0.5695402 0.6617205 0.6617205 0.7136251   4      5
## 6      d_r    Amador 0.4205719 0.5589675 0.6617205 0.6617205 0.7136251   5      6
##   contributing_fraction
## 1             0.7333333
## 2             0.7333333
## 3             0.7333333
## 4             0.7333333
## 5             0.7333333
## 6             0.7333333
head(x.colors)
##   taxonname top bottom       d_r       d_g       d_b      d_L      d_A      d_B clay phfield
## 1    Amador   0      1 0.6617205 0.5656823 0.4098928 61.72703 4.870262 19.10015   17     5.6
## 2    Amador   1      2 0.6617205 0.5656823 0.4098928 61.72703 4.870262 19.10015   17     5.6
## 3    Amador   2      3 0.6617205 0.5656823 0.4098928 61.72703 4.870262 19.10015   17     5.6
## 4    Amador   3      4 0.6617205 0.5656823 0.4098928 61.72703 4.870262 19.10015   17     5.6
## 5    Amador   4      5 0.6617205 0.5656823 0.4098928 61.72703 4.870262 19.10015   17     5.6
## 6    Amador   5      6 0.6617205 0.5656823 0.4098928 61.72703 4.870262 19.10015   17     5.6
##   total_frags_pct
## 1               7
## 2               7
## 3               7
## 4               7
## 5               7
## 6               7
# composite sRGB triplets into an R-compatible color
# note that missing colors must be padded with NA
x.colors$soil_color <- NA
not.na <- which(complete.cases(x.colors[, c('d_L', 'd_A', 'd_B')]))
cols.srgb <- data.frame(coords(as(LAB(x.colors$d_L, x.colors$d_A, x.colors$d_B), 'sRGB')))

x.colors$soil_color[not.na] <- with(cols.srgb[not.na, ], rgb(R, G, B, maxColorValue = 1))

# aggregate bedrock depth probability by taxonname
# at 90% level of confidence
dp <- aggregateSoilDepth(x, groups = 'taxonname', crit.prob=0.9)

# init a new SoilProfileCollection from aggregate data
depths(x.colors) <- taxonname ~ top + bottom
# join-in our depth to contact data
site(x.colors) <- dp

# generate index for new ordering
new.order <- match(c('Whiterock', 'Copperopolis', 'Amador', 'Auburn', 'Loafercreek', 'Argonaut', 'Crimeahouse'), profile_id(x.colors))

# make a table of data for adding depth brackets to profile sketches
# see ?addBracket for details
depth.brackets <- site(x.colors)
names(depth.brackets)[2:3] <- c('top', 'bottom')

# add a label
depth.brackets$label <- 'P(soil >= 90%)'
par(mar=c(1, 0, 3, 1.5))

plotSPC(
  x.colors, 
  divide.hz = FALSE, 
  name = NA, 
  plot.order = new.order, 
  col.label = 'Soil Color',
  lwd = 1.25,
  axis.line.offset = -3,
  cex.depth.axis = 1,
  cex.id = 1
)

addBracket(depth.brackets, col='black', label.cex=0.85)

title('Aggregate Soil Properties')

d <- profile_compare(x.colors, vars = c('d_L', 'd_A', 'd_B'), max_d = 150, k = 0)

plotProfileDendrogram(
  x.colors, 
  clust = diana(d), 
  scaling.factor = 1, 
  name = NA, 
  divide.hz = FALSE,
  y.offset = 10,
  width = 0.2
)

addBracket(depth.brackets, col='black', label.cex=0.6)

title('Aggregate Soil Properties')

par(mar=c(1,0,3,0))

plotSPC(x.colors, divide.hz=FALSE, color='clay', name=NA, plot.order=new.order, col.label='Clay Content (%)', lwd=1.25, axis.line.offset=-3, cex.depth.axis=1, cex.id=1)
addBracket(depth.brackets, col='black', label.cex=0.85)

plot(x.colors, divide.hz=FALSE, color='phfield', name=NA, plot.order=new.order, col.label='pH', lwd=1.25, axis.line.offset=-3, cex.depth.axis=1, cex.id=1)
addBracket(depth.brackets, col='black', label.cex=0.85)

plot(x.colors, divide.hz=FALSE, color='total_frags_pct', name=NA, plot.order=new.order, col.label='Total Rock Fragment Volume (%)', lwd=1.25, axis.line.offset=-3, cex.depth.axis=1, cex.id=1)
addVolumeFraction(x.colors, 'total_frags_pct')
addBracket(depth.brackets, col='black', label.cex=0.85)

d <- profile_compare(x.colors, vars=c('clay', 'phfield', 'total_frags_pct'), max_d=150, k=0)

plotProfileDendrogram(x.colors, diana(d), scaling.factor = 1, name = NA,  divide.hz=FALSE, y.offset=10, width = 0.25)

plotProfileDendrogram(x.colors, diana(d), scaling.factor = 1, name = NA,  divide.hz=FALSE, y.offset=10, width = 0.25, color='clay')
addBracket(depth.brackets, col='black', label.cex=0.6)
addVolumeFraction(x.colors, 'total_frags_pct')

plotProfileDendrogram(x.colors, diana(d), scaling.factor = 1, name = NA,  divide.hz=FALSE, y.offset=10, width = 0.25, color='phfield')

plotProfileDendrogram(x.colors, diana(d), scaling.factor = 1, name = NA,  divide.hz=FALSE, y.offset=10, width = 0.25, color='total_frags_pct')

d <- profile_compare(x.colors, vars=c('d_L', 'clay', 'phfield', 'total_frags_pct'), max_d=150, k=0)

plotProfileDendrogram(x.colors, diana(d), scaling.factor = 1, name = NA, divide.hz=FALSE, y.offset=10, width = 0.25)
addBracket(depth.brackets, col='black', label.cex=0.6, offset = -0.35)