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)

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

head(x.colors)

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

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

plotProfileDendrogram(x.colors, diana(d), scaling.factor = 0.5, name='',  divide.hz=FALSE, y.offset=10, 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 = 0.5, 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, offset = -0.25)

pig <- soilColorSignature(x.colors, r='d_r', g='d_g', b='d_b', RescaleLightnessBy = 1, method = 'colorBucket')
# move row names over for distance matrix
row.names(pig) <- pig$taxonname

head(pig)

d <- daisy(pig[, -1])
dd <- diana(d)

par(mar=c(1,1,1,1))
plotProfileDendrogram(x.colors, dd, dend.y.scale = max(d) * 2, scaling.factor = 0.002, y.offset = 0.04, width=0.15, name='',  divide.hz=FALSE)
addBracket(x.colors$soil.top, x.colors$soil.bottom, col='black', label='P(soil >= 90%)', label.cex=0.6, offset = -0.25)

pig <- soilColorSignature(x.colors, r='d_r', g='d_g', b='d_b', RescaleLightnessBy = 1, method = 'depthSlices')
# move row names over for distance matrix
row.names(pig) <- pig$taxonname

d <- daisy(pig[, -1])
dd <- diana(d)

par(mar=c(1,1,1,1))
plotProfileDendrogram(x.colors, dd, dend.y.scale = max(d) * 2, scaling.factor = 0.5, y.offset = 7, width=0.15, name='',  divide.hz=FALSE)
addBracket(x.colors$soil.top, x.colors$soil.bottom, col='black', label='P(soil >= 90%)', label.cex=0.6, offset = -0.25)

pig <- soilColorSignature(x.colors, r='d_r', g='d_g', b='d_b', RescaleLightnessBy = 1, method = 'pam')
# move row names over for distance matrix
row.names(pig) <- pig$taxonname

d <- daisy(pig[, -1])
dd <- diana(d)

par(mar=c(1,1,1,1))
plotProfileDendrogram(x.colors, dd, dend.y.scale = max(d) * 2, scaling.factor = 0.5, y.offset = 7, width=0.15, name='',  divide.hz=FALSE)
addBracket(x.colors$soil.top, x.colors$soil.bottom, col='black', label='P(soil >= 90%)', label.cex=0.6, offset = -0.25)