Introduction

This is a simple demonstration of how to do some soils-related stuff in R with the aqp package.

A Visual Demonstration of Horizon Depths Measured Vertically vs. Perpendicular to the Surface Slope

The relationship between perpendicular horizon depth (\(T_{p}\)) and vertical horizon depth (\(T_{v}\)) is defined as:

\[T_{p} = T_{v} cos(\theta)\]

where \(\theta\) is the surface slope measured in radians. Surface slope measured in percent (\(s\)) can be converted to radians with \(\theta = atan(s / 100)\).

Setup environment by loading required packages, and user-defined functions

# load required packages, you may have to install these if missing:
# install.packages('Hmisc', dep=TRUE)
library(lattice)
library(Hmisc)
library(Cairo)
library(aqp)
library(xtable)

# convert angle in percent to radians
pct2rad <- function(pct) {
    atan(pct / 100)
}

# compute difference in thickness as a function of slope pct
thickDiff <- function(pct, Tv) {
    Tp <- Tv * cos(pct2rad(pct))
    return(Tv - Tp)
}

# transform vertical to perpendicular depths
# operates on a single SPC object
VtoP <- function(i, which) {
    # get vector of bottom or top depths
    b <- horizons(i)[[which]]
    # convert vertical to perpendicular depths and round to an integer
    b.p <- round(b * cos(pct2rad(i$slope)))
    return(b.p)
}

Generate a sample data set for subsequent demonstrations

# generate a single profile, keep only the relevant stuff
set.seed(54321)
p <- random_profile(1)[, c('id', 'top', 'bottom', 'name')]

# stack 10 copies of these data
pp <- rbind(p, p, p, p, p, p, p, p, p, p)

# add a unique ID to each chunk copy
pp$id <- rep(1:10, each=nrow(p))

# upgrade our data.frame into a SoilProfileCollection object
depths(pp) <- id ~ top + bottom

# generate a sequence of (10) simulated slope values
s <- data.frame(id=1:10, slope=seq(from=1, to=70, length.out=10))
# assign these values into the 'site' attribute table of our pedons
site(pp) <- s

# assume data were generated with vertical horizon depths
# compute corresponding perpendicular horizon bottom depths
pp$Tp <- profileApply(pp, VtoP, which='bottom')

# check our new object:
print(pp)
## Object of class SoilProfileCollection
## Number of profiles: 10
## Depth range: 48-48 cm
## 
## Horizon attributes:
##   id top bottom name Tp
## 1  1   0     17   H1 17
## 2  1  17     26   H2 26
## 3  1  26     38   H3 38
## 4  1  38     48   H4 48
## 5  2   0     17   H1 17
## 6  2  17     26   H2 26
## 
## Sampling site attributes:
##   id     slope
## 1  1  1.000000
## 2  2  8.666667
## 3  3 16.333333
## 4  4 24.000000
## 5  5 31.666667
## 6  6 39.333333

Plot differences in horizon depth as a function of slope

# adjust figure margins
par(mar=c(2,0,0.75,0))
# plot SoilProfileCollection
plot(pp, n.depth.ticks=10)

# add a title
title('Vertical Horizon Depths (black) vs. Perpendicular Horizon Depths (blue)', cex.main=0.75)

# add axis with slope values
axis(side=1, cex.axis=0.75, line=-1.25, at=1:10, labels=round(s$slope))
mtext(side=1, 'Slope (%)', line=0.5, cex=0.75)

# iterate over profiles and add offset hz boundaries
for(i in 1:length(pp)) {
    pp.i <- pp[i, ]
    # use color to show only those offsets that are > 1 cm
    cols <- ifelse(abs(pp.i$bottom - pp.i$Tp) > 1, 'RoyalBlue', NA)
    arrows(x0=i, x1=i, y0=pp.i$bottom, y1=pp.i$Tp, length=0.08, col=cols)
    segments(x0=i-0.2, x1=i+0.2, y0=pp.i$Tp, y1=pp.i$Tp, col=cols, lty=2)
}

plot of chunk plot-pedons-with-arrows

This time plot after converting horizon depths from vertical basis to perpendicular basis

# make a copy of the original profiles
pp.perp <- pp

# convert horizon top and bottom depths from vertical to perpendicular basis
pp.perp$top <- profileApply(pp.perp, VtoP, which='top')
pp.perp$bottom <- profileApply(pp.perp, VtoP, which='bottom')

# setup margins
par(mar=c(2,0,0.75,0))
# plot profiles
plot(pp.perp, n.depth.ticks=10)
# add title
title('Changes in Horizon Depth as a Function of Surface Slope', cex.main=0.75)
# add slope axis
axis(side=1, cex.axis=0.75, line=-1.25, at=1:10, labels=round(s$slope))
mtext(side=1, 'Slope (%)', line=0.5, cex=0.75)
# add horizontal lines to assist with perception of depth
abline(h=seq(from=5, to=75, by=5), lty=3, col='grey')

plot of chunk change-depths

Tabulate the mean and cumulative difference in horizon depths by pedon

# compute difference in horizon thickness
diff.thick <- (pp$bottom - pp$top) - (pp.perp$bottom - pp.perp$top)
theID <- horizons(pp)$id # extract IDs
# compute mean and total difference by pedon
diff.by.pedon <- by(diff.thick, theID, function(i) cbind(mean.diff.cm=mean(i), total.diff.cm=sum(i)))
# composite into a data.frame and add slope values
diff.by.pedon.df <- data.frame(do.call('rbind', diff.by.pedon), slope=round(s$slope))
# convert tabular output into HTML
print(xtable(diff.by.pedon.df), type='html')
mean.diff.cm total.diff.cm slope
1 0.00 0.00 1.00
2 0.00 0.00 9.00
3 0.25 1.00 16.00
4 0.25 1.00 24.00
5 0.50 2.00 32.00
6 0.75 3.00 39.00
7 1.25 5.00 47.00
8 1.50 6.00 55.00
9 1.75 7.00 62.00
10 2.25 9.00 70.00

Another visualization of differences in vertical vs. perpendicular depths as a function of slope

# make a grid of horizon thickness and slope percents
d <- expand.grid(thick=c(5,10,20,50,100), slope=1:70)

# compute the difference between vertical and perpendicular depths
d$thickDiff <- thickDiff(d$slope, d$thick)

# check the results
head(d)
##   thick slope    thickDiff
## 1     5     1 0.0002499813
## 2    10     1 0.0004999625
## 3    20     1 0.0009999250
## 4    50     1 0.0024998125
## 5   100     1 0.0049996250
## 6     5     2 0.0009997001
# plot these data
# note the different syntax-- this is lattice-style plotting
xYplot(thickDiff ~ slope, data=d, groups=factor(thick, labels=paste(c(5,10,20,50,100), 'cm depth')), type='l',
ylab='(Vertical - Perpendicular) Hz Thickness', xlab='Slope (%)', 
main='Difference in Apparent Horizon Thickness as a Function of Slope', 
scales=list(alternating=3, tick.number=10), 
par.settings=list(superpose.line=list(col='black', lwd=1.5)), 
label.curves=list(cex=1), 
panel=function(...) {
    panel.abline(v=seq(0, 70, by=10), h=seq(0, 20, by=2), col=grey(0.75), lty=3)
    panel.xYplot(...)
    })

plot of chunk plot-diff-by-slope