1 Introduction

The semi-automatic creation of soil profile sketches from morphologic data (OSDs, SSURGO components, NASIS pedons, etc.) has become something of an obsession for me, starting with field notes collected during my first season mapping soils in Pinnacles National Park, CA (CA788). These sketches appeared in SoilWeb circa 2007, with approximate colors derived from the OSDs in 2009. The {aqp} package (2011) for R formalized one possible process for converting basic horizon morphology into cartoon-like sketches of soil profiles. If interested, I’d highly recommend the SoilProfileCollection object documentation for a detailed description of soil profile sketch authoring.

1.1 Horizon Boundaries

1.2 Visual Encoding Ideas

1.2.1 Distinctness

  • diagonal lines depict the average depth over which the boundary is defined

    • do we need to use a more rigorous geometry model?
    • may be hard to interpret without detailed explanation / examples
    • visually compact
    • simple geometry
    • readily constrained by adjacent horizon geometry

  • color blending or shaded regions of overlap

  • horizon boundary line thickness

1.2.2 Topography

1.3 Required Packages

# stable version + dependencies
install.packages(c('aqp', 'soilDB'), dependencies = TRUE)

# development version
remotes::install_github("ncss-tech/aqp", dependencies=FALSE, upgrade=FALSE, build=FALSE)

1.4 Setup

# required libraries
library(aqp)
library(soilDB)

1.5 Offset Encoding

# same result
hzDistinctnessCodeToOffset('G')
hzDistinctnessCodeToOffset('gradual')

1.6 Get Some Data

# select soil series names
soil.series <- c('leon', 'musick', 'clarksville', 'pardee', 'amador', 'lucy', 'dylan', 'tristan', 'pierre', 'drummer', 'zook')

# get the data
osds <- fetchOSD(soil.series)

# encode horizon boundary distinctness via vertical offset
osds$hd <- hzDistinctnessCodeToOffset(osds$distinctness)

# encode horizon boundary topography via vertical offset
osds$hzto <- hzTopographyCodeToOffset(osds$topography)

# also encode horizon boundary topography as line type
osds$hzto.lty <- hzTopographyCodeToLineType(osds$topography)

# label data source, used later 
site(osds)$source <- 'OSD'

# concise representation of 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)

2 Horizon Sketches in {aqp}

A typical set of sketches, derived from the OSDs, is presented below. Note that horizon boundaries are depicted as horizontal lines—suggesting very abrupt/smooth (VS) horizon boundaries. We all know that VS horizon boundaries are not that common; can we do better? Examples used in this article have been sourced via fetchOSD() of the {soilDB} R package.

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

plotSPC(osds, width = 0.3, cex.id = 0.66, cex.names = 0.55)

Using the horizon boundary data from the OSDs and encoding via “diagonals” (distinctness) and “chevrons” (topography), we see the following. Note that horizon boundary topography has also been encoded via line type for clarity. What do you think? Is this an improvement or just clutter? Feel free to contact me with questions, suggestions, or a completely different approach to increasing the information density of these sketches. The code used to generate these figures is posted on the {aqp} GitHub site.

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

plotSPC(osds, width = 0.3, hz.distinctness.offset = 'hd', hz.topography.offset = 'hzto', cex.id = 0.66, cex.names = 0.55, hz.boundary.lty = 'hzto.lty') 

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

3 Encoding Horizon Distinctness and Topography

Horizon boundary distinctness describes the vertical distance over which one horizon transitions to the next. Horizon boundary topography describes the complexity of the horizon boundary. In {aqp} these two features are encoded using diagonal lines (distinctness) and “chevrons” (topography). The following figure demonstrates several possible combinations of horizon boundary distinctness and topography. Horizon depths are marked with small green symbols. The magnitude of these effects can be controlled when converting text boundary codes (e.g. CW) into representative “offsets”. The default values have been selected based on an interpretation of our Field Book for Describing and Sampling Soils. The {aqp} manual contains a detailed explanation of usage and assumptions.

# a single series
x <- subset(osds, id == 'DRUMMER')

# # full set of boundary topography codes
# g <- expand.grid(
#   distinctness = c('V', 'A', 'C', 'G', 'D'),
#   topography = c('S', 'W', 'I', 'B')
# )

# simplified set, for this example
g <- expand.grid(
  distinctness = c('A', 'C', 'G', 'D'),
  topography = c('S', 'W', 'I')
)


# storage for intermediate results
l <- list()
# copies of the original OSD
s <- duplicate(x, times = nrow(g))

# iterate over all possible combinations
# creating / modifying a copy of the OSD at each iteration
for(i in 1:nrow(g)) {
  # current profile
  ss <- s[i, ]
  # set for all horizons
  horizons(ss)$hd <- hzDistinctnessCodeToOffset(g$distinctness[i])
  horizons(ss)$ht <- hzTopographyCodeToOffset(g$topography[i])
  horizons(ss)$ht.lty <- hzTopographyCodeToLineType(g$topography[i]) 
  # save modified copy
  l[[i]] <- ss  
}

