This document is based on aqp version 1.10-1 and soilDB version 1.8.1.

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(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')
# 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 Whiterock
19 20 112 22 27 20 85 31
# convert soil colors to 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 contributing_fraction
d_r Amador 0.4823269 0.5789309 0.6316626 0.6675592 0.6826665 0 1 0.6842105
d_r Amador 0.4823269 0.5789309 0.6316626 0.6675592 0.6826665 1 2 0.6842105
d_r Amador 0.4823269 0.5789309 0.6316626 0.6675592 0.6826665 2 3 0.6842105
d_r Amador 0.4823269 0.5789309 0.6316626 0.6675592 0.6826665 3 4 0.6842105
d_r Amador 0.4823692 0.5812506 0.6334787 0.6675894 0.6826666 4 5 0.6842105
d_r Amador 0.4062669 0.5553961 0.6330297 0.6675890 0.6826666 5 6 0.6842105
head(x.colors)
taxonname top bottom d_r d_g d_b d_L d_A d_B clay phfield total_frags_pct
Amador 0 1 0.6316626 0.5264226 0.3892964 57.77242 5.630916 21.63351 13.28053 5.482295 6.647779
Amador 1 2 0.6316626 0.5264226 0.3892964 57.77242 5.630916 21.63351 13.28053 5.482295 6.647779
Amador 2 3 0.6316626 0.5264226 0.3892964 57.77242 5.630916 21.63351 13.28053 5.482295 6.647779
Amador 3 4 0.6316626 0.5264226 0.3892964 57.77242 5.630916 21.63351 13.28053 5.482295 6.647779
Amador 4 5 0.6334787 0.5263167 0.3917715 57.77208 5.641995 21.10100 13.61498 5.502247 6.647779
Amador 5 6 0.6330297 0.5255969 0.3909398 57.75578 5.709881 20.62976 13.63533 5.536855 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))
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(x.colors$soil.top, x.colors$soil.bottom, col='black', label='P(soil >= 90%)', 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(x.colors$soil.top, x.colors$soil.bottom, col='black', label='P(soil >= 90%)', 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(x.colors$soil.top, x.colors$soil.bottom, col='black', label='P(soil >= 90%)', 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(x.colors$soil.top, x.colors$soil.bottom, col='black', label='P(soil >= 90%)', 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(x.colors$soil.top, x.colors$soil.bottom, col='black', label='P(soil >= 90%)', 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(x.colors$soil.top, x.colors$soil.bottom, col='black', label='P(soil >= 90%)', label.cex=0.6)
addVolumeFraction(x.colors, 'total_frags_pct')