This tutorial outlines a simple method for summarizing soil color from a collection of soil profiles using the
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
aggregateColorPlot contain additional details.
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.
library(soilDB) library(aqp) library(sharpshootR) library(cluster)
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 par(mar=c(0,0,0,4)) plotSPC(loafercreek, name = NA, print.id = 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.
# 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) par(mar=c(4.5,2.5,4.5,0)) aggregateColorPlot(aggregateColor(s, 'slice'), label.cex=0.65, main='Loafercreek Dry Colors\nDepth Slices', print.n.hz = TRUE)
# 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)