# combine list elements -> SPC
s <- combine(l)

# associated horizon boundary codes
s$source <- sprintf("%s%s", g$distinctness, g$topography)

# make a figure
par(mar = c(1, 0, 2, 2), bg = 'black', fg = 'white')
plotSPC(s, width = 0.33, print.id=TRUE, hz.distinctness.offset = 'hd', hz.topography.offset = 'ht', label='source', cex.names=0.8, name = NA, color = NA, default.color = 'black', hz.boundary.lty = 'ht.lty', max.depth = 185)
title('Horizon Boundary Types', line = 0, col.main = 'white')

# label original horizon depths
p.seq <- rep(1:length(s), each = nrow(s[1, ]))
points(x = p.seq, y = s$bottom, pch = 15, col = 'green', cex = 0.5)

# line type legend
legend('bottom', horiz = TRUE, legend = c('Smooth', 'Wavy', 'Irregular', 'Broken'), lty = 1:4, inset = 0)

3.1 Examples

Horizon boundary distinctness demonstration for the LEON soil series, data sourced from the OSD. Note that vertical offset of the “diagonal” boundaries is constrained by the depths of adjacent horizons.

# keep a single OSD for the demo
x <- subset(osds, id == 'LEON')

# safely duplicate these data 6 times
# --> new IDs generated, old IDs saved into .oldID
s <- duplicate(x, times = 6)

# storage for intermediate SPC objs
l <- list()

# iterate over horizon distinctness depths
hzdo <- c(0, 0.5, 2, 5, 15, 20)/2

for(i in seq_along(hzdo)) {
  # current profile
  ss <- s[i, ]
  # set for all horizons
  ss$hd <- hzdo[i]
  l[[i]] <- ss
}

# combine copies of the OSD, with hz distinctness examples
s <- combine(l)

# horizon distinctness codes
s$source <- c('none', 'very abrupt', 'abrupt', 'clear', 'gradual', 'diffuse')

# combine OSD with horizon distinctness demo profiles
z <- combine(x, s)

# plot
par(mar = c(0, 0, 2, 2), bg = 'black', fg = 'white')
plotSPC(z, width=0.33, print.id=TRUE, hz.distinctness.offset = 'hd', label='source', cex.names=0.75)
title('Horizon Boundary Distinctness', line = 0, col.main = 'white')

Horizon boundary topography demonstration for the LEON soil series, data sourced from the OSD. Note that vertical offset of the “chevron” is constrained by the depths of adjacent horizons.

# safely duplicate these data 4 times
# --> new IDs generated, old IDs saved into .oldID
s <- duplicate(x, times = 4)

# storage for intermediate SPC objs
l <- list()

# iterate over horizon topography offsets
# theses aren't the default values but maybe should be
hzto <- c(0, 5, 10, 25)

for(i in seq_along(hzto)) {
  # current profile
  ss <- s[i, ]
  # set for all horizons
  ss$hzto <- hzto[i]
  l[[i]] <- ss
}

# combine copies of the OSD, with hz topography examples
s <- combine(l)

# horizon distinctness codes
s$source <- c('smooth', 'wavy', 'irregular', 'broken')


# combine OSD with horizon topography demo profiles
z <- combine(x, s)

par(mar = c(0, 0, 2, 2), bg = 'black', fg = 'white')
plotSPC(z, width=0.25, print.id=TRUE, hz.topography.offset = 'hzto', label='source', cex.names=0.75)
title('Horizon Boundary Topography', line = 0, col.main = 'white')

4 Alternative Approaches

Horizon distinctness via vertical segments with transparency.

# TODO: generalize this
hdn <- horizonDepths(osds)

# extract relevant pieces: profile ID, top, bottom, horizon distinctness information
s <- horizons(osds)[, c(idname(osds), hdn, 'hd')]

# TODO: truncate at adjacent horizon depths
# bracket bounds
s$top <- s$bottom - s$hd
s$bottom <- s$bottom + s$hd

par(mar = c(0, 0, 0, 1.5))
plotSPC(osds, width = 0.25, cex.id = 0.66, cex.names = 0.55, hz.distinctness.offset = 'hd', max.depth = 200) 

# encode with brackets
addBracket(s, tick.length = 0, offset = -0.33, lwd = 5, col=rgb(red=0, green=0, blue=1, alpha=0.5))

par(mar = c(0, 0, 0, 0))
plotSPC(osds, width = 0.25, cex.id = 0.66, cex.names = 0.55, name.style = 'center-center', hz.depths = TRUE, fixLabelCollisions = TRUE, hz.depths.offset = 0.08, depth.axis = FALSE, hz.distinctness.offset = 'hd', max.depth = 200) 

addBracket(s, tick.length = 0, offset = -0.33, lwd = 5, col=rgb(red=0, green=0, blue=1, alpha=0.5))

This document is based on aqp version 2.1.0, and soilDB version 2.8.5.