This tutorial outlines a simple method for summarizing soil color from a collection of soil profiles using the aqp and sharpshootR packages for R. Suppose you have a collection of soil profiles and would like to compare the relative frequencies of soil color among different horizons, depths, hillslope positions, or bedrock. How would you determine the most frequent color or a reasonable range in colors for a given group? One solution to this problem is presented in the figure below. In this case, relative proportions of soil color are arranged by genetic horizon. Larger “color bars” are the most frequently described colors, both in terms of number of observed horizons and the thickness of those horizons. The small numbers in parenthesis denote the number of horizons associated with any given color. See examples below for more ideas. The manual page for aggregateColor and aggregateColorPlot contain additional details.

Setup R Envionment

If you have never used the aqp or sharpshootR packages before, you will likely need to install them. This only needs to be done once.

# stable version from CRAN and all dependencies
install.packages('aqp', dep=TRUE) 
install.packages('soilDB', dep=TRUE)
install.packages('sharpshootR', dep=TRUE)

Now that you have all of the R packages that this document depends on, it would be a good idea to load them. R packages must be installed anytime you change versions of R (e.g. after an upgrade), and loaded anytime you want to access functions from within those packages.


Sample Data

While the methods outlined in this document can be applied to any collection of pedons, it is convenient to work with a standardized set of data. You can follow along with the analysis by copying code from the following blocks and running it in your R session. The sample data used in this document is based on 30 soil profiles that have been correlated to the Loafercreek soil series from the Sierra Nevada Foothill Region of California. Note that the internal structure of the loafercreek data is identical to the structure returned by fetchNASIS() from the soilDB package. All horizon-level values are pulled from the pedon horizon table of the pedons being analyzed.

# load sample data from the soilDB package
data(loafercreek, package = 'soilDB')
# graphical check
plotSPC(loafercreek, name = NA, = FALSE, cex.names = 0.8, axis.line.offset = -1, max.depth = 100, divide.hz = FALSE)


These examples can be readily adapted to your own data. Note that you may need to adjust object (change loafercreek to the name of your SoilProfileCollection object) and column (genhz, etc.) names. Data loaded from NASIS via fetchNASIS can be used without modification (except for object name) in these examples.

Typically color summaries are grouped by some kind of “horizon-level” attribute: depth slice, horizon designation, generalized horizon label, or within a diagnostic feature. However, it is possible to group colors by a “site-level” attribute such as bedrock kind, taxonname, or hillslope position. Note that grouping colors by a site-level attribute will pool colors from all depths.

Depth Slices

# slice color data at select depths
s <- slice(loafercreek, c(5, 10, 15, 25, 50, 75) ~ soil_color, strict = FALSE)

# make horizon labels based on slice depth
s$slice <- paste0(s$hzdept, ' cm')
s$slice <- factor(s$slice, levels=guessGenHzLevels(s, 'slice')$levels)

aggregateColorPlot(aggregateColor(s, 'slice'), label.cex=0.65, main='Loafercreek Dry Colors\nDepth Slices', print.n.hz = TRUE)

Generalized Horizon Labels

# generalize horizon names using REGEX rules
n <- c('Oi', 'A', 'BA','Bt1','Bt2','Bt3','Cr','R')
p <- c('O', '^A$|Ad|Ap|AB','BA$|Bw', 'Bt1$|^B$','^Bt$|^Bt2$','^Bt3|^Bt4|CBt$|BCt$|2Bt|2CB$|^C$','Cr','R')
loafercreek$genhz <- generalize.hz(loafercreek$hzname, n, p)

# remove non-matching generalized horizon names
loafercreek$genhz[loafercreek$genhz == 'not-used'] <- NA
loafercreek$genhz <- factor(loafercreek$genhz)

aggregateColorPlot(aggregateColor(loafercreek, 'genhz'), main='Loafercreek Series Dry Colors\nGeneralized Horizon Labels', print.n.hz = TRUE, label.cex = 0.65)