library(aqp)
library(soilDB)
library(sharpshootR)
library(dendextend)

# soils of interest
s.list <- c('hornitos', 'perkins', 'argonaut', 'inks', 'mokelumne', 'dunstone', 'auburn', 'pentz', 'pardee', 'peters', 'amador', 'laniger')

# fetch and convert data into an SPC
s <- fetchOSD(s.list)
# estimate soil depth based on horizon designation
s$soil.depth <- profileApply(s, estimateSoilDepth, name = 'hzname')

# plot profiles, order by soil depth
par(mar = c(0, 0, 0, 0))
plotSPC(s, name = 'hzname', plot.order = order(s$soil.depth), cex.names = 0.85, axis.line.offset = -5, width = 0.33, id.style = 'top', name.style = 'center-center', plot.depth.axis = FALSE, hz.depths = TRUE)

Divisive clustering, the default method.

# plot dendrogram + profiles
# result is the distance matrix used to generate dendrogram
d <- SoilTaxonomyDendrogram(s, cex.taxon.labels = 0.8, width = 0.33, y.offset = 0.4)

Try agglomerative clustering, need to adjust offset.

# plot dendrogram + profiles
SoilTaxonomyDendrogram(s, cex.taxon.labels = 0.8, width = 0.33, y.offset = 0.5, cluster.method = 'agglomerative')

Approximate sorting of dendrogram according to catenary position (hillslope position proportions).

# TODO: try norforlk and siblings
# http://soilmap2-1.lawr.ucdavis.edu/dylan/soilweb/sde/?series=norfolk

# set of soil series to query
soils <- c('cecil', 'altavista', 'lloyd', 'wickham', 'wilkes',  'chewacla', 'congaree')

# get morphology + extended summaries
s <- fetchOSD(soils, extended = TRUE)

# extract catenary relationships from hillslope proportion summary
# and display
res <- vizHillslopePosition(s$hillpos)
print(res$fig)

Organize ST dendrogram according to catenary position.

par(mar=c(0,0,1,1), mfrow=c(2,1))

SoilTaxonomyDendrogram(s$SPC, width = 0.33, cex.taxon.labels = 0.88, shrink = TRUE)
mtext('default sorting', side = 2, line=-1, font=3, cex=1.25)

# labels used to re-order (if possible) terminal branches of dendrogram
labs <- profile_id(s$SPC)[res$order]

SoilTaxonomyDendrogram(s$SPC, rotationOrder = labs, width = 0.33, cex.taxon.labels = 0.88, shrink = TRUE)
mtext('approx. catenary sorting', side = 2, line=-1, font=3, cex=1.25)

A classic chronosequence from the San Joaquin Valley, CA (Jenny, 1981).

# series of interest, in order from youngest -> oldest
s <- c('tujunga', 'hanford', 'greenfield', 'snelling', 'san joaquin')
osds <- fetchOSD(s)

# ordering vector for sketches
idx <- match(toupper(s), profile_id(osds))

# encode horizon boundary distinctness via vertical offset
osds$hd <- hzDistinctnessCodeToOffset(
  osds$distinctness,
  codes=c('very abrupt', 'abrupt', 'clear', 'gradual', 'diffuse')
)

# encode horizon boundary topography via vertical offset
osds$hzto <- hzTopographyCodeToOffset(
  osds$topography,
  codes = c('smooth', 'wavy', 'irregular', 'broken')
)

# also encode horizon boundary topography las line type
osds$hzto.lty <- hzTopographyCodeToLineType(
  osds$topography,
  codes = c('smooth', 'wavy', 'irregular', 'broken')
)


# concise representation of hz bnd distinctness and topography
# similar to field notes
osds$bnd.code <- sprintf(
  "%s%s",
  substr(osds$distinctness, 1, 1),
  substr(osds$topography, 1, 1)
)

# remove NA
osds$bnd.code <- gsub('NANA', '', osds$bnd.code)

par(mar = c(0, 0, 0, 1), bg = 'black', fg = 'white')

plotSPC(
osds,
plot.order = idx,
width = 0.3,
name.style = 'center-center',
cex.names = 0.66,
plot.depth.axis = FALSE,
hz.depths = TRUE,
fixLabelCollisions = TRUE,
shrink = TRUE,
hz.distinctness.offset = 'hd',
hz.topography.offset = 'hzto',
hz.boundary.lty = 'hzto.lty'
)

legend(
'bottomright',
horiz = TRUE,
legend = c('Smooth', 'Wavy', 'Irregular', 'Broken'),
lty = 1:4,
inset = 0.05,
bty = 'n',
cex = 0.85
)

# note that `rotationOrder` uses the ordering of series names (uppercase to match profile IDs)
# to re-order the terminal branches of the dendrogram
SoilTaxonomyDendrogram(
osds,
rotationOrder = toupper(s),
cex.taxon.labels = 0.85,
width = 0.3,
name.style = 'center-center',
cex.names = 0.66,
plot.depth.axis = FALSE,
hz.depths = TRUE,
fixLabelCollisions = TRUE,
shrink = TRUE,
hz.distinctness.offset = 'hd',
hz.topography.offset = 'hzto',
hz.boundary.lty = 'hzto.lty'
)

legend('bottomright', horiz = TRUE, legend = c('Smooth', 'Wavy', 'Irregular', 'Broken'), lty = 1:4, inset = 0.01, bty = 'n', cex = 0.85)

This document is based on aqp version 2.0, soilDB version 2.7.8, and sharpshootR version 2.2.