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)