Setup R session.

# latest version from github
library(aqp)
library(soilDB)
library(sharpshootR)

library(cluster)
library(MASS)
library(lattice)
library(tactile)

library(vegan)

Get some example data via OSDs and sketch using horizon boundary distinctness.

# example profile
x <- fetchOSD(c('musick', 'drummer', 'sierra'))

# sketch using horizon boundary encoded
par(mar = c(0, 0, 0, 0))
plotSPC(x, hz.distinctness.offset = 'hzd', name.style = 'center-center', plot.depth.axis = FALSE, hz.depths = TRUE, cex.names = 0.75)

Emphasize horizon boundary codes / translation to vertical variation.

x$hz.bdy.txt <- sprintf("%s: %scm", x$distinctness, x$hzd)

par(mar = c(0, 0, 0, 0))
plotSPC(x, hz.distinctness.offset = 'hzd', name.style = 'center-center', plot.depth.axis = FALSE, hz.depths = TRUE, name = 'hz.bdy.txt', cex.names = 0.75)

Simulation from a template: perturbation of horizon boundaries based on horizon boundary distinctness.

# number of new IDs sets the number of realizations
s <- perturb(
  subset(x, id == 'MUSICK'), 
  id = sprintf("Sim. %02d", 1:9),
  boundary.attr = 'hzd', 
  min.thickness = 5
)

# combine with original
o <- subset(x, id == 'MUSICK')

# combine multiple SPCs -> single SPC
s <- combine(
  o,
  s
)

# neat
par(mar=c(0,0,0,0))
plotSPC(s, name.style = 'center-center', plot.depth.axis = FALSE, hz.depths = TRUE, width = 0.3, cex.names = 0.75)

Align many realizations and suppress all labels.

# simulate
s <- perturb(
  subset(x, id == 'MUSICK'), 
  n = 50, 
  boundary.attr = 'hzd', 
  min.thickness = 5
)

# truncate all profiles to 0-150cm
s <- trunc(s, 0, 150)

# no foreground color (normally white)
# looks like a cake
par(mar=c(0,0,0,0), fg = NA)
plotSPC(s, print.id = FALSE, name = NA, divide.hz = FALSE, plot.depth.axis = FALSE, width=0.55, x.idx.offset = 2, y.offset = -15)

s <- perturb(
  subset(x, id == 'DRUMMER'), 
  n = 50, 
  boundary.attr = 'hzd', 
  min.thickness = 5
)

# truncate all profiles to 0-150cm
s <- trunc(s, 0, 150)

par(mar=c(0,0,0,0), fg=NA)
plotSPC(s, print.id=FALSE, name=NA, divide.hz = FALSE, plot.depth.axis = FALSE, width=0.55, x.idx.offset = 2, y.offset = -15)

Sort the profile sketches based on a comparison of generalized horizon designation and how it varies with depth.

musick <- subset(x, id == 'MUSICK')
horizons(musick)$hd <- 5

# simulate
s <- perturb(
  musick, 
  id = sprintf("%03d", 1:50),
  boundary.attr = 'hd', 
  min.thickness = 5
)

# truncate all profiles to 0-150cm
s <- trunc(s, 0, 150)

hzn.id <- hzdesgnname(s)
hzn <- hzDesgn(s[1, ])
s$genhz <- generalize.hz(s[[hzn.id]], new = hzn, pat = hzn)

d <- NCSP(s, vars = c('genhz'), maxDepth = 150, k = 0)
new.order <- as.hclust(diana(d))$order

# looks like a cake
par(mar=c(0,0,0,0), fg=NA)
plotSPC(s, print.id = FALSE, name = NA, divide.hz = FALSE, plot.depth.axis = FALSE, width=0.55, x.idx.offset = 2, plot.order = new.order, y.offset = -15)

Use simulation to demonstration pair-wise distances between profiles. A “spike” is added so that relative distance can be interpreted.

# example data
data(sp4)
depths(sp4) <- id ~ top + bottom

