Introduction

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

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

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'))

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

# tabulate number of profiiles
table(x$taxonname)
## 
##       Amador     Argonaut       Auburn Copperopolis  Crimeahouse   Hennekenot  Loafercreek 
##           19           20          112           22           27           20           85 
##    Whiterock 
##           31
# convert soil colors to CIE LAB colorspace
# moist colors
cols.lab <- as(sRGB(x$m_r, x$m_g, x$m_b), 'LAB')@coords
x$m_L <- cols.lab[, 1]
x$m_A <- cols.lab[, 2]
x$m_B <- cols.lab[, 3]

# dry colors
cols.lab <- as(sRGB(x$d_r, x$d_g, x$d_b), 'LAB')@coords
x$d_L <- cols.lab[, 1]
x$d_A <- cols.lab[, 2]
x$d_B <- cols.lab[, 3]


# 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.4823269 0.5789309 0.6316626 0.6675592 0.6826665   0      1
## 2      d_r    Amador 0.4823269 0.5789309 0.6316626 0.6675592 0.6826665   1      2
## 3      d_r    Amador 0.4823269 0.5789309 0.6316626 0.6675592 0.6826665   2      3
## 4      d_r    Amador 0.4823269 0.5789309 0.6316626 0.6675592 0.6826665   3      4
## 5      d_r    Amador 0.4823692 0.5812506 0.6334787 0.6675894 0.6826666   4      5
## 6      d_r    Amador 0.4062669 0.5553961 0.6330297 0.6675890 0.6826666   5      6
##   contributing_fraction
## 1             0.6842105
## 2             0.6842105
## 3             0.6842105
## 4             0.6842105
## 5             0.6842105
## 6             0.6842105
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.6316626 0.5264226 0.3892964 57.77242 5.630916 21.63351 13.28053 5.482295
## 2    Amador   1      2 0.6316626 0.5264226 0.3892964 57.77242 5.630916 21.63351 13.28053 5.482295
## 3    Amador   2      3 0.6316626 0.5264226 0.3892964 57.77242 5.630916 21.63351 13.28053 5.482295
## 4    Amador   3      4 0.6316626 0.5264226 0.3892964 57.77242 5.630916 21.63351 13.28053 5.482295
## 5    Amador   4      5 0.6334787 0.5263167 0.3917715 57.77208 5.641995 21.10100 13.61498 5.502247
## 6    Amador   5      6 0.6330297 0.5255969 0.3909398 57.75578 5.709882 20.62976 13.63533 5.536855
##   total_frags_pct
## 1        6.647779
## 2        6.647779
## 3        6.647779
## 4        6.647779
## 5        6.647779
## 6        7.401305
# 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(as(LAB(x.colors$d_L, x.colors$d_A, x.colors$d_B), 'sRGB')@coords)

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

# aggregate bedrock depth probabilty by taxonname
# at 90% level of confidence
dp <- aggregateSoilDepth(x, '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', 'Hennekenot'), 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,0))
plot(x.colors, divide.hz=FALSE, name='', plot.order=new.order, col.label='Soil Color', lwd=1.25, axis.line.offset=-6, 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, diana(d), scaling.factor = 1, name='',  divide.hz=FALSE, y.offset=10)
addBracket(depth.brackets, col='black', label.cex=0.6)
title('Aggregate Soil Properties')

par(mar=c(1,0,3,0))
plot(x.colors, divide.hz=FALSE, color='clay', name='', plot.order=new.order, col.label='Clay Content (%)', lwd=1.25, axis.line.offset=-6, 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='', plot.order=new.order, col.label='pH', lwd=1.25, axis.line.offset=-6, 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='', plot.order=new.order, col.label='Total Rock Fragment Volume (%)', lwd=1.25, axis.line.offset=-6, 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 = 0.5, name='',  divide.hz=FALSE, y.offset=10)

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