# simulated based on profile 1
p.idx <- 1

# use profile 6 as a "spike"
spike.idx <- 6

# simulate 
horizons(sp4)$bdy <- 4
set.seed(1010101)
p <- perturb(sp4[p.idx, ], n = 10, boundary.attr = 'bdy', min.thickness = 2)

# remove the old profile ID column
site(p)$id <- NULL

# combine: simulations + original + spike
z <- combine(
    p, 
    sp4[p.idx, ],
    sp4[spike.idx, ]
)

par(mar = c(0,0,3,1))
plotSPC(z, color = 'ex_Ca_to_Mg', width = 0.3, name.style = 'center-center', plot.depth.axis = TRUE, hz.depths = TRUE, cex.names = 0.75, col.label = 'Exchangeable Ca:Mg')

d <- NCSP(z, vars = c('ex_Ca_to_Mg', 'CEC_7'), maxDepth = 40, k = 0)
dd <- diana(d)

plotProfileDendrogram(z, dd, scaling.factor = 1.5, y.offset = 5, color = 'ex_Ca_to_Mg', width = 0.3, name.style = 'center-center', plot.depth.axis = FALSE, hz.depths = TRUE, cex.names = 0.75, col.label = 'Exchangeable Ca:Mg')

# simulate three of each soil
b <- perturb(sp4, n = 3, boundary.attr = 'bdy', min.thickness = 2)

# compare
d <- NCSP(b, vars = c('ex_Ca_to_Mg', 'CEC_7'), maxDepth = 50, k = 0)
dd <- diana(d)

# viz
par(mar = c(1, 1, 0, 0))
plotProfileDendrogram(b, dd, scaling.factor = 1.5, y.offset = 4, color = 'ex_Ca_to_Mg', width = 0.3, name.style = 'center-center', plot.depth.axis = FALSE, hz.depths = TRUE, cex.names = 0.65, col.label = 'Exchangeable Ca:Mg')

musick <- subset(x, id == 'MUSICK')
horizons(musick)$hd <- 5

drummer <- subset(x, id == 'DRUMMER')
horizons(drummer)$hd <- 5

sierra <- subset(x, id == 'SIERRA')
horizons(sierra)$hd <- 5


# simulate
s1 <- perturb(
  musick, 
  id = sprintf("Music-%03d", 1:50),
  boundary.attr = 'hd', 
  min.thickness = 5
)

# simulate
s2 <- perturb(
  drummer, 
  id = sprintf("Drummer-%03d", 1:50),
  boundary.attr = 'hd', 
  min.thickness = 5
)

# simulate
s3 <- perturb(
  sierra, 
  id = sprintf("Sierra-%03d", 1:50),
  boundary.attr = 'hd', 
  min.thickness = 5
)



s <- aqp::combine(s1, s2, s3)

# truncate all profiles to 0-150cm
s <- trunc(s, 0, 150)


s$genhz <- generalize.hz(s[[hzdesgnname(s)]], new = c('A', 'B', 'C'), pat = c('A', 'B', 'C'))
table(s$genhz)
## 
##        A        B        C not-used 
##      200      550      150       50
d <- NCSP(s, vars = 'genhz', maxDepth = 150, k = 0, isColor = FALSE)

## map distance matrix to 2D space via principal coordinates
d.betadisper <- betadisper(d, group = s$.oldID, bias.adjust = TRUE, sqrt.dist = FALSE, type='median')

## fancy plot
par(mar=c(3,3,3,1), mfcol=c(1,2))

boxplot(d.betadisper, varwidth=TRUE, las=1)

plot(
  d.betadisper, hull=FALSE, ellipse=TRUE, conf=0.5, las=1,
  col=c('Royalblue', 'Orange', 'Darkgreen', 'Firebrick'), 
  main='Ordination of Between-Profile Distances\n50% Probability Ellipse',
  xlab='', ylab=''
)

Simulation of plausible soil colors.

# simulateSoilColor()

This document is based on aqp version 2.0.